diff options
Diffstat (limited to 'vo/votools/gasplib')
| -rw-r--r-- | vo/votools/gasplib/calcds.x | 139 | ||||
| -rw-r--r-- | vo/votools/gasplib/ccgseq.x | 94 | ||||
| -rw-r--r-- | vo/votools/gasplib/ccgsxy.x | 200 | ||||
| -rw-r--r-- | vo/votools/gasplib/dcmpsv.f | 230 | ||||
| -rw-r--r-- | vo/votools/gasplib/eqtopix.x | 43 | ||||
| -rw-r--r-- | vo/votools/gasplib/fitsvd.f | 35 | ||||
| -rw-r--r-- | vo/votools/gasplib/gsctab.x | 66 | ||||
| -rw-r--r-- | vo/votools/gasplib/ksbsvd.f | 24 | ||||
| -rw-r--r-- | vo/votools/gasplib/mkpkg | 22 | ||||
| -rw-r--r-- | vo/votools/gasplib/pixtoeq.x | 37 | ||||
| -rw-r--r-- | vo/votools/gasplib/rdaslf.x | 131 | ||||
| -rw-r--r-- | vo/votools/gasplib/rdxy.x | 43 | ||||
| -rw-r--r-- | vo/votools/gasplib/regren.f | 309 | ||||
| -rw-r--r-- | vo/votools/gasplib/treqst.x | 33 | ||||
| -rw-r--r-- | vo/votools/gasplib/trsteq.x | 45 | ||||
| -rw-r--r-- | vo/votools/gasplib/varsvd.f | 21 |
16 files changed, 1472 insertions, 0 deletions
diff --git a/vo/votools/gasplib/calcds.x b/vo/votools/gasplib/calcds.x new file mode 100644 index 00000000..0f434b26 --- /dev/null +++ b/vo/votools/gasplib/calcds.x @@ -0,0 +1,139 @@ +include <math.h> +define SZ_CDMATRIX 4 # CD1_1, CD1_2, CD2_1, CD2_2 + +# CALCDS -- Procedure to calculate the values of the CD matrix from +# information given by a extracted Guide Star file. + +procedure calcds (plt_centre_ra, + plt_centre_dec, + plt_centre_x, + plt_centre_y, + x_corner, + y_corner, + x_pixel_size, + y_pixel_size, + plate_scale, + x_size, + y_size, + im_cen_ra, + im_cen_dec, + amd_x, + amd_y, + cd_matrix) + +double plt_centre_ra # Plate centre RA (radians) +double plt_centre_dec # Plate centre DEC +double plt_centre_x # X center position (microns) +double plt_centre_y # Y center position +int x_corner # x lower left of the extracted subset +int y_corner # y +double x_pixel_size # X scan pixel size (microns) +double y_pixel_size # Y scan pixel size +double plate_scale # Plate scale (arcsec/mm) +int x_size # Extracted image size x_axis (pixel) +int y_size # Extracted image size y_axis (pixel) +double im_cen_ra # Extracted image center RA (radians) +double im_cen_dec # Extracted image center DEC (radians) +double amd_x[ARB] # Xi plate solution coefficients +double amd_y[ARB] # Eta coefficients (arsec/mm) +double cd_matrix[SZ_CDMATRIX] # CD1_1, CD1_2, CD2_1, CD2_2 (degrees/pixel) + +pointer sp, xip, etap, x_arr, y_arr, u, v, w, cvm +int sx, sy, xlim, ylim +int i, j, k, nterms, xi, eta, npts +double x_coeff[2], y_coeff[2], xchisqr, ychisqr +double x_sigma[2], y_sigma[2], x, y, xc, yc +real ww[100] +double new_plt_centre_x, new_plt_centre_y, xref, yref, mag, color +double ra, dec +int nx,ny,nxy + +begin + mag = 0.0 + color= 0.0 + # calculate new plate center in microns + new_plt_centre_x = (x_size/2.)*x_pixel_size + new_plt_centre_y = (y_size/2.)*y_pixel_size + + call smark (sp) + call salloc (xip, 100, TY_DOUBLE) + call salloc (etap, 100, TY_DOUBLE) + call salloc (x_arr, 2*100, TY_DOUBLE) + call salloc (y_arr, 2*100, TY_DOUBLE) + + sx = max (1, x_size/10) + sy = max (1, y_size/10) + + xlim = x_size - mod(x_size, sx) + ylim = y_size - mod(y_size, sy) + nx = xlim/sx + ny = ylim/sy + nxy = nx*ny + xi = xip + eta = etap + k = 0 + do i = sx, xlim, sx { + y = i # x coord. from lower left + do j = sy, ylim, sy { + x =j # y coord. from lower left + xc = x + x_corner + yc = y + y_corner + # Obtain ra and dec from this grid (w/r to the original lower + # left corner) given the original plate center. + call ccgseq (plt_centre_ra, plt_centre_dec, + plt_centre_x, plt_centre_y, + x_pixel_size, y_pixel_size, plate_scale, amd_x, amd_y, + xc, yc, mag, color, ra, dec) + + # Calculate xi and eta given the new plate center + call treqst (im_cen_ra, im_cen_dec, ra, dec, Memd[xi], Memd[eta]) + xi = xi + 1 + eta = eta + 1 + + # pixel to mm from the new plate center, notice x, y are + # w/r to the new lower left corner +# xref = (new_plt_centre_x - x * x_pixel_size) / 1000. + xref = (x * x_pixel_size - new_plt_centre_x) / 1000. + yref = (y * y_pixel_size - new_plt_centre_y) / 1000. + + # form normal equations for the model + # xi = a*xref + b*yref + # eta = c*yref + d*xref + # + Memd[x_arr+k] = xref # XAR(j,1) + Memd[x_arr+k+nxy] = yref # XAR(j,2) + Memd[y_arr+k] = yref # YAR(i,1) + Memd[y_arr+k+nxy] = xref # YAR(i,2) + + k = k + 1 + ww[k] = 1.0 + } + } + + # calculate coefficients now + + npts = k + nterms = 2 + call salloc (u, npts*nterms, TY_DOUBLE) + call salloc (v, nterms*nterms, TY_DOUBLE) + call salloc (w, nterms, TY_DOUBLE) + call salloc (cvm, nterms*nterms, TY_DOUBLE) + call fitsvd (Memd[x_arr], Memd[xip], ww, npts, x_coeff, + nterms, Memd[u], Memd[v], Memd[w], xchisqr) + call varsvd (Memd[v], nterms, Memd[w], Memd[cvm], nterms) + do i =1, nterms + x_sigma[i] = sqrt(Memd[cvm+(i-1)+(i-1)*nterms]) + call fitsvd (Memd[y_arr], Memd[etap], ww, npts, y_coeff, + nterms, Memd[u], Memd[v], Memd[w], ychisqr) + call varsvd (Memd[v], nterms, Memd[w], Memd[cvm], nterms) + do i =1, nterms + y_sigma[i] = sqrt(Memd[cvm+(i-1)+(i-1)*nterms]) + + # degrees/pixel = (arcsec/mm)*(mm/pixel)*(degrees/arcsec) + cd_matrix[1] = x_coeff[1]*(x_pixel_size/1000.)*(1/3600.) + cd_matrix[2] = x_coeff[2]*(y_pixel_size/1000.)*(1/3600.) + cd_matrix[3] = y_coeff[2]*(y_pixel_size/1000.)*(1/3600.) + cd_matrix[4] = y_coeff[1]*(x_pixel_size/1000.)*(1/3600.) + + call sfree (sp) +end diff --git a/vo/votools/gasplib/ccgseq.x b/vo/votools/gasplib/ccgseq.x new file mode 100644 index 00000000..f0e071f3 --- /dev/null +++ b/vo/votools/gasplib/ccgseq.x @@ -0,0 +1,94 @@ +# CCGSEQ -- Routine for computing RA and Dec for a given X,Y pixel +# position on a GSSS plate. + +procedure ccgseq (plate_centre_ra, plate_centre_dec, plate_centre_x, + plate_centre_y, x_pixel_size, y_pixel_size, + plate_scale, amd_x, amd_y, object_x, object_y, + object_mag, object_col, object_ra, object_dec) + +double plate_centre_ra #Plate Right Ascension (radians) +double plate_centre_dec #Plate Declination (radians) +double plate_centre_x #Position used in soln.(microns) +double plate_centre_y # " +double x_pixel_size #Scan pixel size (microns) +double y_pixel_size # " +double plate_scale #Plate scale (arcsec/mm) +double amd_x[20] #Plate model coefficients +double amd_y[20] # " +double object_x #Pixel position for object +double object_y # " +double object_mag #Magnitude +double object_col #Colour +double object_ra #Object RA (radians) +double object_dec # " Dec + +double x #Position from centre (mm) +double y # " +double xi_object #Standard coords (arcsec) +double eta_object # " +double p1,p2,p3,p4 + +begin + # Convert x,y from pixels to mm measured from the plate centre + + x = (plate_centre_x - object_x * x_pixel_size) / 1000.0d0 + y = (object_y * y_pixel_size - plate_centre_y) / 1000.0d0 + + # Compute standard coordinates from x,y and plate model coefficients + + p1 = amd_x(1) *x + + amd_x(2) *y + + amd_x(3) + + amd_x(4) *x**2 + + amd_x(5) *x*y + + amd_x(6) *y**2 + + p2 = amd_x(7) *(x**2+y**2) + + amd_x(8) *x**3 + + amd_x(9) *x**2*y + + amd_x(10) *x*y**2 + + amd_x(11) *y**3 + + p3 = amd_x(12) *x*(x**2+y**2) + + amd_x(13) *x*(x**2+y**2)**2 + + amd_x(14) *object_mag + + amd_x(15) *object_mag**2 + + amd_x(16) *object_mag**3 + + p4 = amd_x(17) *object_mag*x + + amd_x(18) *object_mag*(x**2+y**2) + + amd_x(19) *object_mag*x*(x**2+y**2) + + amd_x(20) *object_col + + xi_object = p1 + p2 + p3 + p4 + + p1 = amd_y(1) *y + + amd_y(2) *x + + amd_y(3) + + amd_y(4) *y**2 + + amd_y(5) *x*y + + amd_y(6) *x**2 + + p2 = amd_y(7) *(x**2+y**2) + + amd_y(8) *y**3 + + amd_y(9) *y**2*x + + amd_y(10) *y*x**2 + + amd_y(11) *x**3 + + p3 = amd_y(12) *y*(x**2+y**2) + + amd_y(13) *y*(x**2+y**2)**2 + + amd_y(14) *object_mag + + amd_y(15) *object_mag**2 + + amd_y(16) *object_mag**3 + + p4 = amd_y(17) *object_mag*y + + amd_y(18) *object_mag*(x**2+y**2) + + amd_y(19) *object_mag*y*(x**2+y**2) + + amd_y(20) *object_col + + eta_object = p1 + p2 + p3 + p4 + + call trsteq (plate_centre_ra, plate_centre_dec, + xi_object, eta_object, object_ra, object_dec) + +end diff --git a/vo/votools/gasplib/ccgsxy.x b/vo/votools/gasplib/ccgsxy.x new file mode 100644 index 00000000..588dacef --- /dev/null +++ b/vo/votools/gasplib/ccgsxy.x @@ -0,0 +1,200 @@ +define MAX_ITERATIONS 20 +# CCGSXY -- Routine for computing X,Y pixel position from a given RA and Dec +# on a GSSS plate. + +procedure ccgsxy (plate_centre_ra, + plate_centre_dec, + plate_centre_x, + plate_centre_y, + x_pixel_size, + y_pixel_size, + plate_scale, + amd_x, + amd_y, + object_x, + object_y, + object_mag, + object_col, + object_ra, + object_dec) + + +double plate_centre_ra #Plate Right Ascension (radians) +double plate_centre_dec #Plate Declination (radians) +double plate_centre_x #Position used in soln.(microns) +double plate_centre_y # " +double x_pixel_size #Scan pixel size (microns) +double y_pixel_size # " +double plate_scale #Plate scale (arcsec/mm) +double amd_x[20] #Plate model coefficients +double amd_y[20] # " +double object_x #Pixel position for object +double object_y # " +double object_mag #Object magnitude +double object_col #Colour +double object_ra #Object RA (radians) +double object_dec # " Dec +double xi_object #Standard coords (arcsec) +double eta_object # " + +double x #Position from centre (mm) +double y # " +double delta_x #correction to x +double delta_y #correction to y +int n_iterations #no.iterations used + +double f #Function for x model +double g #Function for y model +double fx #Deriv. x model wrt x +double gx #Deriv. y model wrt x +double fy #Deriv. x model wrt y +double gy #deriv. y model wrt y + +int ierr #Stats flag + +begin + + # + # Convert RA and Dec to standard coordinates + # + call treqst (plate_centre_ra, plate_centre_dec, + object_ra, object_dec, xi_object, eta_object) + + # Iterate by Newtons method + n_iterations = 0 + ierr = 0 + + # Get initial estimate of x,y + x = xi_object / plate_scale + y = eta_object / plate_scale + + while (ierr == 0) { + n_iterations=n_iterations+1 + + # + # X plate model + # + f = amd_x(1)*x + + amd_x(2)*y + + amd_x(3) + + amd_x(4)*x*x + + amd_x(5)*x*y + + amd_x(6)*y*y + + amd_x(7)*(x*x+y*y) + + amd_x(8)*x*x*x + + amd_x(9)*x*x*y + + amd_x(10)*x*y*y + + amd_x(11)*y*y*y + + amd_x(12)*x*(x*x+y*y) + + amd_x(13)*x*(x*x+y*y)**2 + + amd_x(14)*object_mag + + amd_x(15)*object_mag**2 + + amd_x(16)*object_mag**3 + + amd_x(17)*object_mag*x + + amd_x(18)*object_mag*(x*x+y*y) + + amd_x(19)*object_mag*x*(x*x+y*y) + + amd_x(20)*object_col - + xi_object + # + # Derivative of X model wrt x + # + fx = amd_x(1) + + amd_x(4)*2.0*x + + amd_x(5)*y + + amd_x(7)*2.0*x + + amd_x(8)*3.0*x*x + + amd_x(9)*2.0*x*y + + amd_x(10)*y*y + + amd_x(12)*(3.0*x*x+y*y) + + amd_x(13)*(5.0*x**4+6.0*x*x*y*y+y**4) + + amd_x(17)*object_mag + + amd_x(18)*object_mag*2.0*x + + amd_x(19)*object_mag*(3.0*x*x+y*y) + # + # Derivative of X model wrt y + # + fy = amd_x(2) + + amd_x(5)*x + + amd_x(6)*2.0*y + + amd_x(7)*2.0*y + + amd_x(9)*x*x + + amd_x(10)*x*2.0*y + + amd_x(11)*3.0*y*y + + amd_x(12)*2.0*x*y + + amd_x(13)*4.0*x*y*(x*x+y*y) + + amd_x(18)*object_mag*2.0*y + + amd_x(19)*object_mag*2.0*x*y + # + # Y plate model + # + g = amd_y(1)*y + + amd_y(2)*x + + amd_y(3) + + amd_y(4)*y*y + + amd_y(5)*y*x + + amd_y(6)*x*x + + amd_y(7)*(x*x+y*y) + + amd_y(8)*y*y*y + + amd_y(9)*y*y*x + + amd_y(10)*y*x*x + + amd_y(11)*x*x*x + + amd_y(12)*y*(x*x+y*y) + + amd_y(13)*y*(x*x+y*y)**2 + + amd_y(14)*object_mag + + amd_y(15)*object_mag**2 + + amd_y(16)*object_mag**3 + + amd_y(17)*object_mag*y + + amd_y(18)*object_mag*(x*x+y*y) + + amd_y(19)*object_mag*y*(x*x+y*y) + + amd_y(20)*object_col - + eta_object + # + # Derivative of Y model wrt x + # + gx = amd_y(2) + + amd_y(5)*y + + amd_y(6)*2.0*x + + amd_y(7)*2.0*x + + amd_y(9)*y*y + + amd_y(10)*y*2.0*x + + amd_y(11)*3.0*x*x + + amd_y(12)*2.0*x*y + + amd_y(13)*4.0*x*y*(x*x+y*y) + + amd_y(18)*object_mag*2.0*x + + amd_y(19)*object_mag*y*2.0*x + # + # Derivative of Y model wrt y + # + gy = amd_y(1) + + amd_y(4)*2.0*y + + amd_y(5)*x + + amd_y(7)*2.0*y + + amd_y(8)*3.0*y*y + + amd_y(9)*2.0*y*x + + amd_y(10)*x*x + + amd_y(12)*3.0*y*y + + amd_y(13)*(5.0*y**4+6.0*x*x*y*y+x**4) + + amd_y(17)*object_mag + + amd_y(18)*object_mag*2.0*y + + amd_y(19)*object_mag*(x*x+3.0*y*y) + + delta_x = (-f * gy + g * fy) / (fx * gy - fy * gx) + delta_y = (-g * fx + f * gx) / (fx * gy - fy * gx) + x = x + delta_x + y = y + delta_y + + if (dmax1 (dabs(delta_x), dabs(delta_y), + dabs(f), dabs(g)) < 1.e-5) + ierr=1 + + if (n_iterations == MAX_ITERATIONS) + ierr=2 + + } + + # Convert x,y from mm about plate centre to pixels + + object_x = (plate_centre_x - x*1000.0d0) / x_pixel_size + object_y = (plate_centre_y + y*1000.0d0) / y_pixel_size + +end diff --git a/vo/votools/gasplib/dcmpsv.f b/vo/votools/gasplib/dcmpsv.f new file mode 100644 index 00000000..1e0c9bac --- /dev/null +++ b/vo/votools/gasplib/dcmpsv.f @@ -0,0 +1,230 @@ + subroutine dcmpsv (a,m,n,w,v) + parameter (nmax=1000) + real*8 a(m,n),w(n),v(n,n),rv1(nmax) + real*8 c, g, f, h, s, y, z, x, scale, anorm + + g=0.0 + scale=0.0 + anorm=0.0 + do i=1,n + l=i+1 + rv1(i)=scale*g + g=0.0 + s=0.0 + scale=0.0 + if (i.le.m) then + do k=i,m + scale=scale+dabs(a(k,i)) + enddo + if (scale.ne.0.0) then + do k=i,m + a(k,i)=a(k,i)/scale + s=s+a(k,i)*a(k,i) + enddo + f=a(i,i) + g=-dsign(dsqrt(s),f) + h=f*g-s + a(i,i)=f-g + if (i.ne.n) then + do j=l,n + s=0.0 + do k=i,m + s=s+a(k,i)*a(k,j) + enddo + f=s/h + do k=i,m + a(k,j)=a(k,j)+f*a(k,i) + enddo + enddo + endif + do k= i,m + a(k,i)=scale*a(k,i) + enddo + endif + endif + w(i)=scale *g + g=0.0 + s=0.0 + scale=0.0 + if ((i.le.m).and.(i.ne.n)) then + do k=l,n + scale=scale+dabs(a(i,k)) + enddo + if (scale.ne.0.0) then + do k=l,n + a(i,k)=a(i,k)/scale + s=s+a(i,k)*a(i,k) + enddo + f=a(i,l) + g=-dsign(dsqrt(s),f) + h=f*g-s + a(i,l)=f-g + do k=l,n + rv1(k)=a(i,k)/h + enddo + if (i.ne.m) then + do j=l,m + s=0.0 + do k=l,n + s=s+a(j,k)*a(i,k) + enddo + do k=l,n + a(j,k)=a(j,k)+s*rv1(k) + enddo + enddo + endif + do k=l,n + a(i,k)=scale*a(i,k) + enddo + endif + endif + anorm=dmax1(anorm,(dabs(w(i))+dabs(rv1(i)))) + enddo + do i=n,1,-1 + if (i.lt.n) then + if (g.ne.0.0) then + do j=l,n + v(j,i)=(a(i,j)/a(i,l))/g + enddo + do j=l,n + s=0.0 + do k=l,n + s=s+a(i,k)*v(k,j) + enddo + do k=l,n + v(k,j)=v(k,j)+s*v(k,i) + enddo + enddo + endif + do j=l,n + v(i,j)=0.0 + v(j,i)=0.0 + enddo + endif + v(i,i)=1.0 + g=rv1(i) + l=i + enddo + do i=n,1,-1 + l=i+1 + g=w(i) + if (i.lt.n) then + do j=l,n + a(i,j)=0.0 + enddo + endif + if (g.ne.0.0) then + g=1.0/g + if (i.ne.n) then + do j=l,n + s=0.0 + do k=l,m + s=s+a(k,i)*a(k,j) + enddo + f=(s/a(i,i))*g + do k=i,m + a(k,j)=a(k,j)+f*a(k,i) + enddo + enddo + endif + do j=i,m + a(j,i)=a(j,i)*g + enddo + else + do j= i,m + a(j,i)=0.0 + enddo + endif + a(i,i)=a(i,i)+1.0 + enddo + do k=n,1,-1 + do its=1,30 + do l=k,1,-1 + nm=l-1 + if ((dabs(rv1(l))+anorm).eq.anorm) go to 2 + if ((dabs(w(nm))+anorm).eq.anorm) go to 1 + enddo +1 c=0.0 + s=1.0 + do i=l,k + f=s*rv1(i) + if ((dabs(f)+anorm).ne.anorm) then + g=w(i) + h=dsqrt(f*f+g*g) + w(i)=h + h=1.0/h + c= (g*h) + s=-(f*h) + do j=1,m + y=a(j,nm) + z=a(j,i) + a(j,nm)=(y*c)+(z*s) + a(j,i)=-(y*s)+(z*c) + enddo + endif + enddo +2 z=w(k) + if (l.eq.k) then + if (z.lt.0.0) then + w(k)=-z + do j=1,n + v(j,k)=-v(j,k) + enddo + endif + go to 3 + endif + if (its.eq.30) pause 'nO CONVERGENCE IN 30 ITERATIONS' + x=w(l) + nm=k-1 + y=w(nm) + g=rv1(nm) + h=rv1(k) + f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y) + g=dsqrt(f*f+1.0) + f=((x-z)*(x+z)+h*((y/(f+dsign(g,f)))-h))/x + c=1.0 + s=1.0 + do j=l,nm + i=j+1 + g=rv1(i) + y=w(i) + h=s*g + g=c*g + z=dsqrt(f*f+h*h) + rv1(j)=z + c=f/z + s=h/z + f= (x*c)+(g*s) + g=-(x*s)+(g*c) + h=y*s + y=y*c + do nm=1,n + x=v(nm,j) + z=v(nm,i) + v(nm,j)= (x*c)+(z*s) + v(nm,i)=-(x*s)+(z*c) + enddo + z=sqrt(f*f+h*h) + w(j)=z + if (z.ne.0.0) then + z=1.0/z + c=f*z + s=h*z + endif + f= (c*g)+(s*y) + x=-(s*g)+(c*y) + do nm=1,m + y=a(nm,j) + z=a(nm,i) + a(nm,j)= (y*c)+(z*s) + a(nm,i)=-(y*s)+(z*c) + enddo + enddo + rv1(l)=0.0 + rv1(k)=f + w(k)=x + enddo +3 continue + enddo + return + end diff --git a/vo/votools/gasplib/eqtopix.x b/vo/votools/gasplib/eqtopix.x new file mode 100644 index 00000000..569f8699 --- /dev/null +++ b/vo/votools/gasplib/eqtopix.x @@ -0,0 +1,43 @@ +include <math.h> +define SZ_CDMTX 4 + +# EQTOPIX -- procedure to obtain pixel coordinates from the equatorial +# ones, giving the center position and scale of the frame. + +procedure eqtopix (plt_ra_cen, plt_dec_cen, plt_x_cen, plt_y_cen, + cdmtx, ra, dec, x, y) + +double plt_ra_cen # Plate centre (radians) +double plt_dec_cen # " " +double plt_x_cen # Plate center in x (pixels) +double plt_y_cen # Plate center in y (pixels) +double cdmtx[SZ_CDMTX] # CD Matrix (cd11, cd12, cd21, cd22) +double ra # Objects RA position (radians) +double dec # Objects DEC position (radians) +double x # Position from centre (pixels) +double y # " + +double xi # Standard coordinate (degrees) +double eta # " +double det, cd11, cd12, cd21, cd22 +double xi_arcs, eta_arcs # Standard coord. in arc seconds. + +begin + + cd11 = cdmtx[1] + cd12 = cdmtx[2] + cd21 = cdmtx[3] + cd22 = cdmtx[4] + + det = (cd11*cd22 - cd21*cd12) + call treqst (plt_ra_cen, plt_dec_cen, ra, dec, + xi_arcs, eta_arcs) + xi = xi_arcs/3600.0 # To degrees + eta = eta_arcs/3600.0 # To degrees + # Calculate x[i] and y[i] by solving the linear equations: + # xi = cd11(X-Xo) + cd12(Y-Yo) + # eta = cd21(X-Xo) + cd22(Y-Yo) + + x = (cd22*xi - cd12*eta)/det + plt_x_cen + y = (cd11*eta - cd21*xi)/det + plt_y_cen +end diff --git a/vo/votools/gasplib/fitsvd.f b/vo/votools/gasplib/fitsvd.f new file mode 100644 index 00000000..a4bb4350 --- /dev/null +++ b/vo/votools/gasplib/fitsvd.f @@ -0,0 +1,35 @@ + subroutine fitsvd (x, y, wg, npts, coef, nterms, + * u, v, w, chisq) + parameter(nmax=1000,mmax=50,tol=1.d-14) + + real wg(npts) + real*8 x(npts,nterms), y(npts), coef(nterms), v(nterms,nterms), + * u(npts,nterms), w(nterms), b(nmax) + real*8 wmax, thresh, chisq, sum + + do i=1,npts + do j=1,nterms + u(i,j)=x(i,j)*wg(i) + enddo + b(i)=y(i)*wg(i) + enddo + call dcmpsv (u,npts,nterms,w,v) + wmax=0. + do j=1,nterms + if(w(j).gt.wmax) wmax=w(j) + enddo + thresh=tol*wmax + do j=1,nterms + if(w(j).lt.thresh) w(j)=0. + enddo + call ksbsvd (u, w, v, npts, nterms, b, coef) + chisq=0. + do i=1,npts + sum=0. + do j=1,nterms + sum=sum+coef(j)*x(i,j) + enddo + chisq=chisq+((y(i)-sum)*wg(i))**2 + enddo + return + end diff --git a/vo/votools/gasplib/gsctab.x b/vo/votools/gasplib/gsctab.x new file mode 100644 index 00000000..0f48f8dd --- /dev/null +++ b/vo/votools/gasplib/gsctab.x @@ -0,0 +1,66 @@ +include <tbset.h> + +define NUM_COLS 5 # Number of table columns + +procedure gsctab (tp, cd, nrows) + +pointer tp # GSC table descriptor +pointer cd[NUM_COLS] # Table column descriptors +int nrows # Number of rows in table + +pointer sp, tname +char colname[SZ_COLNAME,NUM_COLS] +pointer null + +pointer tbtopn() +int tbpsta() + +begin + call smark (sp) + call salloc (tname, SZ_FNAME, TY_CHAR) + call clgstr ("gsctable", Memc[tname], SZ_FNAME) + + # Open the GSC table + tp = tbtopn (Memc[tname], READ_WRITE, 0) + + # Number of rows (coordinates) in the table + nrows = tbpsta (tp, TBL_NROWS) + + call strcpy ("RA_DEG", colname[1,1], SZ_COLNAME) + call strcpy ("DEC_DEG", colname[1,2], SZ_COLNAME) + call strcpy ("x_pix", colname[1,3], SZ_COLNAME) + call strcpy ("y_pix", colname[1,4], SZ_COLNAME) + call strcpy ("valid", colname[1,5], SZ_COLNAME) + + # Find the table columns + call tbcfnd (tp, colname, cd, NUM_COLS) + + if (cd[1] <= 0) + call error (0, "No R.A. column found in table") + + if (cd[2] <= 0) + call error (0, "No Dec. column found in table") + + if (cd[3] <= 0) + # Create the output (pixel coordinate) column + call tbcdef (tp, cd[3], colname[1,3], "pixels", + "%6.1f", TY_REAL, 0, 1) + + if (cd[4] <= 0) + # Create the output (pixel coordinate) column + call tbcdef (tp, cd[4], colname[1,4], "pixels", + "%6.1f", TY_REAL, 0, 1) + + if (cd[5] <= 0) { + # Create the output (valid flag) column + call tbcdef (tp, cd[5], colname[1,5], EOS, + "%1b", TY_INT, 0, 1) + + # Initialize to false + call salloc (null, nrows, TY_INT) + call amovki (NO, Memi[null], nrows) + call tbcpti (tp, cd[5], Memi[null], 1, nrows) + } + + call sfree (sp) +end diff --git a/vo/votools/gasplib/ksbsvd.f b/vo/votools/gasplib/ksbsvd.f new file mode 100644 index 00000000..d93609e9 --- /dev/null +++ b/vo/votools/gasplib/ksbsvd.f @@ -0,0 +1,24 @@ + subroutine ksbsvd (u,w,v,m,n,b,x) + parameter (nmax=1000) + real*8 u(m,n),w(n),v(n,n),b(m),x(n),tmp(nmax) + real*8 s + + do j=1,n + s=0. + if(w(j).ne.0.)then + do i=1,m + s=s+u(i,j)*b(i) + enddo + s=s/w(j) + endif + tmp(j)=s + enddo + do j=1,n + s=0. + do jj=1,n + s=s+v(j,jj)*tmp(jj) + enddo + x(j)=s + enddo + return + end diff --git a/vo/votools/gasplib/mkpkg b/vo/votools/gasplib/mkpkg new file mode 100644 index 00000000..0da7b00c --- /dev/null +++ b/vo/votools/gasplib/mkpkg @@ -0,0 +1,22 @@ +$checkout libpkg.a ../ +$update libpkg.a +$checkin libpkg.a ../ +$exit + +libpkg.a: + ccgseq.x + ccgsxy.x + eqtopix.x <math.h> + gsctab.x <tbset.h> + pixtoeq.x <math.h> + rdaslf.x + rdxy.x + trsteq.x <math.h> + treqst.x + calcds.x <math.h> + regren.f + fitsvd.f + dcmpsv.f + ksbsvd.f + varsvd.f + ; diff --git a/vo/votools/gasplib/pixtoeq.x b/vo/votools/gasplib/pixtoeq.x new file mode 100644 index 00000000..123403f0 --- /dev/null +++ b/vo/votools/gasplib/pixtoeq.x @@ -0,0 +1,37 @@ +include <math.h> +define SZ_CDMTX 4 + +# PHYSTOPIX -- procedure to obtain pixel coordinates from the equatorial +# ones, giving the center position and scale of the frame. + +procedure pixtoeq (plt_ra_cen, plt_dec_cen, plt_x_cen, plt_y_cen, + cdmtx, x, y, ra, dec) + +double plt_ra_cen # Plate centre (radians) +double plt_dec_cen # " " +double plt_x_cen # Plate center in x (pixels) +double plt_y_cen # Plate center in y (pixels) +double cdmtx[SZ_CDMTX] # CD Matrix (cd11, cd12, cd21, cd22) +double ra # Objects RA position (radians) +double dec # Objects DEC position (radians) +double x # Position from centre (pixels) +double y # " + +double xi # Standard coordinate (degrees) +double eta # " +double cd11, cd12, cd21, cd22 +double xi_arcs, eta_arcs # Standard coord. in arc seconds. + +begin + + cd11 = cdmtx[1] + cd12 = cdmtx[2] + cd21 = cdmtx[3] + cd22 = cdmtx[4] + + xi = cd11*(x-plt_x_cen) + cd12*(y-plt_y_cen) + eta = cd21*(x-plt_x_cen) + cd22*(y-plt_y_cen) + xi_arcs = xi*3600.0 # to arcs + eta_arcs = eta*3600.0 + call trsteq (plt_ra_cen, plt_dec_cen, xi_arcs, eta_arcs, ra, dec) +end diff --git a/vo/votools/gasplib/rdaslf.x b/vo/votools/gasplib/rdaslf.x new file mode 100644 index 00000000..3cd57c0b --- /dev/null +++ b/vo/votools/gasplib/rdaslf.x @@ -0,0 +1,131 @@ +# code: Nelson Zarate January 89 + +define SZ_KEYWORD 8 +define NTERMS_MODEL 20 + +# RD_ASLF -- Read astrometry solution parameters obtained with program +# 'pltsol' + +procedure rd_aslf (im, fd, namd, icd, crpix1, crpix2, crval1, + crval2, x_pixel_size, y_pixel_size, plate_scale, + namdx, namdy) + +pointer im # image file descriptor +pointer fd # Ascii file descriptor +bool namd # Are new coefficient from pltsol in header? +bool icd # Use CD values? +double crpix1, crpix2 # Reference pixel +double crval1, crval2 # Reference equatorial position (degrees) +double x_pixel_size # Y pixel size (microns) +double y_pixel_size # X pixel size (microns) +double plate_scale # [arc_sec/mm] +double namdx[NTERMS_MODEL] # New coeff in x direction +double namdy[NTERMS_MODEL] # New coeff in y direction + +char buf[SZ_KEYWORD] +char outstr[SZ_LINE] +int nch, fscan(), getline(), i, j, k, ic, ctoi(), ip, ii +int strncmp(), ctowrd(), ctod(), idb_find() +double imgetd() +pointer pp, rp + +begin + + # Clear the coefficient arrays. + call aclrd (namdx, NTERMS_MODEL) + call aclrd (namdy, NTERMS_MODEL) + + if (im == -1) { + nch = getline (fd, outstr) # Title line + nch = getline (fd, outstr) # Blank line + nch = fscan (fd) + call gargwrd (buf, SZ_KEYWORD) + call gargd (crpix1) + call gargstr (outstr, SZ_LINE) + nch = fscan (fd) + call gargwrd (buf, SZ_KEYWORD) + call gargd (crpix2) + call gargstr (outstr, SZ_LINE) + nch = fscan (fd) + call gargwrd (buf, SZ_KEYWORD) + call gargd (crval1) + call gargstr (outstr, SZ_LINE) + nch = fscan (fd) + call gargwrd (buf, SZ_KEYWORD) + call gargd (crval2) + call gargstr (outstr, SZ_LINE) + nch = fscan (fd) + call gargwrd (buf, SZ_KEYWORD) + call gargd (x_pixel_size) + call gargstr (outstr, SZ_LINE) + nch = fscan (fd) + call gargwrd (buf, SZ_KEYWORD) + call gargd (y_pixel_size) + call gargstr (outstr, SZ_LINE) + nch = fscan (fd) + call gargwrd (buf, SZ_KEYWORD) + call gargd (plate_scale) + call gargstr (outstr, SZ_LINE) + + do i = 1, NTERMS_MODEL { # read XCOEFF + k = 7 # starting character number for digits + nch = fscan (fd) + call gargwrd (buf, SZ_KEYWORD) + nch = ctoi (buf, k, ic) + if (strncmp (buf, "YCOEFF", 6) == 0) { + call gargd (namdy[ic]) + break + } + call gargd (namdx[ic]) + call gargstr (outstr, SZ_LINE) + } + do j = 2, NTERMS_MODEL { + k = 7 # starting character number for digits + nch = fscan (fd) + if (nch == EOF) + break + call gargwrd (buf, SZ_KEYWORD) + nch = ctoi (buf, k, ic) + # in case there is more than one set of coeeficients, take only + # the first one. + if (strncmp (buf, "XCOEFF", 6) == 0) + break + call gargd (namdy[ic]) + call gargstr (outstr, SZ_LINE) + } + } else { + crpix1 = imgetd (im, "CRPIX1") + crpix2 = imgetd (im, "CRPIX2") + crval1 = imgetd (im, "CRVAL1") + crval2 = imgetd (im, "CRVAL2") + if (namd) { + x_pixel_size = imgetd (im, "XPIXELSZ") + y_pixel_size = imgetd (im, "YPIXELSZ") + plate_scale = imgetd (im, "PLTSCALE") + } + + # Are new coefficients from 'pltsol' task in header + if (namd) { + do k = 1, NTERMS_MODEL { + namdx[k] = 0.0 + namdy[k] = 0.0 + } + call strcpy ("NAMDX?", buf, SZ_KEYWORD) + if (idb_find (im, buf, rp) == 0) + call error (13,"New solution coefficients (NAMDX*) not found") + for (pp=rp; Memc[pp] != EOS; pp=pp+81) { + k = 6 + ii = 11 # starting char number of keyword value + ip = 1 + # get keyword + nch = ctowrd (Memc[pp], ip, buf, SZ_KEYWORD) + # get keyword sequence + nch = ctoi (buf, k, ic) + if (strncmp (buf, "NAMDX", 5) == 0) + nch = ctod (Memc[pp], ii, namdx[ic]) + else if (strncmp (buf, "NAMDY", 5) == 0) + nch = ctod (Memc[pp], ii, namdy[ic]) + } + } + } +end diff --git a/vo/votools/gasplib/rdxy.x b/vo/votools/gasplib/rdxy.x new file mode 100644 index 00000000..360f9cd0 --- /dev/null +++ b/vo/votools/gasplib/rdxy.x @@ -0,0 +1,43 @@ +include <mach.h> + +# RDXY -- procedure to get (x,y) pixel values from a line. +# The line can have any number of columns but only up to 162 +# characters per line. +# +procedure rdxy (line, icol, x, y) + +char line[SZ_LINE] # input line with (x,y) in it +int icol[2] # X, Y column numbers within the line +double x # X pixel value +double y # Y pixel value + +pointer sp, pp +int i, ip, ic, ctowrd(), ctod() +int maxcol, nchar + +begin + + # find the right most column to read from buffer + maxcol = max (icol[1], icol[2]) + + call smark (sp) + call salloc (pp, maxcol*MAX_DIGITS, TY_CHAR) + + ip = 1 + ic = pp + do i = 1, maxcol { + nchar = ctowrd (line, ip, Memc[ic], MAX_DIGITS) + # store word in equal length array + call amovkc (" ", Memc[ic+nchar], MAX_DIGITS-nchar) + Memc[ic+MAX_DIGITS] = EOS + ic = ic + MAX_DIGITS+1 + } + + # get the output parameters + + nchar = ctod (Memc[pp], (icol(1)-1)*(MAX_DIGITS+1)+1, x) + nchar = ctod (Memc[pp], (icol(2)-1)*(MAX_DIGITS+1)+1, y) + + call sfree (sp) + +end diff --git a/vo/votools/gasplib/regren.f b/vo/votools/gasplib/regren.f new file mode 100644 index 00000000..7433b15f --- /dev/null +++ b/vo/votools/gasplib/regren.f @@ -0,0 +1,309 @@ +c The following is the IDL version of the regression routine by Bevington +c Some modification have been made to the book version in regard to the input +c X array and the data type of some of the arrays. + +c FUNCTION REGRESS,X,Y,W,YFIT,A0,SIGMA,FTEST,R,RMUL,CHISQ +c; +c;+ +c; NAME: +c; REGRESS +c; PURPOSE: +c; Multiple linear regression fit. +c; Fit the function: +c; Y(i) = A0 + A(0)*X(0,i) + A(1)*X(1,i) + ... + +c; A(Nterms-1)*X(Nterms-1,i) +c; CATEGORY: +c; G2 - Correlation and regression analysis. +c; CALLING SEQUENCE: +c; Coeff = REGRESS(X,Y,W,YFIT,A0,SIGMA,FTEST,RMUL,CHISQ) +c; INPUTS: +c; X = array of independent variable data. X must +c; be dimensioned (Nterms, Npoints) where there are Nterms +c; coefficients to be found (independent variables) and +c; Npoints of samples. +c; Y = vector of dependent variable points, must +c; have Npoints elements. +c; W = vector of weights for each equation, must +c; be a Npoints elements vector. For no +c; weighting, set w(i) = 1., for instrumental weighting +c; w(i) = 1/standard_deviation(Y(i)), for statistical +c; weighting w(i) = 1./Y(i) +c; +c; OUTPUTS: +c; Function result = coefficients = vector of +c; Nterms elements. Returned as a column +c; vector. +c; +c; OPTIONAL OUTPUT PARAMETERS: +c; Yfit = array of calculated values of Y, Npoints +c; elements. +c; A0 = Constant term. +c; Sigma = Vector of standard deviations for +c; coefficients. +c; Ftest = value of F for test of fit. +c; Rmul = multiple linear correlation coefficient. +c; R = Vector of linear correlation coefficient. +c; Chisq = Reduced weighted chi squared. +c; COMMON BLOCKS: +c; None. +c; SIDE EFFECTS: +c; None. +c; RESTRICTIONS: +c; None. +c; PROCEDURE: +c; Adapted from the program REGRES, Page 172, +c; Bevington, Data Reduction and Error Analysis for the +c; Physical Sciences, 1969. +c; +c; MODIFICATION HISTORY: +c; Written, DMS, RSI, September, 1982. +c;- +c; +c SY = SIZE(Y) ;GET DIMENSIONS OF X AND Y. +c SX = SIZE(X) +c IF (N_ELEMENTS(W) NE SY(1)) OR (SX(0) NE 2) OR (SY(1) NE SX(2)) THEN BEGIN +c PRINT,'REGRESS - Incompatible arrays' +c RETURN,0 +c ENDIF +c; +c NTERM = SX(1) ;# OF TERMS +c NPTS = SY(1) ;# OF OBSERVATIONS +c ; +c SW = TOTAL(W) ;SUM OF WEIGHTS +c YMEAN = TOTAL(Y*W)/SW ;Y MEAN +c XMEAN = (X * (REPLICATE(1.,NTERM) # W)) # REPLICATE(1./SW,NPTS) +c WMEAN = SW/NPTS +c WW = W/WMEAN +c ; +c NFREE = NPTS-1 ;DEGS OF FREEDOM +c SIGMAY = SQRT(TOTAL(WW * (Y-YMEAN)^2)/NFREE) ;W*(Y(I)-YMEAN) +c XX = X- XMEAN # REPLICATE(1.,NPTS) ;X(J,I) - XMEAN(I) +c WX = REPLICATE(1.,NTERM) # WW * XX ;W(I)*(X(J,I)-XMEAN(I)) +c SIGMAX = SQRT( XX*WX # REPLICATE(1./NFREE,NPTS)) ;W(I)*(X(J,I)-XM)*(X(K,I)-XM) +c R = WX #(Y - YMEAN) / (SIGMAX * SIGMAY * NFREE) +c ARRAY = INVERT((WX # TRANSPOSE(XX))/(NFREE * SIGMAX #SIGMAX)) +c A = (R # ARRAY)*(SIGMAY/SIGMAX) ;GET COEFFICIENTS +c YFIT = A # X ;COMPUTE FIT +c A0 = YMEAN - TOTAL(A*XMEAN) ;CONSTANT TERM +c YFIT = YFIT + A0 ;ADD IT IN +c FREEN = NPTS-NTERM-1 > 1 ;DEGS OF FREEDOM, AT LEAST 1. +c CHISQ = TOTAL(WW*(Y-YFIT)^2)*WMEAN/FREEN ;WEIGHTED CHI SQUARED +c SIGMA = SQRT(ARRAY(INDGEN(NTERM)*(NTERM+1))/WMEAN/(NFREE*SIGMAX^2)) ;ERROR TERM +c RMUL = TOTAL(A*R*SIGMAX/SIGMAY) ;MULTIPLE LIN REG COEFF +c IF RMUL LT 1. THEN FTEST = RMUL/NTERM / ((1.-RMUL)/FREEN) ELSE FTEST=1.E6 +c RMUL = SQRT(RMUL) +c RETURN,A +c END + + subroutine regren (x, ndim1, ndim2, y, weight, npts, nterms, yfit, + *a0, a, sigmaa, chisqr) + double precision sum, ymean, sigma, chisq + integer npts, nterms + double precision x(ndim1,ndim2),y(1),yfit(1) + double precision r(20), array(20,20), sigmax(20), xmean(20) + double precision chisqr, a0, a(1) + real weight(1), sigmaa(1) + real sigma0, ftest, freen, free1, rmul + real fnpts, det, varnce, wmean, freej + integer i, j, k +c +c initialize sums and arrays +c +11 sum=0. + ymean=0. + sigma=0. + chisq=0. + rmul=0. + do 17 i=1,npts +17 yfit(i)=0. +21 do 28 j=1,nterms + xmean(j)=0. + sigmax(j)=0. + r(j)=0. + a(j)=0. + sigmaa(j)=0. + do 28 k=1,nterms +28 array(j,k)=0. +c +c accumulate weighted sums +c +30 do 50 i=1,npts + sum=sum+weight(i) + ymean=ymean+weight(i)*y(i) + do 44 j=1,nterms +44 xmean(j)=xmean(j)+weight(i)*x(j,i) +50 continue +51 ymean=ymean/sum + do 53 j=1,nterms +53 xmean(j)=xmean(j)/sum + fnpts=npts + wmean=sum/fnpts + do 57 i=1,npts +57 weight(i)=weight(i)/wmean +c +c accumulate matrices r and array +c +61 do 67 i=1,npts + sigma=sigma+weight(i)*(y(i)-ymean)**2 + do 67 j=1,nterms + sigmax(j)=sigmax(j)+weight(i)*(x(j,i)-xmean(j))**2 + r(j)=r(j)+weight(i)*(x(j,i)-xmean(j))*(y(i)-ymean) + do 67 k=1,j +67 array(j,k)=array(j,k)+weight(i)*(x(j,i)-xmean(j))* + *(x(k,i)-xmean(k)) +71 free1=npts-1 +72 sigma=dsqrt(sigma/free1) + do 78 j=1,nterms +74 sigmax(j)=dsqrt(sigmax(j)/free1) + r(j)=r(j)/(free1*sigmax(j)*sigma) + do 78 k=1,j + array(j,k)=array(j,k)/(free1*sigmax(j)*sigmax(k)) +78 array(k,j)=array(j,k) +c +c invert symmetric matrix +c +81 call minv20 (array,nterms,det) + if (det) 101,91,101 +91 a0=0. + sigma0=0. + rmul=0. + chisqr=0. + ftest=0. + goto 150 +c +c calculate coefficients, fit, and chi square +c +101 a0=ymean +102 do 108 j=1,nterms + do 104 k=1,nterms +104 a(j)=a(j)+r(k)*array(j,k) +105 a(j)=a(j)*sigma/sigmax(j) +106 a0=a0-a(j)*xmean(j) +107 do 108 i=1,npts +108 yfit(i)=yfit(i)+a(j)*x(j,i) +111 do 113 i=1,npts + yfit(i)=yfit(i)+a0 +113 chisq=chisq+weight(i)*(y(i)-yfit(i))**2 + freen=npts-nterms-1 +115 chisqr=chisq*wmean/freen +c +c calculate uncertainties +c +124 varnce=chisqr +131 do 133 j=1,nterms +132 sigmaa(j)=array(j,j)*varnce/(free1*sigmax(j)**2) + sigmaa(j)=sqrt(sigmaa(j)) +133 rmul=rmul+a(j)*r(j)*sigmax(j)/sigma + freej=nterms +c +noao: When rmul = 1, the following division (stmt 135) would blow up. +c It has been changed so ftest is set to -99999. in this case. + if (1. - rmul) 135, 935, 135 +135 ftest=(rmul/freej)/((1.-rmul)/freen) + goto 136 +935 ftest = -99999. +c -noao +136 rmul=sqrt(rmul) +141 sigma0=varnce/fnpts + do 145 j=1,nterms + do 145 k=1,nterms +145 sigma0=sigma0+varnce*xmean(j)*xmean(k)*array(j,k)/ + *(free1*sigmax(j)*sigmax(k)) +146 sigma0=sqrt(sigma0) +150 return + end +c subroutine matinv.f +c +c source +c bevington, pages 302-303. +c +c purpose +c invert a symmetric matrix and calculate its determinant +c +c usage +c call matinv (array, norder, det) +c +c description of parameters +c array - input matrix which is replaced by its inverse +c norder - degree of matrix (order of determinant) +c det - determinant of input matrix +c +c subroutines and function subprograms required +c none +c +c comment +c dimension statement valid for norder up to 20 +c + subroutine minv20 (array,norder,det) + double precision array,amax,save + dimension array(20,20),ik(20),jk(20) +c +10 det=1. +11 do 100 k=1,norder +c +c find largest element array(i,j) in rest of matrix +c + amax=0. +21 do 30 i=k,norder + do 30 j=k,norder +23 if (dabs(amax)-dabs(array(i,j))) 24,24,30 +24 amax=array(i,j) + ik(k)=i + jk(k)=j +30 continue +c +c interchange rows and columns to put amax in array(k,k) +c +31 if (amax) 41,32,41 +32 det=0. + goto 140 +41 i=ik(k) + if (i-k) 21,51,43 +43 do 50 j=1,norder + save=array(k,j) + array(k,j)=array(i,j) +50 array(i,j)=-save +51 j=jk(k) + if (j-k) 21,61,53 +53 do 60 i=1,norder + save=array(i,k) + array(i,k)=array(i,j) +60 array(i,j)=-save +c +c accumulate elements of inverse matrix +c +61 do 70 i=1,norder + if (i-k) 63,70,63 +63 array(i,k)=-array(i,k)/amax +70 continue +71 do 80 i=1,norder + do 80 j=1,norder + if (i-k) 74,80,74 +74 if (j-k) 75,80,75 +75 array(i,j)=array(i,j)+array(i,k)*array(k,j) +80 continue +81 do 90 j=1,norder + if (j-k) 83,90,83 +83 array(k,j)=array(k,j)/amax +90 continue + array(k,k)=1./amax +100 det=det*amax +c +c restore ordering of matrix +c +101 do 130 l=1,norder + k=norder-l+1 + j=ik(k) + if (j-k) 111,111,105 +105 do 110 i=1,norder + save=array(i,k) + array(i,k)=-array(i,j) +110 array(i,j)=save +111 i=jk(k) + if (i-k) 130,130,113 +113 do 120 j=1,norder + save=array(k,j) + array(k,j)=-array(i,j) +120 array(i,j)=save +130 continue +140 return + end diff --git a/vo/votools/gasplib/treqst.x b/vo/votools/gasplib/treqst.x new file mode 100644 index 00000000..9e6b6548 --- /dev/null +++ b/vo/votools/gasplib/treqst.x @@ -0,0 +1,33 @@ +define ARCSEC_PER_RADIAN 206264.8062470964d0 +# TREQST -- Procedure to convert RA and Dec to standard coordinates +# given the plate centre. + +procedure treqst (plate_centre_ra, plate_centre_dec, + object_ra, object_dec, xi_object, eta_object) + +double plate_centre_ra #i: (radians) +double plate_centre_dec #i: (radians) +double object_ra #i: (radians) +double object_dec #i: (radians) +double xi_object #o: standard coordinate (rad) +double eta_object #o: standard coordinate (rad) + +double div + +begin + # Find divisor + div=(dsin(object_dec)*dsin(plate_centre_dec)+ + dcos(object_dec)*dcos(plate_centre_dec)* + dcos(object_ra-plate_centre_ra)) + + # Compute standard coords and convert to arcsec + + xi_object = dcos(object_dec)*dsin(object_ra-plate_centre_ra)* + ARCSEC_PER_RADIAN/div + + eta_object = (dsin(object_dec)*dcos(plate_centre_dec)- + dcos(object_dec)*dsin(plate_centre_dec)* + dcos(object_ra-plate_centre_ra))* + ARCSEC_PER_RADIAN/div + +end diff --git a/vo/votools/gasplib/trsteq.x b/vo/votools/gasplib/trsteq.x new file mode 100644 index 00000000..68ce7997 --- /dev/null +++ b/vo/votools/gasplib/trsteq.x @@ -0,0 +1,45 @@ +define ARCSEC_PER_RADIAN 206264.8062470964d0 + +# TRSTEQ -- procedure to compute the RA and Dec from Standard coordinates, +# given the plate centre. +# +# In rectangular coordinates the vector (1, xi, eta) points toward +# the object; the origin is the observer's location, the x-axis points +# toward the reference pixel (plate centre), the y-axis is in the direction +# of increasing right ascension, and the z-axis is in the direction of +# increasing declination. The coordinate system is then rotated by the +# declination so the x-axis passes through the equator at the RA of the +# reference pixel; the components of the vector in this coordinate system +# are used to compute (RA - reference_RA) and declination. +# +# original, written by ??? +# Phil Hodge, 15-Nov-1999 Rewrite, based on tables$lib/stxtools/xtwcs.x. + +procedure trsteq (ra0, dec0, xi, eta, ra, dec) + +double ra0, dec0 # i: RA & Dec at plate centre (radians) +double xi, eta # i: standard coordinates (arcsec) +double ra, dec # o: right ascension & declination (radians) +#-- +double xi_r, eta_r # xi & eta in radians +double x, y, z # vector (not unit length) pointing toward object +double dra # RA at (xi,eta) - RA at plate centre + +begin + xi_r = xi / ARCSEC_PER_RADIAN + eta_r = eta / ARCSEC_PER_RADIAN + + # Rotate the rectangular coordinate system of the vector (1, xi, eta) + # by the declination so the x-axis will pass through the equator. + x = cos (dec0) - eta_r * sin (dec0) + y = xi_r + z = sin (dec0) + eta_r * cos (dec0) + + if (x == 0.d0 && y == 0.d0) + dra = 0.d0 + else + dra = atan2 (y, x) + ra = ra0 + dra + + dec = atan2 (z, sqrt (x*x + y*y)) +end diff --git a/vo/votools/gasplib/varsvd.f b/vo/votools/gasplib/varsvd.f new file mode 100644 index 00000000..81e800f3 --- /dev/null +++ b/vo/votools/gasplib/varsvd.f @@ -0,0 +1,21 @@ + subroutine varsvd (v,ma,w,cvm,ncvm) + parameter (mmax=20) + real*8 v(ma,ma),w(ma),cvm(ncvm,ncvm),wti(mmax) + real*8 sum + + do i=1,ma + wti(i)=0. + if(w(i).ne.0.0d0) wti(i)=1./(w(i)*w(i)) + enddo + do i=1,ma + do j=1,i + sum=0. + do k=1,ma + sum=sum+v(i,k)*v(j,k)*wti(k) + enddo + cvm(i,j)=sum + cvm(j,i)=sum + enddo + enddo + return + end |