From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- noao/astcat/src/awcs/atmwshow.x | 129 ++++++++++ noao/astcat/src/awcs/calcds.x | 128 ++++++++++ noao/astcat/src/awcs/ccqseq.x | 95 ++++++++ noao/astcat/src/awcs/dbwcs.x | 522 ++++++++++++++++++++++++++++++++++++++++ noao/astcat/src/awcs/dcmpsv.f | 233 ++++++++++++++++++ noao/astcat/src/awcs/dsswcs.x | 300 +++++++++++++++++++++++ noao/astcat/src/awcs/fitsvd.f | 38 +++ noao/astcat/src/awcs/ksbsvd.f | 27 +++ noao/astcat/src/awcs/mkpkg | 22 ++ noao/astcat/src/awcs/parswcs.x | 251 +++++++++++++++++++ noao/astcat/src/awcs/treqst.x | 49 ++++ noao/astcat/src/awcs/trsteq.x | 64 +++++ noao/astcat/src/awcs/varsvd.f | 24 ++ 13 files changed, 1882 insertions(+) create mode 100644 noao/astcat/src/awcs/atmwshow.x create mode 100644 noao/astcat/src/awcs/calcds.x create mode 100644 noao/astcat/src/awcs/ccqseq.x create mode 100644 noao/astcat/src/awcs/dbwcs.x create mode 100644 noao/astcat/src/awcs/dcmpsv.f create mode 100644 noao/astcat/src/awcs/dsswcs.x create mode 100644 noao/astcat/src/awcs/fitsvd.f create mode 100644 noao/astcat/src/awcs/ksbsvd.f create mode 100644 noao/astcat/src/awcs/mkpkg create mode 100644 noao/astcat/src/awcs/parswcs.x create mode 100644 noao/astcat/src/awcs/treqst.x create mode 100644 noao/astcat/src/awcs/trsteq.x create mode 100644 noao/astcat/src/awcs/varsvd.f (limited to 'noao/astcat/src/awcs') diff --git a/noao/astcat/src/awcs/atmwshow.x b/noao/astcat/src/awcs/atmwshow.x new file mode 100644 index 00000000..cc39bbf3 --- /dev/null +++ b/noao/astcat/src/awcs/atmwshow.x @@ -0,0 +1,129 @@ +# AT_MWSHOW -- Print a quick summary of the current wcs. + +procedure at_mwshow (mwim, ltv, ltm, w, r, cd, ndim) + +pointer mwim # pointer to the current wcs +double ltv[ARB] # the lterm offsets +double ltm[ndim,ARB] # the lterm rotation matrix +double w[ARB] # the fits crval parameters +double r[ARB] # the fits crpix parameters +double cd[ndim,ARB] # the fits rotation matrix +int ndim # the dimension of the wcs + +int i,j +pointer sp, str +errchk mw_gwattrs() + +begin + # Allocate working space. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Print the axis banner. + call printf (" AXIS ") + do i = 1, ndim { + call printf ("%8d ") + call pargi (i) + } + call printf ("\n") + + # Print the crval parameters. + call printf (" CRVAL ") + do i = 1, ndim { + call printf ("%8g ") + call pargd (w[i]) + } + call printf ("\n") + + # Print the crpix parameters. + call printf (" CRPIX ") + do i = 1, ndim { + call printf ("%8g ") + call pargd (r[i]) + } + call printf ("\n") + + # Print the cd matrix. + do i = 1, ndim { + call printf (" CD %d ") + call pargi (i) + do j = 1, ndim { + call printf ("%8g ") + call pargd (cd[j,i]) + } + call printf ("\n") + } + + # Print the ltv parameters. + call printf (" LTV ") + do i = 1, ndim { + call printf ("%8g ") + call pargd (ltv[i]) + } + call printf ("\n") + + # Print the ltm matrix. + do i = 1, ndim { + call printf (" LTM %d ") + call pargi (i) + do j = 1, ndim { + call printf ("%8g ") + call pargd (ltm[i,j]) + } + call printf ("\n") + } + + # Print the transformation type. + call printf (" WTYPE ") + do i = 1, ndim { + iferr (call mw_gwattrs (mwim, i, "wtype", Memc[str], SZ_LINE)) + Memc[str] = EOS + call printf ("%8s ") + call pargstr (Memc[str]) + } + call printf ("\n") + + # Print the axis type. + call printf (" AXTYPE ") + do i = 1, ndim { + iferr (call mw_gwattrs (mwim, i, "axtype", Memc[str], SZ_LINE)) + Memc[str] = EOS + call printf ("%8s ") + call pargstr (Memc[str]) + } + call printf ("\n") + + # Print the units. + call printf (" UNITS ") + do i = 1, ndim { + iferr (call mw_gwattrs (mwim, i, "units", Memc[str], SZ_LINE)) + Memc[str] = EOS + call printf ("%8s ") + call pargstr (Memc[str]) + } + call printf ("\n") + + # Print the label. + call printf (" LABEL ") + do i = 1, ndim { + iferr (call mw_gwattrs (mwim, i, "label", Memc[str], SZ_LINE)) + Memc[str] = EOS + call printf ("%8s ") + call pargstr (Memc[str]) + } + call printf ("\n") + + # Print the format. + call printf (" FORMAT ") + do i = 1, ndim { + iferr (call mw_gwattrs (mwim, i, "format", Memc[str], SZ_LINE)) + Memc[str] = EOS + call printf ("%8s ") + call pargstr (Memc[str]) + } + call printf ("\n") + + call printf ("\n") + + call sfree (sp) +end diff --git a/noao/astcat/src/awcs/calcds.x b/noao/astcat/src/awcs/calcds.x new file mode 100644 index 00000000..de7ebcf2 --- /dev/null +++ b/noao/astcat/src/awcs/calcds.x @@ -0,0 +1,128 @@ +include + +define SZ_GRID 10 +define SZ_NTERMS 2 + +# CALCDS -- Procedure to calculate the values of the CD matrix from the +# GSSS plate solution and a grid of 100 tie points. This routine was +# adapted from one in stsdas$pkg/analysis/gasp/gasplib/. See the routine +# stsdas$copyright.stsdas. + +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 #I plate centre RA (radians) +double plt_centre_dec #I plate centre DEC (radians) +double plt_centre_x #I x center position (microns) +double plt_centre_y #I y center position (microns) +int x_corner #I x lower left of the extracted image +int y_corner #I y lower left of the extracted image +double x_pixel_size #I y scan pixel size (microns) +double y_pixel_size #I y scan pixel size (microns) +double plate_scale #I plate scale (arcsec / mm) +int x_size #I extracted image size x_axis (pixel) +int y_size #I extracted image size y_axis (pixel) +double im_cen_ra #I extracted image center RA (radians) +double im_cen_dec #I extracted image center DEC (radians) +double amd_x[ARB] #I XI plate solution coefficients +double amd_y[ARB] #I ETA coefficients (arsec / mm) +double cd_matrix[ARB] #O CD1_1, CD1_2, CD2_1, CD2_2 (degrees / pixel) + +double ra, dec, new_plt_centre_x, new_plt_centre_y, xref, yref, mag, color +double x_coeff[SZ_NTERMS], y_coeff[SZ_NTERMS], xchisqr, ychisqr +double x_sigma[SZ_NTERMS], y_sigma[SZ_NTERMS], x, y, xc, yc +pointer sp, xip, etap, x_arr, y_arr, ww, u, v, w, cvm +int sx, sy, xlim, ylim, nx, ny, nxy +int i, j, nterms, xi, eta, npts + +begin + # Initialize color and magnitude. + mag = 0.0d0 + color = 0.0d0 + + # Calculate new plate center in microns. + new_plt_centre_x = (x_size / 2.0d0) * x_pixel_size + new_plt_centre_y = (y_size / 2.0d0) * y_pixel_size + + call smark (sp) + call salloc (xip, SZ_GRID * SZ_GRID, TY_DOUBLE) + call salloc (etap, SZ_GRID * SZ_GRID, TY_DOUBLE) + call salloc (x_arr, SZ_NTERMS * SZ_GRID * SZ_GRID, TY_DOUBLE) + call salloc (y_arr, SZ_NTERMS * SZ_GRID * SZ_GRID, TY_DOUBLE) + call salloc (ww, SZ_GRID * SZ_GRID, TY_REAL) + + sx = max (1, x_size / SZ_GRID) + sy = max (1, y_size / SZ_GRID) + 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 + + # Compute the grid points. + npts = 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+npts] = xref # XAR(j,1) + Memd[x_arr+npts+nxy] = yref # XAR(j,2) + Memd[y_arr+npts] = yref # YAR(i,1) + Memd[y_arr+npts+nxy] = xref # YAR(i,2) + Memr[ww+npts] = 1.0 + npts = npts + 1 + } + } + + # Calculate the coefficients. + nterms = SZ_NTERMS + 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], Memr[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], Memr[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.0d0 / 3600.0d0) + cd_matrix[2] = x_coeff[2] * (y_pixel_size / 1000.0d0 / 3600.0d0) + cd_matrix[3] = y_coeff[2] * (y_pixel_size / 1000.0d0 / 3600.0d0) + cd_matrix[4] = y_coeff[1] * (x_pixel_size / 1000.0d0 / 3600.0d0) + + call sfree (sp) +end diff --git a/noao/astcat/src/awcs/ccqseq.x b/noao/astcat/src/awcs/ccqseq.x new file mode 100644 index 00000000..542ea05e --- /dev/null +++ b/noao/astcat/src/awcs/ccqseq.x @@ -0,0 +1,95 @@ +# CCGSEQ -- Routine for computing RA and Dec for a given X,Y pixel +# position on a GSSS image. Adapted from stsdas$pkg/analysis/gasp/gasplib/. +# See stsdas$copyright.stdas for copyright restrictions. + +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 #I plate Right Ascension (radians) +double plate_centre_dec #I plate Declination (radians) +double plate_centre_x #I x position used in solution (microns) +double plate_centre_y #I y position used in solution (microns) +double x_pixel_size #I scan pixel size in x (microns) +double y_pixel_size #I scan pixel size in y (microns) +double plate_scale #I plate scale (arcsec / mm) +double amd_x[ARB] #I ra plate model coefficients +double amd_y[ARB] #I dec plate model coefficinets +double object_x #I x pixel position for object +double object_y #I y pixel positions for objects +double object_mag #I object magnitude +double object_col #I object colour +double object_ra #O object ra (radians) +double object_dec #O object dec (radians) + +double x # position from center (mm) +double y # position from center (mm) +double xi_object # xi standard coordinate (arcsec) +double eta_object # eta standard coordinate (arcsec) +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/noao/astcat/src/awcs/dbwcs.x b/noao/astcat/src/awcs/dbwcs.x new file mode 100644 index 00000000..f3109050 --- /dev/null +++ b/noao/astcat/src/awcs/dbwcs.x @@ -0,0 +1,522 @@ +include +include +include +include +include + +# These should probably go into aimpars.h. + +define IMDB_WCSDICT "|wxref|wyref|wxmag|wymag|wxrot|wyrot|wraref|wdecref|\ +wproj|wsystem|" + +define IMDB_WCS_WXREF 1 +define IMDB_WCS_WYREF 2 +define IMDB_WCS_WXMAG 3 +define IMDB_WCS_WYMAG 4 +define IMDB_WCS_WXROT 5 +define IMDB_WCS_WYROT 6 +define IMDB_WCS_WLNGREF 7 +define IMDB_WCS_WLATREF 8 +define IMDB_WCS_WPROJ 9 +define IMDB_WCS_WSYSTEM 10 + + +# AT_DBWCS -- Compute a FITS WCS from an image using WCS definitions +# stored in the image surveys configuration file and transferred to the +# image query results structure. At the moment I am going to keep this +# routine simple by not worrying about the units of any quantities but the +# world coordinates of the reference point. This routine can be made more +# sophisticated later as time permits. The information is there ... + +int procedure at_dbwcs (im, res, update, verbose) + +pointer im #I the input image descriptor +pointer res #I the image query results descriptor +bool update #I update rather than list the wcs +bool verbose #I verbose mode + +double xref, yref, xmag, ymag, xrot, yrot, lngref, latref, dval +pointer sp, kfield, kname, kvalue, kunits, wtype, ctype, coo, mw +int i, ip, stat, coostat, ktype, nwcs, wkey, lngunits, latunits +double imgetd(), at_imhms() +int cq_istati(), cq_winfon(), strdic(), ctod(), ctowrd(), sk_decwcs() +bool streq() +errchk imgetd() + +begin + # Return if the input is not 2D. + + if (IM_NDIM(im) != 2) + return (ERR) + + # Allocate working space. + + call smark (sp) + call salloc (kfield, CQ_SZ_QPNAME, TY_CHAR) + call salloc (kname, CQ_SZ_QPNAME, TY_CHAR) + call salloc (kvalue, CQ_SZ_QPVALUE, TY_CHAR) + call salloc (kunits, CQ_SZ_QPUNITS, TY_CHAR) + call salloc (wtype, SZ_FNAME, TY_CHAR) + call salloc (ctype, SZ_FNAME, TY_CHAR) + + # Assume some sensible defaults, e.g. the reference point is at + # the center of the image, the orientation is the standard astronomical + # orientation with ra increasing to the left and declination increasing + # to the top, the projection is tan, the coordinate system is J2000. + + xref = (IM_LEN(im,1) + 1.0d0)/ 2.0d0 + yref = (IM_LEN(im,2) + 1.0d0)/ 2.0d0 + xmag = INDEFD + ymag = INDEFD + xrot = 180.0d0 + yrot = 0.0d0 + lngref = INDEFD + latref = INDEFD + call strcpy ("tan", Memc[wtype], SZ_FNAME) + call strcpy ("J2000", Memc[ctype], SZ_FNAME) + + # Loop over the mwcs database quantities. + + nwcs = cq_istati (res, CQNWCS) + do i = 1, nwcs { + + # Get the keyword information. + if (cq_winfon (res, i, Memc[kfield], CQ_SZ_QPNAME, Memc[kname], + CQ_SZ_QPNAME, Memc[kvalue], CQ_SZ_QPVALUE, ktype, Memc[kunits], + CQ_SZ_QPUNITS) != i) + next + + # Which keyword have we got ? + wkey = strdic (Memc[kfield], Memc[kfield], CQ_SZ_QPNAME, + IMDB_WCSDICT) + ip = 1 + switch (wkey) { + + # Get the x reference point in pixels. + case IMDB_WCS_WXREF: + if (streq (Memc[kname], "INDEF")) { + if (streq (Memc[kvalue], "INDEF")) + dval = INDEFD + else if (ctod (Memc[kvalue], ip, dval) <= 0) + dval = INDEFD + } else iferr (dval = imgetd (im, Memc[kname])) + dval = INDEFD + if (! IS_INDEFD(dval)) + xref = dval + + case IMDB_WCS_WYREF: + if (streq (Memc[kname], "INDEF")) { + if (streq (Memc[kvalue], "INDEF")) + dval = INDEFD + else if (ctod (Memc[kvalue], ip, dval) <= 0) + dval = INDEFD + } else iferr (dval = imgetd (im, Memc[kname])) + dval = INDEFD + if (! IS_INDEFD(dval)) + yref = dval + + case IMDB_WCS_WXMAG: + if (streq (Memc[kname], "INDEF")) { + if (streq (Memc[kvalue], "INDEF")) + dval = INDEFD + else if (ctod (Memc[kvalue], ip, dval) <= 0) + dval = INDEFD + } else iferr (dval = imgetd (im, Memc[kname])) + dval = INDEFD + if (! IS_INDEFD(dval)) + xmag = dval + + case IMDB_WCS_WYMAG: + if (streq (Memc[kname], "INDEF")) { + if (streq (Memc[kvalue], "INDEF")) + dval = INDEFD + else if (ctod (Memc[kvalue], ip, dval) <= 0) + dval = INDEFD + } else iferr (dval = imgetd (im, Memc[kname])) + dval = INDEFD + if (! IS_INDEFD(dval)) + ymag = dval + + case IMDB_WCS_WXROT: + if (streq (Memc[kname], "INDEF")) { + if (streq (Memc[kvalue], "INDEF")) + dval = INDEFD + else if (ctod (Memc[kvalue], ip, dval) <= 0) + dval = INDEFD + } else iferr (dval = imgetd (im, Memc[kname])) + dval = INDEFD + if (! IS_INDEFD(dval)) + xrot = dval + + case IMDB_WCS_WYROT: + if (streq (Memc[kname], "INDEF")) { + if (streq (Memc[kvalue], "INDEF")) + dval = INDEFD + else if (ctod (Memc[kvalue], ip, dval) <= 0) + dval = INDEFD + } else iferr (dval = imgetd (im, Memc[kname])) + dval = INDEFD + if (! IS_INDEFD(dval)) + yrot = dval + + case IMDB_WCS_WLNGREF: + if (streq (Memc[kname], "INDEF")) { + if (streq (Memc[kvalue], "INDEF")) + dval = INDEFD + else if (ctod (Memc[kvalue], ip, dval) <= 0) + dval = INDEFD + } else { + dval = at_imhms (im, Memc[kname]) + if (IS_INDEFD(dval)) { + iferr (dval = imgetd (im, Memc[kname])) + dval = INDEFD + } + } + if (! IS_INDEFD(dval)) + lngref = dval + lngunits = strdic (Memc[kunits], Memc[kunits], CQ_SZ_QPUNITS, + SKY_LNG_UNITLIST) + + case IMDB_WCS_WLATREF: + if (streq (Memc[kname], "INDEF")) { + if (streq (Memc[kvalue], "INDEF")) + dval = INDEFD + else if (ctod (Memc[kvalue], ip, dval) <= 0) + dval = INDEFD + } else { + dval = at_imhms (im, Memc[kname]) + if (IS_INDEFD(dval)) { + iferr (dval = imgetd (im, Memc[kname])) + dval = INDEFD + } + } + if (! IS_INDEFD(dval)) + latref = dval + latunits = strdic (Memc[kunits], Memc[kunits], CQ_SZ_QPUNITS, + SKY_LAT_UNITLIST) + + case IMDB_WCS_WPROJ: + if (streq (Memc[kname], "INDEF")) { + if (streq (Memc[kvalue], "INDEF")) + call strcpy ("tan", Memc[wtype], SZ_FNAME) + else if (ctowrd (Memc[kvalue], ip, Memc[wtype], + SZ_FNAME) <= 0) + call strcpy ("tan", Memc[wtype], SZ_FNAME) + } else iferr (call imgstr (im, Memc[kname], Memc[wtype], + SZ_FNAME)) + call strcpy ("tan", Memc[wtype], SZ_FNAME) + + case IMDB_WCS_WSYSTEM: + if (streq (Memc[kname], "INDEF")) { + if (streq (Memc[kvalue], "INDEF")) + call strcpy ("J2000", Memc[ctype], SZ_FNAME) + else if (ctowrd (Memc[kvalue], ip, Memc[ctype], + SZ_FNAME) <= 0) + call strcpy ("J2000", Memc[ctype], SZ_FNAME) + } else iferr (call imgstr (im, Memc[kname], Memc[ctype], + SZ_FNAME)) + call strcpy ("J2000", Memc[ctype], SZ_FNAME) + + default: + ; + } + } + + # Check to see of the critical quantities image scale and reference + # point are defined. Quit if they are not, otherwise update the + # header. + + if (IS_INDEFD(xmag) || IS_INDEFD(ymag) || IS_INDEFD(lngref) || + IS_INDEFD(latref)) { + + stat = ERR + + } else { + + # Open the coordinate system structure. + coostat = sk_decwcs (Memc[ctype], mw, coo, NULL) + + # Update hte header. + if (coostat == ERR || mw != NULL) { + if (mw != NULL) + call mw_close (mw) + stat = ERR + } else { + if (lngunits > 0) + call sk_seti (coo, S_NLNGUNITS, lngunits) + if (latunits > 0) + call sk_seti (coo, S_NLATUNITS, latunits) + if (verbose) + call printf (" Writing FITS wcs using image survey db\n") + call at_uwcs (im, coo, Memc[wtype], lngref, latref, xref, + yref, xmag, ymag, xrot, yrot, false, update) + stat = OK + } + + # Close the coordinate structure + if (coo != NULL) + call sk_close (coo) + } + + call sfree (sp) + return (stat) +end + + +define NEWCD Memd[ncd+(($2)-1)*ndim+($1)-1] + +# AT_UWCS -- Compute the image wcs from the user parameters. + +procedure at_uwcs (im, coo, projection, lngref, latref, xref, yref, + xscale, yscale, xrot, yrot, transpose, update) + +pointer im #I pointer to the input image +pointer coo #I pointer to the coordinate structure +char projection[ARB] #I the sky projection geometry +double lngref, latref #I the world coordinates of the reference point +double xref, yref #I the reference point in pixels +double xscale, yscale #I the x and y scale in arcsec / pixel +double xrot, yrot #I the x and y axis rotation angles in degrees +bool transpose #I transpose the wcs +bool update #I update rather than list the wcs + + +double tlngref, tlatref +int l, i, ndim, naxes, axmap, wtype, ax1, ax2, szatstr +pointer mw, sp, r, w, cd, ltm, ltv, iltm, nr, ncd, axes, axno, axval +pointer projstr, projpars, wpars, mwnew, atstr +int mw_stati(), sk_stati(), strdic(), strlen(), itoc() +pointer mw_openim(), mw_open() +errchk mw_newsystem(), mw_gwattrs() + +begin + mw = mw_openim (im) + ndim = mw_stati (mw, MW_NPHYSDIM) + + # Allocate working memory for the vectors and matrices. + call smark (sp) + call salloc (projstr, SZ_FNAME, TY_CHAR) + call salloc (projpars, SZ_LINE, TY_CHAR) + call salloc (wpars, SZ_LINE, TY_CHAR) + call salloc (r, ndim, TY_DOUBLE) + call salloc (w, ndim, TY_DOUBLE) + call salloc (cd, ndim * ndim, TY_DOUBLE) + call salloc (ltm, ndim * ndim, TY_DOUBLE) + call salloc (ltv, ndim, TY_DOUBLE) + call salloc (iltm, ndim * ndim, TY_DOUBLE) + call salloc (nr, ndim, TY_DOUBLE) + call salloc (ncd, ndim * ndim, TY_DOUBLE) + call salloc (axes, IM_MAXDIM, TY_INT) + call salloc (axno, IM_MAXDIM, TY_INT) + call salloc (axval, IM_MAXDIM, TY_INT) + + # Open the new wcs + mwnew = mw_open (NULL, ndim) + call mw_gsystem (mw, Memc[projstr], SZ_FNAME) + iferr { + call mw_newsystem (mw, "image", ndim) + } then { + call mw_newsystem (mwnew, Memc[projstr], ndim) + } else { + call mw_newsystem (mwnew, "image", ndim) + } + + # Set the LTERM. + call mw_gltermd (mw, Memd[ltm], Memd[ltv], ndim) + call mw_sltermd (mwnew, Memd[ltm], Memd[ltv], ndim) + + # Store the old axis map for later use. + call mw_gaxmap (mw, Memi[axno], Memi[axval], ndim) + + # Get the 2 logical axes. + call mw_gaxlist (mw, 03B, Memi[axes], naxes) + axmap = mw_stati (mw, MW_USEAXMAP) + ax1 = Memi[axes] + ax2 = Memi[axes+1] + + # Set the axes and projection type. + if (projection[1] == EOS) { + call mw_swtype (mwnew, Memi[axes], ndim, "linear", "") + } else { + call sscan (projection) + call gargwrd (Memc[projstr], SZ_FNAME) + call gargstr (Memc[projpars], SZ_LINE) + call sprintf (Memc[wpars], SZ_LINE, + "axis 1: axtype = ra %s axis 2: axtype = dec %s") + call pargstr (Memc[projpars]) + call pargstr (Memc[projpars]) + call mw_swtype (mwnew, Memi[axes], ndim, Memc[projstr], Memc[wpars]) + } + + # Copy in the atrributes of the other axes. + szatstr = SZ_LINE + call malloc (atstr, szatstr, TY_CHAR) + do l = 1, ndim { + if (l == ax1 || l == ax2) + next + iferr { + call mw_gwattrs (mw, l, "wtype", Memc[projpars], SZ_LINE) + } then { + call mw_swtype (mwnew, l, 1, "linear", "") + } else { + call mw_swtype (mwnew, l, 1, Memc[projpars], "") + } + for (i = 1; ; i = i + 1) { + if (itoc (i, Memc[projpars], SZ_LINE) <= 0) + Memc[atstr] = EOS + repeat { + iferr (call mw_gwattrs (mw, l, Memc[projpars], + Memc[atstr], szatstr)) + Memc[atstr] = EOS + if (strlen (Memc[atstr]) < szatstr) + break + szatstr = szatstr + SZ_LINE + call realloc (atstr, szatstr, TY_CHAR) + } + if (Memc[atstr] == EOS) + break + call mw_swattrs (mwnew, 1, Memc[projpars], Memc[atstr]) + } + } + call mfree (atstr, TY_CHAR) + + # Compute the referemce point world coordinates. + switch (sk_stati(coo, S_NLNGUNITS)) { + case SKY_DEGREES: + tlngref = lngref + case SKY_RADIANS: + tlngref = RADTODEG(lngref) + case SKY_HOURS: + tlngref = 15.0d0 * lngref + default: + tlngref = lngref + } + switch (sk_stati(coo, S_NLATUNITS)) { + case SKY_DEGREES: + tlatref = latref + case SKY_RADIANS: + tlatref = RADTODEG(latref) + case SKY_HOURS: + tlatref = 15.0d0 * latref + default: + tlatref = latref + } + + if (! transpose) { + Memd[w+ax1-1] = tlngref + Memd[w+ax2-1] = tlatref + } else { + Memd[w+ax2-1] = tlngref + Memd[w+ax1-1] = tlatref + } + + # Compute the reference point pixel coordinates. + Memd[nr+ax1-1] = xref + Memd[nr+ax2-1] = yref + + # Compute the new CD matrix. + if (! transpose) { + NEWCD(ax1,ax1) = xscale * cos (DEGTORAD(xrot)) / 3600.0d0 + NEWCD(ax2,ax1) = -yscale * sin (DEGTORAD(yrot)) / 3600.0d0 + NEWCD(ax1,ax2) = xscale * sin (DEGTORAD(xrot)) / 3600.0d0 + NEWCD(ax2,ax2) = yscale * cos (DEGTORAD(yrot)) / 3600.0d0 + } else { + NEWCD(ax1,ax1) = xscale * sin (DEGTORAD(xrot)) / 3600.0d0 + NEWCD(ax2,ax1) = yscale * cos (DEGTORAD(yrot)) / 3600.0d0 + NEWCD(ax1,ax2) = xscale * cos (DEGTORAD(xrot)) / 3600.0d0 + NEWCD(ax2,ax2) = -yscale * sin (DEGTORAD(yrot)) / 3600.0d0 + } + + if (! update) + call at_mwshow (mwnew, Memd[ltv], Memd[ltm], Memd[w], Memd[nr], + Memd[ncd], ndim) + + # Reset the axis map. + call mw_seti (mw, MW_USEAXMAP, axmap) + + # Recompute and store the new wcs if update is enabled. + call mw_saxmap (mwnew, Memi[axno], Memi[axval], ndim) + if (sk_stati (coo, S_PIXTYPE) == PIXTYPE_PHYSICAL) { + call mw_swtermd (mwnew, Memd[nr], Memd[w], Memd[ncd], ndim) + } else { + call mwmmuld (Memd[ncd], Memd[ltm], Memd[cd], ndim) + call mwinvertd (Memd[ltm], Memd[iltm], ndim) + call asubd (Memd[nr], Memd[ltv], Memd[r], ndim) + call mwvmuld (Memd[iltm], Memd[r], Memd[nr], ndim) + call mw_swtermd (mwnew, Memd[nr], Memd[w], Memd[cd], ndim) + } + + # Save the fit. + if (! transpose) { + call sk_seti (coo, S_PLNGAX, ax1) + call sk_seti (coo, S_PLATAX, ax2) + } else { + call sk_seti (coo, S_PLNGAX, ax2) + call sk_seti (coo, S_PLATAX, ax1) + } + if (update) { + call sk_saveim (coo, mwnew, im) + call mw_saveim (mwnew, im) + } + + # Close the wcs, + call mw_close (mwnew) + call mw_close (mw) + + # Force the CDELT keywords to update. This will be unecessary when + # mwcs is updated to deal with non-quoted and / or non left-justified + # CTYPE keywords.. + wtype = strdic (Memc[projstr], Memc[projstr], SZ_FNAME, WTYPE_LIST) + if (wtype > 0) + call sk_seti (coo, S_WTYPE, wtype) + call sk_ctypeim (coo, im) + + # Reset the fit. This will be unecessary when wcs is updated to deal + # with non-quoted and / or non left-justified CTYPE keywords. + call sk_seti (coo, S_WTYPE, 0) + call sk_seti (coo, S_PLNGAX, 0) + call sk_seti (coo, S_PLATAX, 0) + + call sfree (sp) +end + + +# AT_IMHMS -- Fetch a quantity form the image header that is in hms or dms +# format, e.g. in the form "+/-hh mm ss.x" or "+/-dd mm ss.s". + +double procedure at_imhms (im, kname) + +pointer im #I the image descriptor +char kname[ARB] #I the image keyword name + +double dval, hours, minutes, seconds +pointer sp, value +int nscan() +errchk imgstr() + +begin + call smark (sp) + call salloc (value, SZ_FNAME, TY_CHAR) + + iferr { + call imgstr (im, kname, Memc[value], SZ_FNAME) + } then { + dval = INDEFD + } else { + call sscan (Memc[value]) + call gargd (hours) + call gargd (minutes) + call gargd (seconds) + if (nscan() != 3) + dval = INDEFD + else if (hours >= 0.0d0) + dval = hours + (minutes / 60.0d0) + (seconds / 3600.0d0) + else + dval = -(abs(hours) + (minutes / 60.0d0) + (seconds / 3600.0d0)) + + } + + call sfree (sp) + + return (dval) +end diff --git a/noao/astcat/src/awcs/dcmpsv.f b/noao/astcat/src/awcs/dcmpsv.f new file mode 100644 index 00000000..5326e098 --- /dev/null +++ b/noao/astcat/src/awcs/dcmpsv.f @@ -0,0 +1,233 @@ +C This routine was copied from the stsdas$pkg/analysis/gasp/gasplib/ +C directory. See stsdas$copyright.stsdas for copyright restrictions. +C + 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/noao/astcat/src/awcs/dsswcs.x b/noao/astcat/src/awcs/dsswcs.x new file mode 100644 index 00000000..0fe3b169 --- /dev/null +++ b/noao/astcat/src/awcs/dsswcs.x @@ -0,0 +1,300 @@ +include +include +include + +define SZ_KEYWORD 8 +define SZ_PLATECOEFF 20 +define SZ_CDMATX 4 + +# AT_MKDSS -- Compute the FITS WCS from the general plate solution for a +# DSS image. This routine assumes that the geometry of the DSS image has not +# been modified since it was extracted, i.e. it has not been shifed,rotated, +# scaled, transposed etc. This routine has been adapted from one in the STSDAS +# GASP package, whose maim limitation for IRAF purposes was that it bypassed +# the IRAF MWCS routines. Return OK it the header is successfully updated, +# ERR otherwise. + +int procedure at_mkdss (im, update, verbose) + +pointer im #I the DSS image descriptor +bool update #I update rather than list the wcs ? +bool verbose #I verbose mode ? + +double amdx[SZ_PLATECOEFF] # the RA plate solution coefficients +double amdy[SZ_PLATECOEFF] # the DEC plate solution coefficients +double plate_cen_x # the x center position in microns +double plate_cen_y # the y center position in microns +double plate_cen_ra # the RA plate center in radians +double plate_cen_dec # the DEC plate center in radians +double x_pixel_size # the x step size in microns +double y_pixel_size # the y step size in microns +double plate_scale # the plate sclae in arcsec / mm +double im_x_center_pix # the x of ll corner of scanned plate +double im_y_center_pix # the y of ll corner of scanned plate +double ra_s # the plate center RA in seconds +double dec_s # the plate center DEC in seconds +double object_mag # the object magnitude +double object_col # the object color +int ra_h, ra_m # the plate center RA hours, minutes +int dec_d, dec_m # the plate center DEC degrees, minutes +char dec_sign # the plate center DEC sign +/- +int xcorner # the ll x of image w/r to plate +int ycorner # the ll y of image w/r to plate +int xsize # naxis1 +int ysize # naxis2 + +double crpix1, crpix2, crval1, crval2, cdmatx[SZ_CDMATX] +int i +char parname[SZ_KEYWORD] + +double imgetd() +real imgetr() +int imaccf(), imgeti() +errchk imgetr(), imgetd(), imgeti() + +begin + # Check that the image is 2D, if not it is not a DSS image. + if (IM_NDIM(im) != 2) + return (ERR) + + # See if image header contains the general plate solution. + if (imaccf (im,"PPO3 ") == NO) + return (ERR) + + # If we have an old DSS image, i.e. the one with the CRPIX rather + # than CNPIX keywords, rename CRPIX to CNPIX and proceed. + # this keyword to CNPIX and proceed. + + if (imaccf (im,"CRPIX1") == YES || imaccf (im, "CRPIX2") == YES) { + if (imaccf (im,"CRVAL1") == YES || imaccf (im, "CRVAL2") == YES) { + if (imaccf (im,"CD1_1") == NO && imaccf (im, "CD1_2") == NO && + imaccf (im, "CD2_1") == NO && imaccf (im, "CD2_2") == NO) { + # This is the case when we have CRPIX, CRVAL and no CD + # so, proceed to calculate the WCS again. + iferr (crpix1 = imgetr (im, "CRPIX1")) + return (ERR) + iferr (crpix2 = imgetr (im, "CRPIX2")) + return (ERR) + call imdelf (im, "CRPIX1") + call imaddr (im, "CNPIX1", real (crpix1)) + call imdelf (im, "CRPIX2") + call imaddr (im, "CNPIX2", real (crpix2)) + } + } else { + iferr (crpix1 = imgetr (im, "CRPIX1")) + return (ERR) + iferr (crpix2 = imgetr (im, "CRPIX2")) + return (ERR) + call imdelf (im, "CRPIX1") + call imaddr (im, "CNPIX1", real (crpix1)) + call imdelf (im, "CRPIX2") + call imaddr (im, "CNPIX2", real (crpix2)) + } + } + if (imaccf (im,"CNPIX1") == NO || imaccf (im, "CNPIX2") == NO ) + return (ERR) + + # Get the plate solution. + iferr { + + # Get the plate center parameters. + plate_cen_x = imgetd (im, "PPO3 ") + plate_cen_y = imgetd (im, "PPO6 ") + x_pixel_size = imgetd (im, "XPIXELSZ") + y_pixel_size = imgetd (im, "YPIXELSZ") + plate_scale = imgetd (im, "PLTSCALE") + ra_h = imgeti (im, "PLTRAH ") + ra_m = imgeti (im, "PLTRAM ") + ra_s = imgetd (im, "PLTRAS ") + call imgstr (im, "PLTDECSN", dec_sign, 1) + dec_d = imgeti (im, "PLTDECD ") + dec_m = imgeti (im, "PLTDECM ") + dec_s = imgetd (im, "PLTDECS ") + plate_cen_ra = DDEGTORAD ((ra_h + ra_m / 60.0d0 + ra_s / + 3600.0d0) * 15.0d0) + plate_cen_dec = DDEGTORAD (dec_d + dec_m / 60.0d0 + dec_s / + 3600.0d0) + if (dec_sign == '-') + plate_cen_dec = -plate_cen_dec + + # Get general plate solution coefficients + do i = 1, SZ_PLATECOEFF { + call sprintf (parname, SZ_KEYWORD, "AMDX%d") + call pargi(i) + amdx[i] = imgetd (im, parname) + } + do i = 1, SZ_PLATECOEFF { + call sprintf (parname, SZ_KEYWORD, "AMDY%d") + call pargi(i) + amdy[i] = imgetd (im, parname) + } + xcorner = imgetr (im, "CNPIX1") + ycorner = imgetr (im, "CNPIX2") + object_mag = 0.0d0 + object_col = 0.0d0 + } then + return (ERR) + + xsize = IM_LEN(im,1) + ysize = IM_LEN(im,2) + crpix1 = xsize / 2.0d0 + crpix2 = ysize / 2.0d0 + + # Center of image w/r to original lower left corner of scanned plate. + im_x_center_pix = xcorner + (xsize / 2.0d0) - 0.5d0 + im_y_center_pix = ycorner + (ysize / 2.0d0) - 0.5d0 + + # Calculate equatorial coordinates for the center of subset giving + # the complete plate solution w/r to the original lower left corner. + call ccgseq (plate_cen_ra, + plate_cen_dec, + plate_cen_x, + plate_cen_y, + x_pixel_size, + y_pixel_size, + plate_scale, + amdx, + amdy, + im_x_center_pix, + im_y_center_pix, + object_mag, + object_col, + crval1, + crval2) + + + # Calculate CD matrix values for the input subset from the original + # plate solution. + call calcds (plate_cen_ra, + plate_cen_dec, + plate_cen_x, + plate_cen_y, + xcorner, + ycorner, + x_pixel_size, + y_pixel_size, + plate_scale, + xsize, + ysize, + crval1, + crval2, + amdx, + amdy, + cdmatx) + + # Update the image header. + crval1 = DRADTODEG (crval1) + crval2 = DRADTODEG (crval2) + + if (verbose || ! update) + call printf (" Converting DSS wcs to FITS wcs\n") + call at_dmwcs (im, crpix1, crpix2, crval1, crval2, cdmatx, update) + + return (OK) +end + +define NEWCD Memd[ncd+(($2)-1)*ndim+($1)-1] + +# AT_DMWCS -- Create new image WCS from the approximation to the DSS plate +# solution. This routine assumes that the geometry of the DSS image has +# not been changed since since the image has been exracted from the image +# survey. + +procedure at_dmwcs (im, xref, yref, lngref, latref, cdmatx, update) + +pointer im #I pointer to the input image +double xref, yref #I the reference point in pixels +double lngref, latref #I the reference point in degrees +double cdmatx[ARB] #I CD1_1, CD1_2, CD2_1, CD2_2 +bool update #I update rather than list the wcs ? + +pointer mw, mwnew +pointer sp, projstr, r, w, cd, ltm, ltv, iltm, nr, ncd, axno, axval, axes +int ndim, ax1, ax2, naxes +pointer mw_openim(), mw_open() +int mw_stati() +errchk mw_newsystem() + +begin + mw = mw_openim (im) + ndim = mw_stati (mw, MW_NPHYSDIM) + + # Allocate working memory for the vectors and matrices. + call smark (sp) + call salloc (projstr, SZ_LINE, TY_CHAR) + call salloc (r, ndim, TY_DOUBLE) + call salloc (w, ndim, TY_DOUBLE) + call salloc (cd, ndim * ndim, TY_DOUBLE) + call salloc (ltm, ndim * ndim, TY_DOUBLE) + call salloc (ltv, ndim, TY_DOUBLE) + call salloc (iltm, ndim * ndim, TY_DOUBLE) + call salloc (nr, ndim, TY_DOUBLE) + call salloc (ncd, ndim * ndim, TY_DOUBLE) + call salloc (axno, IM_MAXDIM, TY_INT) + call salloc (axval, IM_MAXDIM, TY_INT) + call salloc (axes, IM_MAXDIM, TY_INT) + + # Open the new wcs. + mwnew = mw_open (NULL, ndim) + call mw_gsystem (mw, Memc[projstr], SZ_LINE) + iferr { + call mw_newsystem (mw, "image", ndim) + } then { + call mw_newsystem (mwnew, Memc[projstr], ndim) + } else { + call mw_newsystem (mwnew, "image", ndim) + } + + # Set the LTERM. + call mw_gltermd (mw, Memd[ltm], Memd[ltv], ndim) + call mw_sltermd (mwnew, Memd[ltm], Memd[ltv], ndim) + + # Store the old axis map for later use. + call mw_gaxmap (mw, Memi[axno], Memi[axval], ndim) + + # Get the 2 logical axes. + call mw_gaxlist (mw, 03B, Memi[axes], naxes) + ax1 = Memi[axes] + ax2 = Memi[axes+1] + + # Set the axes and projection type. + call sprintf (Memc[projstr], SZ_LINE, + "axis 1: axtype = ra axis 2: axtype = dec ") + call mw_swtype (mwnew, Memi[axes], ndim, "tan", Memc[projstr]) + + # Set the reference point world coordinates. + Memd[w+ax1-1] = lngref + Memd[w+ax2-1] = latref + + # Set the reference point pixel coordinates. + Memd[nr+ax1-1] = xref + Memd[nr+ax2-1] = yref + + # Compute the new CD matrix. + NEWCD(ax1,ax1) = cdmatx[1] # xscale * cos (DEGTORAD(xrot)) / 3600.0d0 + NEWCD(ax2,ax1) = cdmatx[2] # -yscale * sin (DEGTORAD(yrot)) / 3600.0d0 + NEWCD(ax1,ax2) = cdmatx[3] # xscale * sin (DEGTORAD(xrot)) / 3600.0d0 + NEWCD(ax2,ax2) = cdmatx[4] # yscale * cos (DEGTORAD(yrot)) / 3600.0d0 + + # List the new wcs. + if (! update) + call at_mwshow (mwnew, Memd[ltv], Memd[ltm], Memd[w], Memd[nr], + Memd[ncd], ndim) + + # Recompute and store the new wcs. + call mw_saxmap (mwnew, Memi[axno], Memi[axval], ndim) + call mwmmuld (Memd[ncd], Memd[ltm], Memd[cd], ndim) + call mwinvertd (Memd[ltm], Memd[iltm], ndim) + call asubd (Memd[nr], Memd[ltv], Memd[r], ndim) + call mwvmuld (Memd[iltm], Memd[r], Memd[nr], ndim) + call mw_swtermd (mwnew, Memd[nr], Memd[w], Memd[cd], ndim) + + # Update the image wcs. + if (update) + call mw_saveim (mwnew, im) + + call mw_close (mwnew) + call mw_close (mw) + + call sfree (sp) +end diff --git a/noao/astcat/src/awcs/fitsvd.f b/noao/astcat/src/awcs/fitsvd.f new file mode 100644 index 00000000..bb8d0f4e --- /dev/null +++ b/noao/astcat/src/awcs/fitsvd.f @@ -0,0 +1,38 @@ +C This routine was copied from the stsdas$pkg/analysis/gasp/gasplib/ +C directory. See stsdas$copyright.stsdas for copyright restrictions. +C + 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/noao/astcat/src/awcs/ksbsvd.f b/noao/astcat/src/awcs/ksbsvd.f new file mode 100644 index 00000000..3c78ec23 --- /dev/null +++ b/noao/astcat/src/awcs/ksbsvd.f @@ -0,0 +1,27 @@ +C This routines was copied from the stsdas$pkg/analysis/gasp/gasplib/ +C directory. See the file stsdas$copyright.stsdas for copyright +C restrictions. + 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/noao/astcat/src/awcs/mkpkg b/noao/astcat/src/awcs/mkpkg new file mode 100644 index 00000000..29160f2a --- /dev/null +++ b/noao/astcat/src/awcs/mkpkg @@ -0,0 +1,22 @@ +# AWCS task subdirectory + +$checkout libpkg.a "../" +$update libpkg.a +$checkin libpkg.a "../" +$exit + +libpkg.a: + dsswcs.x + dbwcs.x + parswcs.x "../../lib/astrom.h" "../../lib/aimpars.h" \ + + atmwshow.x + ccqseq.x + calcds.x + trsteq.x + treqst.x + dcmpsv.f + fitsvd.f + ksbsvd.f + varsvd.f + ; diff --git a/noao/astcat/src/awcs/parswcs.x b/noao/astcat/src/awcs/parswcs.x new file mode 100644 index 00000000..56cf33f2 --- /dev/null +++ b/noao/astcat/src/awcs/parswcs.x @@ -0,0 +1,251 @@ +include +include "../../lib/astrom.h" +include "../../lib/aimpars.h" +include + + +# AT_PARWCS -- Compute a FITS WCS from an image using WCS definitions +# read from the AWCSPARS parameter file and stored in the astromery +# package descriptor. At the moment I am going to keep this routine simple +# by not worrying about the units of any quantities but the world coordinates +# of the reference point. This routine can be made more sophisticated later +# as time permits. The information is there ... + +int procedure at_parwcs (im, at, update, verbose) + +pointer im #I the input image descriptor +pointer at #I the astrometry package descriptor +bool update #I update rather than list the wcs +bool verbose #I verbose mode ? + +double xref, yref, xmag, ymag, xrot, yrot, lngref, latref, dval +pointer sp, wfield, wtype, ctype, wcst, sym, coo, mw +int i, wkey, lngunits, latunits, coostat, stat +double at_imhms(), imgetd(), at_statd() +pointer at_statp(), stfind() +int at_wrdstr(), at_stati(), sk_decwcs() +bool streq() +errchk imgetd() + +begin + # Return if the input is not 2D. + if (IM_NDIM(im) != 2) + return (ERR) + + # Return if the wcs pointer is undefined. + if (at_statp (at, PWCS) == NULL) + return (ERR) + + # Return if the keyword symbol table is undefined. + wcst = at_statp (at, WCST) + if (wcst == NULL) + return (ERR) + + # Allocate working space. + call smark (sp) + call salloc (wfield, SZ_FNAME, TY_CHAR) + call salloc (wtype, SZ_FNAME, TY_CHAR) + call salloc (ctype, SZ_FNAME, TY_CHAR) + + # Initialize. + xref = (1.0d0 + IM_LEN(im,1)) / 2.0d0 + yref = (1.0d0 + IM_LEN(im,2)) / 2.0d0 + xmag = INDEFD + ymag = INDEFD + xrot = 180.0d0 + yrot= 0.0d0 + lngref = INDEFD + latref = INDEFD + lngunits = 0 + latunits = 0 + call strcpy ("tan", Memc[wtype], SZ_FNAME) + call strcpy ("J2000", Memc[ctype], SZ_FNAME) + + do i = 1, AT_NWFIELDS { + + # Which keyword have we got ? + wkey = at_wrdstr (i, Memc[wfield], SZ_FNAME, AT_WFIELDS) + + switch (wkey) { + + # Get the x reference point in pixels. + case WCS_WXREF: + sym = stfind (wcst, Memc[wfield]) + if (sym != NULL) { + if (streq (AT_WCSTKVAL(sym), "INDEF")) + dval = at_statd (at, WXREF) + else iferr (dval = imgetd (im, AT_WCSTKVAL(sym))) + dval = at_statd (at, WXREF) + } else + dval = at_statd (at, WXREF) + if (! IS_INDEFD(dval)) + xref = dval + + case WCS_WYREF: + sym = stfind (wcst, Memc[wfield]) + if (sym != NULL) { + if (streq (AT_WCSTKVAL(sym), "INDEF")) + dval = at_statd (at, WYREF) + else iferr (dval = imgetd (im, AT_WCSTKVAL(sym))) + dval = at_statd (at, WYREF) + } else + dval = at_statd (at, WYREF) + if (! IS_INDEFD(dval)) + yref = dval + + case WCS_WXMAG: + sym = stfind (wcst, Memc[wfield]) + if (sym != NULL) { + if (streq (AT_WCSTKVAL(sym), "INDEF")) + dval = at_statd (at, WXMAG) + else iferr (dval = imgetd (im, AT_WCSTKVAL(sym))) + dval = at_statd (at, WXMAG) + } else + dval = at_statd (at, WXMAG) + if (! IS_INDEFD(dval)) + xmag = dval + + case WCS_WYMAG: + sym = stfind (wcst, Memc[wfield]) + if (sym != NULL) { + if (streq (AT_WCSTKVAL(sym), "INDEF")) + dval = at_statd (at, WYMAG) + else iferr (dval = imgetd (im, AT_WCSTKVAL(sym))) + dval = at_statd (at, WYMAG) + } else + dval = at_statd (at, WYMAG) + if (! IS_INDEFD(dval)) + ymag = dval + + case WCS_WXROT: + sym = stfind (wcst, Memc[wfield]) + if (sym != NULL) { + if (streq (AT_WCSTKVAL(sym), "INDEF")) + dval = at_statd (at, WXROT) + else iferr (dval = imgetd (im, AT_WCSTKVAL(sym))) + dval = at_statd (at, WXROT) + } else + dval = at_statd (at, WXROT) + if (! IS_INDEFD(dval)) + xrot = dval + + case WCS_WYROT: + sym = stfind (wcst, Memc[wfield]) + if (sym != NULL) { + if (streq (AT_WCSTKVAL(sym), "INDEF")) + dval = at_statd (at, WYROT) + else iferr (dval = imgetd (im, AT_WCSTKVAL(sym))) + dval = at_statd (at, WYROT) + } else + dval = at_statd (at, WYROT) + if (! IS_INDEFD(dval)) + yrot = dval + + case WCS_WRAREF: + sym = stfind (wcst, Memc[wfield]) + if (sym != NULL) { + if (streq (AT_WCSTKVAL(sym), "INDEF")) + dval = at_statd (at, WRAREF) + else { + dval = at_imhms (im, AT_WCSTKVAL(sym)) + if (IS_INDEFD(dval)) { + iferr (dval = imgetd (im, AT_WCSTKVAL(sym))) + dval = at_statd (at, WRAREF) + } + } + } else + dval = at_statd (at, WRAREF) + if (! IS_INDEFD(dval)) + lngref = dval + + case WCS_WDECREF: + sym = stfind (wcst, Memc[wfield]) + if (sym != NULL) { + if (streq (AT_WCSTKVAL(sym), "INDEF")) + dval = at_statd (at, WDECREF) + else { + dval = at_imhms (im, AT_WCSTKVAL(sym)) + if (IS_INDEFD(dval)) { + iferr (dval = imgetd (im, AT_WCSTKVAL(sym))) + dval = at_statd (at, WDECREF) + } + } + } else + dval = at_statd (at, WDECREF) + if (! IS_INDEFD(dval)) + latref = dval + + case WCS_WRAUNITS: + lngunits = at_stati (at, WRAUNITS) + + case WCS_WDECUNITS: + latunits = at_stati (at, WDECUNITS) + + case WCS_WPROJ: + sym = stfind (wcst, Memc[wfield]) + if (sym != NULL) { + if (streq (AT_WCSTKVAL(sym), "INDEF")) + call at_stats (at, WPROJ, Memc[wtype], SZ_FNAME) + else iferr (call imgstr (im, AT_WCSTKVAL(sym), Memc[wtype], + SZ_FNAME)) + call at_stats (at, WPROJ, Memc[wtype], SZ_FNAME) + } else + call at_stats (at, WPROJ, Memc[wtype], SZ_FNAME) + if (streq (Memc[wtype], "INDEF")) + call strcpy ("tan", Memc[wtype], SZ_FNAME) + + case WCS_WSYSTEM: + sym = stfind (wcst, Memc[wfield]) + if (sym != NULL) { + if (streq (AT_WCSTKVAL(sym), "INDEF")) + call at_stats (at, WSYSTEM, Memc[ctype], SZ_FNAME) + else iferr (call imgstr (im, AT_WCSTKVAL(sym), Memc[ctype], + SZ_FNAME)) + call at_stats (at, WSYSTEM, Memc[ctype], SZ_FNAME) + } else + call at_stats (at, WSYSTEM, Memc[ctype], SZ_FNAME) + if (streq (Memc[ctype], "INDEF")) + call strcpy ("J2000", Memc[ctype], SZ_FNAME) + + default: + ; + } + } + + # Update the header. + if (IS_INDEFD(xmag) || IS_INDEFD(ymag) || IS_INDEFD(lngref) || + IS_INDEFD(latref)) { + + stat = ERR + + } else { + + # Open coordinate system struct + coostat = sk_decwcs (Memc[ctype], mw, coo, NULL) + + if (coostat == ERR || mw != NULL) { + if (mw != NULL) + call mw_close (mw) + stat = ERR + } else { + if (verbose) + call printf ( + " Writing FITS wcs using default parameters\n") + if (lngunits > 0) + call sk_seti (coo, S_NLNGUNITS, lngunits) + if (latunits > 0) + call sk_seti (coo, S_NLATUNITS, latunits) + call at_uwcs (im, coo, Memc[wtype], lngref, latref, xref, + yref, xmag, ymag, xrot, yrot, false, update) + stat = OK + } + + # Close the coordinate structure + if (coo != NULL) + call sk_close (coo) + } + + call sfree (sp) + + return (stat) +end diff --git a/noao/astcat/src/awcs/treqst.x b/noao/astcat/src/awcs/treqst.x new file mode 100644 index 00000000..c4cd27e5 --- /dev/null +++ b/noao/astcat/src/awcs/treqst.x @@ -0,0 +1,49 @@ +# This routine was copied from the stsdas$pkg/analysis/gasp/gasplib/ +# directory. See stsdas$copyright.stsdas for copyright restrictions. +# + +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 plate ra center (radians) +double plate_centre_dec #I plate dec center (radians) +double object_ra #I object ra center (radians) +double object_dec #I object dec center (radians) +double xi_object #O object xi standard coordinate (arcsecs) +double eta_object #O object eta standard coordinate (arcsecs) + +#double div +double ra, cosra, sinra, cosdec, sindec, cosd0, sind0, cosdist + +begin + ra = object_ra - plate_centre_ra + cosra = cos (ra) + sinra = sin (ra) + cosdec = cos (object_dec) + sindec = sin (object_dec) + cosd0 = cos (plate_centre_dec) + sind0 = sin (plate_centre_dec) + cosdist = sindec * sind0 + cosdec * cosd0 * cosra + xi_object = cosdec * sinra * ARCSEC_PER_RADIAN / cosdist + eta_object = (sindec * cosd0 - cosdec * sind0 * cosra) * + ARCSEC_PER_RADIAN / cosdist + +# # Find the divisor. +# div = (sin(object_dec) * sin(plate_centre_dec) + +# cos(object_dec) * cos(plate_centre_dec) * +# cos(object_ra -plate_centre_ra)) +# +# # Compute standard coords and convert to arcsec +# xi_object = cos(object_dec) * sin(object_ra-plate_centre_ra) * +# ARCSEC_PER_RADIAN/div +# eta_object = (sin(object_dec) * cos(plate_centre_dec) - +# cos(object_dec) * dsin(plate_centre_dec) * +# cos(object_ra - plate_centre_ra)) * +# ARCSEC_PER_RADIAN/div + +end diff --git a/noao/astcat/src/awcs/trsteq.x b/noao/astcat/src/awcs/trsteq.x new file mode 100644 index 00000000..5fa5ea8e --- /dev/null +++ b/noao/astcat/src/awcs/trsteq.x @@ -0,0 +1,64 @@ +# This routine was copied from stsdas$pkg/asnalysis/gasp/gasplib/. See +# stsdas$copyright.stsdas for copyright restrictions. +# +include + +define ARCSEC_PER_RADIAN 206264.8062470964d0 + +# TRSTEQ -- Procedure to compute the RA and DEC from the standard coordinates +# given the plate centre. + +procedure trsteq (plate_centre_ra, plate_centre_dec, xi, eta, ra, dec) + +double plate_centre_ra #I plate center ra (radians) +double plate_centre_dec #I plate center dec (radians) +double xi #I xi standard coordinate (arcsec) +double eta #I eta standard coordinate (arcsec) +double ra #O ra (radians) +double dec #O dec (radians) + +#double object_xi, object_eta, numerator, denominator +double object_xi, object_eta, x, y, z + +begin + # Convert from arcseconds to radians. + object_xi = xi/ARCSEC_PER_RADIAN + object_eta = eta/ARCSEC_PER_RADIAN + + # Convert to RA and Dec + x = cos (plate_centre_dec) - object_eta * sin (plate_centre_dec) + y = object_xi + z = sin (plate_centre_dec) + object_eta * cos (plate_centre_dec) + + if (x == 0.0d0 && y == 0.0d0) + ra = 0.0d0 + else + ra = atan2 (y, x) + dec = atan2 (z, sqrt (x * x + y * y)) + ra = ra + plate_centre_ra + if (ra < 0.0d0) + ra = ra + DTWOPI + else if (ra > DTWOPI) + ra = ra - DTWOPI + +## numerator = object_xi / dcos(plate_centre_dec) +# numerator = object_xi +# +## denominator = 1.0d0 - object_eta * dtan(plate_centre_dec) +# denominator = cos (plate_centre_dec) - +# object_eta * sin (plate_centre_dec) +# ra = atan2 (numerator,denominator) + plate_centre_ra +# if (ra < 0.0d0) +# ra = ra + DTWOPI +# else if (ra > DTWOPI) +# ra = ra - DTWOPI +# +## numerator = dcos(ra-plate_centre_ra) * +## (object_eta + dtan(plate_centre_dec)) +# numerator = cos (ra - plate_centre_ra) * +# (cos (plate_centre_dec) * object_eta + sin (plate_centre_dec)) +## denominator = 1.0d0 - object_eta * dtan(plate_centre_dec) +# denominator = cos (plate_centre_dec) - object_eta * +# sin (plate_centre_dec) +# dec = atan2 (numerator, denominator) +end diff --git a/noao/astcat/src/awcs/varsvd.f b/noao/astcat/src/awcs/varsvd.f new file mode 100644 index 00000000..b779a1e0 --- /dev/null +++ b/noao/astcat/src/awcs/varsvd.f @@ -0,0 +1,24 @@ +C This routine was copied from stsdas$pkg/analysis/gasp/gasplib/. +C See stsdas$copyright.stsdas for copyright restrictions. +C + 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 -- cgit