Initial commit

13 files changed, 1882 insertions, 0 deletions

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 <math.h>
+
+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 <imhdr.h>
+include <math.h>
+include <mwset.h>
+include <pkg/skywcs.h>
+include <pkg/cq.h>
+
+# 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 <imhdr.h>
+include <mwset.h>
+include <math.h>
+
+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	<imhdr.h> <mwset.h> <math.h>
+	dbwcs.x		<imhdr.h> <mwset.h> <math.h> <pkg/skywcs.h> <pkg/cq.h>
+	parswcs.x	<imhdr.h> "../../lib/astrom.h" "../../lib/aimpars.h" \
+	                <pkg/skywcs.h>
+	atmwshow.x
+	ccqseq.x
+	calcds.x	<math.h>
+	trsteq.x	<math.h>
+	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 <imhdr.h>
+include "../../lib/astrom.h"
+include "../../lib/aimpars.h"
+include <pkg/skywcs.h>
+
+
+# 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 <math.h>
+
+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