Repatch (from linux) of OSX IRAF

1 files changed, 1286 insertions, 0 deletions

diff --git a/pkg/utilities/nttools/stxtools/xtwcs.x b/pkg/utilities/nttools/stxtools/xtwcs.x
new file mode 100644
index 00000000..aa8b2798
--- /dev/null
+++ b/pkg/utilities/nttools/stxtools/xtwcs.x
@@ -0,0 +1,1286 @@
+include <imhdr.h>
+include <math.h>
+
+# This file contains the following high-level routines for converting
+# between world coordinates and pixel coordinates:
+#
+# xt_wcs_init		initialize struct for world coordinate system
+# xt_wcs_init_c		initialize from input cdelt, crota, etc
+# xt_wcs_init_cd	initialize from input CD matrix, etc
+# xt_wcs_free		deallocate wcs struct
+# xt_wc_pix		convert from world coordinates to pixel coordinates
+# xt_pix_wc		convert from pixel coordinates to world coordinates
+#
+# Phil Hodge, 27-Sept-1988  Created, based on code by Nelson & Zolt.
+# Phil Hodge, 6-April-1990  CD matrix mult. was transposed in xt_pix_wc.
+# Phil Hodge, 26-July-1991  In xt_e_ctype, change GBS to GLS (global sine).
+
+define	LEN_WCS		136	# size of wcs struct for naxis <= 7
+
+define	W_VALID		Memi[$1]	# coordinates valid, YES or NO?
+define	W_NAXIS		Memi[$1+1]	# number of axes
+define	W_RA_AX		Memi[$1+2]	# which axis is RA?  zero if none
+define	W_DEC_AX	Memi[$1+3]	# which axis is Dec?  zero if none
+define	W_PROJECTION	Memi[$1+4]	# projection type
+
+#   6 is currently not used
+
+#   7 -  55:  full CD matrix (7x7); units = e.g. degrees
+#  56 - 104:  LU decomposition of CD matrix
+# 105 - 111:  index returned by ludcmp for use by lubksb
+# 112 - 118:  reference pixel location
+# 119 - 122:  cosine & sine of declination at the reference pixel
+# 123 - 136:  coordinates at crpix; units = e.g. degrees
+
+define	W_CD		Memr[P2R($1+6 +($2-1)+($3-1)*7)]
+define	W_CDLU		Memr[P2R($1+55 +($2-1)+($3-1)*7)]
+define	W_CDINDX	Memr[P2R($1+104)]		# this is an array of 7
+define	W_CRPIX		Memr[P2R($1+110+$2)]
+define	W_COSDEC	Memd[P2D($1+118)]
+define	W_SINDEC	Memd[P2D($1+120)]
+define	W_CRVAL		Memd[P2D($1+120)+$2]
+
+# Projection types.
+
+define	W_LINEAR	0
+define	W_GNOMONIC	1	# TAN
+define	W_SINE		2	# SIN
+define	W_ARC		3	# ARC
+define	W_NORTH_POLAR	4	# NCP, north celestial pole (Westerbork)
+define	W_STEREOGRAPHIC	5	# STG (conformal)
+define	W_AITOFF	6	# AIT (equal-area)
+define	W_GLOBAL_SINE	7	# GLS (equal-area)
+define	W_MERCATOR	8	# MER (conformal)
+
+
+# xt_wcs_init -- initialize wcs struct
+# This routine allocates space for a structure describing the world
+# coordinate system for an image, fills in the values or defaults, and
+# returns a pointer to that structure.
+
+procedure xt_wcs_init (im, wcs)
+
+pointer im			# i: pointer to image descriptor
+pointer wcs			# o: pointer to world coord system struct
+#--
+real	dummy			# returned by ludcmp and ignored
+int	ira, idec		# index of RA, Dec axes
+int	j, k			# loop indexes
+errchk	xt_load_ctstruct
+
+begin
+	call calloc (wcs, LEN_WCS, TY_STRUCT)
+
+	W_VALID(wcs) = YES			# initial value
+	W_NAXIS(wcs) = IM_NDIM(im)
+
+	call xt_load_wcsstruct (im, wcs)	# get CRVAL, etc from image
+
+	if (W_NAXIS(wcs) >= 2) {
+
+	    ira = W_RA_AX(wcs)
+	    idec = W_DEC_AX(wcs)
+
+	    if (idec > 0) {
+		W_COSDEC(wcs) = cos (DEGTORAD(W_CRVAL(wcs,idec)))
+		W_SINDEC(wcs) = sin (DEGTORAD(W_CRVAL(wcs,idec)))
+	    } else {
+		W_COSDEC(wcs) = 1.d0
+		W_SINDEC(wcs) = 0.d0
+	    }
+
+	    # Copy the CD matrix to W_CDLU, and do the LU decomposition
+	    # on W_CDLU in-place.
+	    do k = 1, IM_MAXDIM
+		do j = 1, IM_MAXDIM
+		    W_CDLU(wcs,j,k) = W_CD(wcs,j,k)
+
+	    iferr {
+		call ludcmp (W_CDLU(wcs,1,1), W_NAXIS(wcs), IM_MAXDIM,
+			W_CDINDX(wcs), dummy)
+	    } then {
+		call mfree (wcs, TY_STRUCT)
+		call error (0, "xt_wcs_init:  cd matrix is singular")
+	    }
+	}
+end
+
+
+# xt_wcs_free -- deallocate wcs struct
+# This routine deallocates space for a wcs structure.
+
+procedure xt_wcs_free (wcs)
+
+pointer wcs		# io: pointer to world coord system struct
+#--
+
+begin
+	if (wcs != NULL)
+	    call mfree (wcs, TY_STRUCT)
+end
+
+
+# xt_wcs_init_c -- initialize wcs struct
+# xt_wcs_init_c and xt_wcs_init_cd allocate space for a structure
+# describing the world coordinate system for an image, fill in the values
+# or defaults, and return a pointer to that structure.  They differ from
+# xt_wcs_init in that these take the coordinate parameters as arguments
+# rather than getting them from the image.
+# xt_wcs_init_c takes cdelt & crota, and xt_wcs_init_cd takes the CD matrix.
+
+procedure xt_wcs_init_c (crval, crpix, cdelt, crota, ctype, naxis, wcs)
+
+double	crval[naxis]		# i: coordinate values at reference pixel
+real	crpix[naxis]		# i: reference pixel
+real	cdelt[naxis]		# i: pixel spacing
+real	crota			# i: rotation angle (if 2-D)
+char	ctype[SZ_CTYPE,naxis]	# i: e.g. "RA---TAN"
+int	naxis			# i: size of arrays
+pointer wcs			# o: pointer to world coord system struct
+#--
+real	dummy			# returned by ludcmp and ignored
+int	ira, idec		# index of RA, Dec axes
+int	j, k			# loop indexes
+errchk	ludcmp
+
+begin
+	do k = 1, naxis
+	    if (cdelt[k] == 0.)
+		call error (0, "xt_wcs_init_c:  zero value of CDELT")
+
+	call calloc (wcs, LEN_WCS, TY_STRUCT)
+
+	W_NAXIS(wcs) = naxis
+	W_VALID(wcs) = YES		# initial value
+
+	# Examine ctype to get ira, idec, proj_type.
+	call xt_e_ctype (ctype, naxis, ira, idec, W_PROJECTION(wcs))
+	W_RA_AX(wcs) = ira
+	W_DEC_AX(wcs) = idec
+
+	do k = 1, naxis {
+	    W_CRVAL(wcs,k) = crval[k]
+	    W_CRPIX(wcs,k) = crpix[k]
+	}
+	do k = naxis+1, IM_MAXDIM {
+	    W_CRVAL(wcs,k) = 0.d0
+	    W_CRPIX(wcs,k) = 1.
+	}
+
+	if (naxis == 1) {
+
+	    W_CD(wcs,1,1) = cdelt[1]
+
+	} else if (naxis >= 2) {
+
+	    if (idec > 0) {
+		W_COSDEC(wcs) = cos (DEGTORAD(W_CRVAL(wcs,idec)))
+		W_SINDEC(wcs) = sin (DEGTORAD(W_CRVAL(wcs,idec)))
+	    } else {
+		W_COSDEC(wcs) = 1.d0
+		W_SINDEC(wcs) = 0.d0
+	    }
+
+	    # Convert cdelt & crota to the CD matrix.
+	    call xt_to_cd (wcs, cdelt, crota, naxis)
+
+	    # Copy the CD matrix, and do the LU decomposition on W_CDLU.
+	    do k = 1, IM_MAXDIM
+		do j = 1, IM_MAXDIM
+		    W_CDLU(wcs,j,k) = W_CD(wcs,j,k)
+
+	    call ludcmp (W_CDLU(wcs,1,1), naxis, IM_MAXDIM,
+			W_CDINDX(wcs), dummy)
+	}
+end
+
+
+# xt_wcs_init_cd -- initialize wcs struct (CD)
+
+procedure xt_wcs_init_cd (crval, crpix, cd, ctype, naxis, wcs)
+
+double	crval[naxis]		# i: coordinate values at reference pixel
+real	crpix[naxis]		# i: reference pixel
+real	cd[naxis,naxis]		# i: CD matrix
+char	ctype[SZ_CTYPE,naxis]	# i: e.g. "RA---TAN"
+int	naxis			# i: size of arrays
+pointer wcs			# o: pointer to world coord system struct
+#--
+real	dummy			# returned by ludcmp and ignored
+int	ira, idec		# index of RA, Dec axes
+int	j, k			# loop indexes
+
+begin
+	call calloc (wcs, LEN_WCS, TY_STRUCT)
+
+	W_NAXIS(wcs) = naxis
+	W_VALID(wcs) = YES		# initial value
+
+	# Examine ctype to get ira, idec, proj_type.
+	call xt_e_ctype (ctype, naxis, ira, idec, W_PROJECTION(wcs))
+	W_RA_AX(wcs) = ira
+	W_DEC_AX(wcs) = idec
+
+	do k = 1, naxis {
+	    W_CRVAL(wcs,k) = crval[k]
+	    W_CRPIX(wcs,k) = crpix[k]
+	}
+	do k = naxis+1, IM_MAXDIM {
+	    W_CRVAL(wcs,k) = 0.d0
+	    W_CRPIX(wcs,k) = 1.
+	}
+
+	if (naxis == 1) {
+
+	    W_CD(wcs,1,1) = cd[1,1]
+
+	} else if (naxis >= 2) {
+
+	    if (idec > 0) {
+		W_COSDEC(wcs) = cos (DEGTORAD(W_CRVAL(wcs,idec)))
+		W_SINDEC(wcs) = sin (DEGTORAD(W_CRVAL(wcs,idec)))
+	    } else {
+		W_COSDEC(wcs) = 1.d0
+		W_SINDEC(wcs) = 0.d0
+	    }
+
+	    # Assign initial values to the CD matrix.
+	    do k = 1, IM_MAXDIM {
+		do j = 1, IM_MAXDIM {
+		    if (j == k) {
+			W_CD(wcs,k,k) = 1.
+			W_CDLU(wcs,k,k) = 1.
+		    } else {
+			W_CD(wcs,j,k) = 0.
+			W_CDLU(wcs,j,k) = 0.
+		    }
+		}
+	    }
+
+	    # Copy the CD matrix, and do the LU decomposition on W_CDLU.
+	    do k = 1, naxis {
+		do j = 1, naxis {
+		    W_CD(wcs,j,k) = cd[j,k]
+		    W_CDLU(wcs,j,k) = cd[j,k]
+		}
+	    }
+
+	    iferr {
+		call ludcmp (W_CDLU(wcs,1,1), naxis, IM_MAXDIM,
+			W_CDINDX(wcs), dummy)
+	    } then {
+		call mfree (wcs, TY_STRUCT)
+		call error (0, "xt_wcs_init_cd:  cd matrix is singular")
+	    }
+	}
+end
+
+# xt_to_cd -- from cdelt & crota to cd matrix
+# This routine computes the CD matrix from CDELT and CROTA.
+
+procedure xt_to_cd (wcs, cdelt, crota, naxis)
+
+pointer wcs			# i: pointer to world coord system struct
+real	cdelt[naxis]		# i: pixel spacing
+real	crota			# i: rotation angle (if 2-D)
+int	naxis			# i: size of arrays
+#--
+real	cosrota, sinrota	# cosine & sine of crota
+real	sign_cdelt[2]		# one, with sign of cdelt1 or cdelt2
+int	ira, idec		# index of RA, Dec axes
+int	j, k			# loop indexes
+
+begin
+	ira = W_RA_AX(wcs)
+	idec = W_DEC_AX(wcs)
+
+	if ( ! IS_INDEFD(crota) ) {
+	    cosrota = cos (DEGTORAD(crota))
+	    sinrota = sin (DEGTORAD(crota))
+	} else {
+	    cosrota = 1.d0
+	    sinrota = 0.d0
+	}
+
+	# Initial values for CD matrix.
+	do k = 1, IM_MAXDIM {
+	    do j = 1, IM_MAXDIM {
+		if (j == k)
+		    W_CD(wcs,k,k) = 1.
+		else
+		    W_CD(wcs,j,k) = 0.
+	    }
+	}
+	do k = 1, naxis
+	    W_CD(wcs,k,k) = cdelt[k]
+
+	if (ira > 0 && idec > 0) {
+
+	    if (cdelt[ira] >= 0.)
+		sign_cdelt[1] = 1.
+	    else
+		sign_cdelt[1] = -1.
+
+	    if (cdelt[idec] >= 0.)
+		sign_cdelt[2] = 1.
+	    else
+		sign_cdelt[2] = -1.
+
+	    W_CD(wcs,ira,ira) = cdelt[ira] * cosrota
+	    W_CD(wcs,ira,idec) = abs (cdelt[idec]) * sign_cdelt[1] * sinrota
+	    W_CD(wcs,idec,ira) = -abs (cdelt[ira]) * sign_cdelt[2] * sinrota
+	    W_CD(wcs,idec,idec) = cdelt[idec] * cosrota
+	}
+end
+
+# xt_e_ctype -- examine ctype
+# Examine each element of the ctype array to find which axes (if any)
+# are RA & Dec (or glon & glat, etc).  Also get the projection type,
+# such as gnomonic, if this was specified in ctype.
+
+procedure xt_e_ctype (ctype, naxis, ra_axis, dec_axis, proj_type)
+
+char	ctype[SZ_CTYPE,naxis]	# i: coordinate type, e.g. "RA---TAN"
+int	naxis			# i: dimension
+int	ra_axis			# o: which axis is RA (or glon, etc)?
+int	dec_axis		# o: which axis is Dec (or glat, etc)?
+int	proj_type		# o: type of projection
+#--
+char	lctype[SZ_CTYPE]	# local copy of an element of ctype
+char	dash			# '-'
+int	k
+int	index			# index of '-' in ctype
+int	strncmp(), strldx()
+
+begin
+	# Assign defaults.
+	ra_axis = 0
+	dec_axis = 0
+	if (naxis == 1)
+	    proj_type = W_LINEAR
+	else
+	    proj_type = W_GNOMONIC
+
+	# Search for "RA", "DEC", etc.
+	do k = 1, naxis {
+	    # Make a local copy of ctype & make sure it's upper case.
+	    call strcpy (ctype[1,k], lctype, SZ_CTYPE)
+	    call strupr (lctype)
+
+	    if (strncmp (lctype, "RA", 2) == 0)
+		ra_axis = k
+	    else if (strncmp (lctype, "DEC", 3) == 0)
+		dec_axis = k
+
+	    else if (strncmp (lctype, "GLON", 4) == 0)
+		ra_axis = k
+	    else if (strncmp (lctype, "LL", 2) == 0)
+		ra_axis = k
+	    else if (strncmp (lctype, "UU", 2) == 0)
+		ra_axis = k
+	    else if (strncmp (lctype, "ELON", 4) == 0)
+		ra_axis = k
+
+	    else if (strncmp (lctype, "GLAT", 4) == 0)
+		dec_axis = k
+	    else if (strncmp (lctype, "MM", 2) == 0)
+		dec_axis = k
+	    else if (strncmp (lctype, "VV", 2) == 0)
+		dec_axis = k
+	    else if (strncmp (lctype, "ELAT", 4) == 0)
+		dec_axis = k
+	}
+
+	if (ra_axis > 0)
+	    k = ra_axis
+	else if (dec_axis > 0)
+	    k = dec_axis
+	else
+	    k = 0
+
+	# If at least one of the axes is like RA or Dec, check to see
+	# whether a projection type was specified.
+	if (k > 0) {
+	    dash = '-'
+	    index = strldx (dash, lctype)
+	    if (index > 0) {
+		index = index + 1
+		if (strncmp (lctype[index], "TAN", 3) == 0)
+		    proj_type = W_GNOMONIC
+		else if (strncmp (lctype[index], "SIN", 3) == 0)
+		    proj_type = W_SINE
+		else if (strncmp (lctype[index], "ARC", 3) == 0)
+		    proj_type = W_ARC
+		else if (strncmp (lctype[index], "NCP", 3) == 0)
+		    proj_type = W_NORTH_POLAR
+		else if (strncmp (lctype[index], "STG", 3) == 0)
+		    proj_type = W_STEREOGRAPHIC
+		else if (strncmp (lctype[index], "AIT", 3) == 0)
+		    proj_type = W_AITOFF
+		else if (strncmp (lctype[index], "GLS", 3) == 0)
+		    proj_type = W_GLOBAL_SINE
+		else if (strncmp (lctype[index], "MER", 3) == 0)
+		    proj_type = W_MERCATOR
+	    }
+	}
+end
+
+
+define	SZ_PNAME	8
+
+# xt_load_wcsstruct -- load coordinate information
+# Get the coordinate information from the image, and load
+# that info into the wcs structure.
+
+procedure xt_load_wcsstruct (im, wcs)
+
+pointer im		# i: pointer to image header struct
+pointer wcs		# i: pointer to world coord system struct
+#--
+char	pname[SZ_PNAME]
+char	ctype[SZ_CTYPE,IM_MAXDIM]
+int	naxis, iax		# dimension of image; loop index for axis
+bool	cdm_found		# true if CD matrix present in image
+int	imaccf()
+double	imgetd()
+real	imgetr()
+errchk	imgstr, imgetd, imgetr, xt_g_cd_matrix, xt_c_cd_matrix
+
+begin
+	naxis = IM_NDIM(im)
+
+	# Get the coordinate info.  If anything is missing set W_VALID to NO.
+	do iax = 1, naxis {
+
+	    # CTYPE for each axis.
+	    call sprintf (pname, SZ_PNAME, "ctype%d")
+		call pargi (iax)
+	    if (imaccf (im, pname) == YES) {
+		call imgstr (im, pname, ctype[1,iax], SZ_CTYPE)
+	    } else {
+		call strcpy ("PIXEL", ctype[1,iax], SZ_CTYPE)
+		W_VALID(wcs) = NO
+	    }
+
+	    # CRVAL for each axis
+	    call sprintf (pname, SZ_PNAME, "crval%d")
+		call pargi (iax)
+	    if (imaccf (im, pname) == YES) {
+		W_CRVAL(wcs,iax) = imgetd (im, pname)
+	    } else {
+		W_CRVAL(wcs,iax) = 0.d0
+		W_VALID(wcs) = NO
+	    }
+
+	    # CRPIX for each axis
+	    call sprintf (pname, SZ_PNAME, "crpix%d")
+		call pargi (iax)
+	    if (imaccf (im, pname) == YES) {
+		W_CRPIX(wcs,iax) = imgetr (im, pname)
+	    } else {
+		W_CRPIX(wcs,iax) = 1.
+		W_VALID(wcs) = NO
+	    }
+	}
+	# Assign reasonable values to the unused elements.
+	do iax = naxis+1, IM_MAXDIM {
+	    W_CRVAL(wcs,iax) = 0.d0
+	    W_CRPIX(wcs,iax) = 1.
+	}
+
+	# Examine ctype array.
+	call xt_e_ctype (ctype, naxis,
+		W_RA_AX(wcs), W_DEC_AX(wcs), W_PROJECTION(wcs))
+
+	# First try to get the CD matrix, and if it isn't there
+	# get CDELT and CROTA and convert to CD.
+
+	call xt_g_cd_matrix (im, wcs, naxis, cdm_found)
+
+	if ( ! cdm_found )
+	    call xt_c_cd_matrix (im, wcs, naxis)
+end
+
+
+# xt_g_cd_matrix -- get CD matrix
+# If the CD matrix is present, get the values and place them into the
+# wcs structure.  Note that we assume that if *any* of the CD matrix
+# parameters are there, they are *all* there.
+
+define	TOLER	1.e-5
+
+procedure xt_g_cd_matrix (im, wcs, naxis, cdm_found)
+
+pointer	im			# i: image pointer
+pointer wcs			# i: pointer to wcs structure
+int	naxis			# i: number of axes in image
+bool	cdm_found		# o: true if CD matrix found
+#--
+real	cd_matrix[IM_MAXDIM,IM_MAXDIM]	# the CD matrix
+char	pname[SZ_PNAME]
+int	i, j
+int	imaccf()
+real	imgetr()
+errchk	imgetr
+
+begin
+	# This is reset below if any element of the CD matrix is found.
+	cdm_found = false
+
+	# Assign default values.
+	do j = 1, IM_MAXDIM
+	    do i= 1, IM_MAXDIM
+		if (i == j)
+		    cd_matrix[i,j] = 1.
+		else
+		    cd_matrix[i,j] = 0.
+
+	# Get each element of the CD matrix.
+	do j = 1, naxis {
+	    do i = 1, naxis {
+		call sprintf (pname, SZ_PNAME, "cd%d_%d")
+		    call pargi (i)
+		    call pargi (j)
+		if (imaccf (im, pname) == YES) {
+		    cd_matrix[i,j] = imgetr (im, pname)
+		    cdm_found = true
+		}
+	    }
+	}
+
+	# Copy to the wcs structure.
+	do j = 1, IM_MAXDIM
+	    do i = 1, IM_MAXDIM
+		W_CD(wcs,i,j) = cd_matrix[i,j]
+end
+
+
+# xt_c_cd_matrix -- create CD matrix
+# If the CD matrix is not present, get the values of CDELT & CROTA,
+# convert to the CD matrix, and store the values in the wcs structure.
+# Since this is called after trying unsuccessfully to get the CD matrix,
+# if cdelt or crota is not present W_VALID will be reset to NO.
+
+procedure xt_c_cd_matrix (im, wcs, naxis)
+
+pointer	im			# i: image pointer
+pointer wcs			# i: pointer to wcs structure
+int	naxis			# i: number of axes in image
+#--
+char	pname[SZ_PNAME]		# parameter name (e.g. "cdelt1")
+real	cdelt[IM_MAXDIM]	# pixel spacing
+real	crota			# rotation angle in degrees
+int	k			# loop index for axis
+int	imaccf()
+real	imgetr()
+errchk	imgetr
+
+begin
+	do k = 1, naxis {
+
+	    # CDELT for each axis.
+	    call sprintf (pname, SZ_PNAME, "cdelt%d")
+		call pargi (k)
+	    if (imaccf (im, pname) == YES) {
+		cdelt[k] = imgetr (im, pname)
+		if (cdelt[k] == 0.)
+		    call error (0, "xt_c_cd_matrix:  cdelt is zero")
+	    } else {
+		cdelt[k] = 1.
+		W_VALID(wcs) = NO
+	    }
+	}
+
+	# For a 1-D image, assign CD1_1 and return.
+	if (naxis == 1) {
+	    W_CD(wcs,1,1) = cdelt[1]
+	    return
+	}
+
+	# CROTA (only one).
+	call strcpy ("crota1", pname, SZ_PNAME)
+	if (imaccf (im, pname) == YES) {
+	    crota = imgetr (im, pname)
+	} else {
+	    crota = 0.
+	    W_VALID(wcs) = NO
+	}
+
+	# Compute CD matrix from CDELT & CROTA.
+	call xt_to_cd (wcs, cdelt, crota, naxis)
+end
+
+
+# xt_wc_pix -- wcs to pixels
+# This routine converts world coordinates to pixel coordinates.
+#
+# In the 1-D case, CRVAL is subtracted from the coordinate, the
+# result is divided by CDELT (same as CD1_1), and CRPIX is added.
+#
+# For 2-D or higher dimension, if two of the axes are like RA and Dec,
+# the input coordinates are converted to standard coordinates Xi
+# and Eta.  The (Xi, Eta) vector is then multiplied on the left by
+# the inverse of the CD matrix, and CRPIX is added.
+# The units for axes like Ra & Dec are degrees, not hours or radians.
+# For linear axes the conversion is the same as for 1-D.
+
+procedure xt_wc_pix (wcs, phys, pix, naxis)
+
+pointer wcs		# i: pointer to world coord system struct
+double	phys[naxis]	# i: physical (world) coordinates (e.g. degrees)
+real	pix[naxis]	# o: pixel coordinates
+int	naxis		# i: size of arrays
+#--
+double	delta_ra		# RA of object - RA at reference pixel
+double	dra_r, dec_r		# delta_ra & declination in radians
+double	xi_r, eta_r		# xi & eta in radians
+real	dphys[IM_MAXDIM]	# phys coord - reference coord
+int	ira, idec		# index of RA, Dec axes
+int	k			# loop index
+errchk	xt_wp_ncp, xt_wp_mer
+
+begin
+	do k = 1, naxis
+	    dphys[k] = phys[k] - W_CRVAL(wcs,k)
+
+	if (naxis == 1) {
+
+	    pix[1] = dphys[1] / W_CD(wcs,1,1) + W_CRPIX(wcs,1)
+
+	} else {
+
+	    ira = W_RA_AX(wcs)
+	    idec = W_DEC_AX(wcs)
+
+	    # Convert RA & Dec to Xi & Eta (standard coordinates).
+	    if (ira > 0 && idec > 0) {
+
+		delta_ra = phys[ira] - W_CRVAL(wcs,ira)		# double prec
+		dra_r = DEGTORAD (delta_ra)
+		dec_r = DEGTORAD (phys[idec])
+
+		switch (W_PROJECTION(wcs)) {
+		case W_GNOMONIC:
+		    call xt_wp_tan (wcs, dra_r, dec_r, xi_r, eta_r)
+		case W_SINE:
+		    call xt_wp_sin (wcs, dra_r, dec_r, xi_r, eta_r)
+		case W_ARC:
+		    call xt_wp_arc (wcs, dra_r, dec_r, xi_r, eta_r)
+		case W_NORTH_POLAR:
+		    call xt_wp_ncp (wcs, dra_r, dec_r, xi_r, eta_r)
+		case W_STEREOGRAPHIC:
+		    call xt_wp_stg (wcs, dra_r, dec_r, xi_r, eta_r)
+		case W_AITOFF:
+		    call xt_wp_ait (wcs, dra_r, dec_r, xi_r, eta_r)
+		case W_GLOBAL_SINE:
+		    call xt_wp_gls (wcs, dra_r, dec_r, xi_r, eta_r)
+		case W_MERCATOR:
+		    call xt_wp_mer (wcs, dra_r, dec_r, xi_r, eta_r)
+		}
+
+		dphys[ira] = RADTODEG (xi_r)		# xi, eta in degrees
+		dphys[idec] = RADTODEG (eta_r)
+	    }
+
+	    # Use LU backsubstitution to get pixel coords from physical coords.
+	    call lubksb (W_CDLU(wcs,1,1), naxis, IM_MAXDIM,
+			W_CDINDX(wcs), dphys)	# dphys is modified in-place
+	    do k = 1, naxis
+		pix[k] = dphys[k] + W_CRPIX(wcs,k)	# copy to output
+	}
+end
+
+
+# xt_pix_wc -- pixels to wcs
+# This routine converts pixel coordinates to world coordinates.
+#
+# In the 1-D case, CRPIX is subtracted from the pixel coordinate,
+# the result is multiplied by CDELT (same as CD1_1), and CRVAL is added.
+#
+# For 2-D or higher dimension, CRPIX is subtracted, and the result is
+# multiplied on the left by the CD matrix.  If two of the axes are like
+# RA and Dec, the pixel coordinates are converted to standard coordinates
+# Xi and Eta.  The (xi, eta) vector is then converted to differences
+# between RA and Dec and CRVAL, and then CRVAL is added to each coordinate.
+
+procedure xt_pix_wc (wcs, pix, phys, naxis)
+
+pointer wcs		# i: pointer to world coord system struct
+real	pix[naxis]	# i: pixel coordinates
+double	phys[naxis]	# o: physical (world) coordinates
+int	naxis		# i: size of arrays
+#--
+double	dpix[IM_MAXDIM]		# pix coord - crpix
+double	sum			# for matrix multiplication
+double	dra_r, dec_r		# delta_ra & declination in radians
+double	xi_r, eta_r		# xi & eta in radians
+int	ira, idec		# index of RA, Dec axes
+int	j, k			# loop indexes
+
+begin
+	do k = 1, naxis
+	    dpix[k] = pix[k] - W_CRPIX(wcs,k)
+
+	if (naxis == 1) {
+
+	    phys[1] = dpix[1] * W_CD(wcs,1,1) + W_CRVAL(wcs,1)
+
+	} else {
+
+	    do j = 1, naxis {
+		sum = 0.d0
+		do k = 1, naxis
+		    sum = sum + W_CD(wcs,j,k) * dpix[k]
+		phys[j] = sum
+	    }
+
+	    ira = W_RA_AX(wcs)
+	    idec = W_DEC_AX(wcs)
+
+	    # Convert Xi & Eta (standard coordinates) to RA & Dec.
+	    if (ira > 0 && idec > 0) {
+		xi_r = DEGTORAD (phys[ira])
+		eta_r = DEGTORAD (phys[idec])
+
+		switch (W_PROJECTION(wcs)) {
+		case W_GNOMONIC:
+		    call xt_pw_tan (wcs, xi_r, eta_r, dra_r, dec_r)
+		case W_SINE:
+		    call xt_pw_sin (wcs, xi_r, eta_r, dra_r, dec_r)
+		case W_ARC:
+		    call xt_pw_arc (wcs, xi_r, eta_r, dra_r, dec_r)
+		case W_NORTH_POLAR:
+		    call xt_pw_ncp (wcs, xi_r, eta_r, dra_r, dec_r)
+		case W_STEREOGRAPHIC:
+		    call xt_pw_stg (wcs, xi_r, eta_r, dra_r, dec_r)
+		case W_AITOFF:
+		    call xt_pw_ait (wcs, xi_r, eta_r, dra_r, dec_r)
+		case W_GLOBAL_SINE:
+		    call xt_pw_gls (wcs, xi_r, eta_r, dra_r, dec_r)
+		case W_MERCATOR:
+		    call xt_pw_mer (wcs, xi_r, eta_r, dra_r, dec_r)
+		}
+
+		phys[idec] = RADTODEG (dec_r)
+		phys[ira] = RADTODEG (dra_r) + W_CRVAL(wcs,ira)
+		if (phys[ira] < 0.d0)
+		    phys[ira] = phys[ira] + 360.d0
+	    }
+	    do k = 1, naxis
+		if (k != ira && k != idec)
+		    phys[k] = phys[k] + W_CRVAL(wcs,k)
+	}
+end
+
+
+# xt_wp_tan -- convert from ra & dec using gnomonic projection
+
+procedure xt_wp_tan (wcs, dra_r, dec_r, xi_r, eta_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	dra_r		# i: RA of object - RA at reference pixel (radians)
+double	dec_r		# i: declination of object (radians)
+double	xi_r		# o: standard coordinate (radians)
+double	eta_r		# o: standard coordinate (radians)
+#--
+double	cosdra, sindra	# cos & sin of dra_r
+double	cosdec, sindec	# cos & sin of object declination
+double	cosdist		# cos of dist from ref pixel to object
+
+begin
+	cosdra = cos (dra_r)
+	sindra = sin (dra_r)
+
+	cosdec = cos (dec_r)
+	sindec = sin (dec_r)
+
+	cosdist = sindec * W_SINDEC(wcs) + cosdec * W_COSDEC(wcs) * cosdra
+
+	xi_r = cosdec * sindra / cosdist
+	eta_r = (sindec * W_COSDEC(wcs) -
+			cosdec * W_SINDEC(wcs) * cosdra) / cosdist
+end
+
+
+# xt_pw_tan -- convert to ra & dec using gnomonic projection
+# In rectangular coordinates the vector (1, xi, eta) points toward
+# the object; the origin is the observer's location, the x-axis points
+# toward the reference pixel, the y-axis is in the direction of increasing
+# right ascension, and the z-axis is in the direction of increasing
+# declination.  The coordinate system is then rotated by the declination so
+# the x-axis passes through the equator at the RA of the reference pixel;
+# the components of the vector in this coordinate system are used to
+# compute (RA - reference_RA) and declination.
+
+procedure xt_pw_tan (wcs, xi_r, eta_r, dra_r, dec_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	xi_r		# i: standard coordinate (radians)
+double	eta_r		# i: standard coordinate (radians)
+double	dra_r		# o: RA of object - RA at reference pixel (radians)
+double	dec_r		# o: declination of object (radians)
+#--
+double	x, y, z		# vector (not unit length) pointing toward object
+
+begin
+	# Rotate the rectangular coordinate system of the vector (1, xi, eta)
+	# by the declination so the x-axis will pass through the equator.
+	x = W_COSDEC(wcs) - eta_r * W_SINDEC(wcs)
+	y = xi_r
+	z = W_SINDEC(wcs) + eta_r * W_COSDEC(wcs)
+
+	if (x == 0.d0 && y == 0.d0)
+	    dra_r = 0.d0
+	else
+	    dra_r = atan2 (y, x)
+	dec_r = atan2 (z, sqrt (x*x + y*y))
+end
+
+
+# xt_wp_sin -- convert from ra & dec using sine projection
+#
+# Reference:  AIPS Memo No. 27 by Eric W. Greisen
+
+procedure xt_wp_sin (wcs, dra_r, dec_r, xi_r, eta_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	dra_r		# i: RA of object - RA at reference pixel (radians)
+double	dec_r		# i: declination of object (radians)
+double	xi_r		# o: standard coordinate (radians)
+double	eta_r		# o: standard coordinate (radians)
+#--
+double	cosdra, sindra	# cos & sin of delta_ra
+double	cosdec, sindec	# cos & sin of object declination
+
+begin
+	cosdra = cos (dra_r)
+	sindra = sin (dra_r)
+
+	cosdec = cos (dec_r)
+	sindec = sin (dec_r)
+
+	xi_r = cosdec * sindra
+	eta_r = sindec * W_COSDEC(wcs) - cosdec * W_SINDEC(wcs) * cosdra
+end
+
+
+# xt_pw_sin -- convert to ra & dec using sine projection
+# In rectangular coordinates the vector (v1, xi, eta), where
+# v1 = sqrt (1 - xi**2 - eta**2), is the location of the object on the
+# unit celestial sphere.  The x-axis points toward the reference pixel,
+# the y-axis is in the direction of increasing right ascension, and the
+# z-axis is in the direction of increasing declination.  The coordinate
+# system is then rotated (around the y-axis) by the declination so the
+# x-axis passes through the equator at the RA of the reference pixel;
+# the components of the vector in this coordinate system are used to
+# compute (RA - reference_RA) and declination.
+
+procedure xt_pw_sin (wcs, xi_r, eta_r, dra_r, dec_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	xi_r		# i: standard coordinate (radians)
+double	eta_r		# i: standard coordinate (radians)
+double	dra_r		# o: RA of object - RA at reference pixel (radians)
+double	dec_r		# o: declination of object (radians)
+#--
+double	v1		# x component of unit vector
+double	x, y, z		# unit vector with x[1] pointing toward equator
+
+begin
+	v1 = sqrt (1.d0 - xi_r*xi_r - eta_r*eta_r)
+
+	# Rotate the rectangular coordinate system of the vector (v1, xi, eta)
+	# by the declination so the x-axis will pass through the equator.
+	x = v1 * W_COSDEC(wcs) - eta_r * W_SINDEC(wcs)
+	y = xi_r
+	z = v1 * W_SINDEC(wcs) + eta_r * W_COSDEC(wcs)
+
+	if (x == 0.d0 && y == 0.d0)
+	    dra_r = 0.d0
+	else
+	    dra_r = atan2 (y, x)
+	dec_r = atan2 (z, sqrt (x*x + y*y))
+end
+
+
+# xt_wp_arc -- convert from ra & dec using arc projection
+#
+# Reference:  AIPS Memo No. 27 by Eric W. Greisen
+
+procedure xt_wp_arc (wcs, dra_r, dec_r, xi_r, eta_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	dra_r		# i: RA of object - RA at reference pixel (radians)
+double	dec_r		# i: declination of object (radians)
+double	xi_r		# o: standard coordinate (radians)
+double	eta_r		# o: standard coordinate (radians)
+#--
+double	cosdra, sindra	# cos & sin of delta_ra
+double	cosdec, sindec	# cos & sin of object declination
+double	theta		# distance (radians) from ref pixel to object
+double	r		# theta / sin (theta)
+
+begin
+	cosdra = cos (dra_r)
+	sindra = sin (dra_r)
+
+	cosdec = cos (dec_r)
+	sindec = sin (dec_r)
+
+	theta = acos (sindec * W_SINDEC(wcs) + cosdec * W_COSDEC(wcs) * cosdra)
+	if (theta == 0.d0)
+	    r = 1.d0
+	else
+	    r = theta / sin (theta)
+
+	xi_r = r * cosdec * sindra
+	eta_r = r * (sindec * W_COSDEC(wcs) - cosdec * W_SINDEC(wcs) * cosdra)
+end
+
+
+# xt_pw_arc -- convert to ra & dec using arc projection
+# The rectangular coordinates of the pixel on a unit celestial sphere
+# are computed in a coordinate system such that the x-axis points toward
+# the reference pixel, the y-axis is in the direction of increasing right
+# ascension, and the z-axis is in the direction of increasing declination.
+# The coordinate system is then rotated (around the y-axis) by the
+# declination so the x-axis passes through the equator at the RA of the
+# reference pixel; the components of the vector in this coordinate system
+# are used to compute (RA - reference_RA) and declination.
+
+procedure xt_pw_arc (wcs, xi_r, eta_r, dra_r, dec_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	xi_r		# i: standard coordinate (radians)
+double	eta_r		# i: standard coordinate (radians)
+double	dra_r		# o: RA of object - RA at reference pixel (radians)
+double	dec_r		# o: declination of object (radians)
+#--
+double	theta		# arc length, i.e. sqrt (xi**2 + eta**2)
+double	v[3]		# unit vector with v[1] pointing toward ref pixel
+double	x, y, z		# vector with x[1] pointing toward equator
+
+begin
+	theta = sqrt (xi_r*xi_r + eta_r*eta_r)
+	if (theta == 0.d0) {
+	    v[1] = 1.d0
+	    v[2] = 0.d0
+	    v[3] = 0.d0
+	} else {
+	    v[1] = cos (theta)
+	    v[2] = sin (theta) / theta * xi_r
+	    v[3] = sin (theta) / theta * eta_r
+	}
+
+	# Rotate the rectangular coordinate system of the vector v by the
+	# declination so the x-axis will pass through the equator.
+	x = v[1] * W_COSDEC(wcs) - v[3] * W_SINDEC(wcs)
+	y = v[2]
+	z = v[1] * W_SINDEC(wcs) + v[3] * W_COSDEC(wcs)
+
+	if (x == 0.d0 && y == 0.d0)
+	    dra_r = 0.d0
+	else
+	    dra_r = atan2 (y, x)
+	dec_r = atan2 (z, sqrt (x*x + y*y))
+end
+
+
+# xt_wp_ncp -- convert from ra & dec using ncp projection
+#
+# References:
+#	AIPS Memo No. 27 by Eric W. Greisen
+#	Data Processing for the Westerbork Synthesis Radio Telescope
+#		by W. N. Brouw
+
+procedure xt_wp_ncp (wcs, dra_r, dec_r, xi_r, eta_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	dra_r		# i: RA of object - RA at reference pixel (radians)
+double	dec_r		# i: declination of object (radians)
+double	xi_r		# o: standard coordinate (radians)
+double	eta_r		# o: standard coordinate (radians)
+#--
+double	cosdra, sindra	# cos & sin of delta_ra
+double	cosdec		# cos of object declination
+
+begin
+	if (W_SINDEC(wcs) == 0.)
+	    call error (1, "NCP projection:  dec is zero")
+
+	cosdra = cos (dra_r)
+	sindra = sin (dra_r)
+
+	cosdec = cos (dec_r)
+
+	xi_r = - cosdec * sindra
+	eta_r = (W_COSDEC(wcs) - cosdec * cosdra) / W_SINDEC(wcs)
+end
+
+
+# xt_pw_ncp -- convert to ra & dec using ncp projection
+#
+# References:
+#	AIPS Memo No. 27 by Eric W. Greisen
+#	Data Processing for the Westerbork Synthesis Radio Telescope
+#		by W. N. Brouw
+
+procedure xt_pw_ncp (wcs, xi_r, eta_r, dra_r, dec_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	xi_r		# i: standard coordinate (radians)
+double	eta_r		# i: standard coordinate (radians)
+double	dra_r		# o: RA of object - RA at reference pixel (radians)
+double	dec_r		# o: declination of object (radians)
+#--
+double	temp
+
+begin
+	temp = W_COSDEC(wcs) - eta_r * W_SINDEC(wcs)
+
+	dra_r = atan2 (-xi_r, temp)
+	dec_r = acos (temp / cos (dra_r))
+	if (W_SINDEC(wcs) < 0)
+	    dec_r = -dec_r
+end
+
+
+# xt_wp_gls -- convert from ra & dec using global-sine projection
+#
+# Reference:  AIPS Memo No. 46 by Eric W. Greisen
+
+procedure xt_wp_gls (wcs, dra_r, dec_r, xi_r, eta_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	dra_r		# i: RA of object - RA at reference pixel (radians)
+double	dec_r		# i: declination of object (radians)
+double	xi_r		# o: standard coordinate (radians)
+double	eta_r		# o: standard coordinate (radians)
+#--
+double	cosdec		# cos of object declination
+double	temp		# delta RA
+int	idec		# which axis is declination axis
+
+begin
+	cosdec = cos (dec_r)
+	idec = W_DEC_AX(wcs)
+
+	temp = dra_r
+
+	# Put dra_r on the interval (-180,+180] degrees.
+	if (temp <= -PI)
+	    temp = temp + TWOPI
+	if (temp > PI)
+	    temp = temp - TWOPI
+
+	xi_r = temp * cosdec
+
+	if (idec > 0)
+	    eta_r = dec_r - DEGTORAD (W_CRVAL(wcs,idec))
+	else
+	    eta_r = dec_r
+end
+
+
+# xt_pw_gls -- convert to ra & dec using global-sine projection
+#
+# Reference:  AIPS Memo No. 46 by Eric W. Greisen
+
+procedure xt_pw_gls (wcs, xi_r, eta_r, dra_r, dec_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	xi_r		# i: standard coordinate (radians)
+double	eta_r		# i: standard coordinate (radians)
+double	dra_r		# o: RA of object - RA at reference pixel (radians)
+double	dec_r		# o: declination of object (radians)
+#--
+double	cosdec		# cosine of object declination
+int	idec		# which axis is declination axis
+
+begin
+	idec = W_DEC_AX(wcs)
+	if (idec > 0)
+	    dec_r = eta_r + DEGTORAD (W_CRVAL(wcs,idec))
+	else
+	    dec_r = eta_r
+
+	cosdec = cos (dec_r)
+	if (cosdec > 0.d0)
+	    dra_r = xi_r / cosdec
+	else
+	    dra_r = 0.d0
+end
+
+# xt_wp_stg -- convert from ra & dec using stereographic projection
+#
+# Reference:  AIPS Memo No. 46 by Eric W. Greisen
+
+procedure xt_wp_stg (wcs, dra_r, dec_r, xi_r, eta_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	dra_r		# i: RA of object - RA at reference pixel (radians)
+double	dec_r		# i: declination of object (radians)
+double	xi_r		# o: standard coordinate (radians)
+double	eta_r		# o: standard coordinate (radians)
+#--
+double	cosdra, sindra	# cos & sin of dra_r
+double	cosdec, sindec	# cos & sin of object declination
+double	cosdist		# cos of dist from ref pixel to object
+double	sincos		# sin (theta) * cos (phi)
+
+begin
+	cosdra = cos (dra_r)
+	sindra = sin (dra_r)
+
+	cosdec = cos (dec_r)
+	sindec = sin (dec_r)
+
+	cosdist = sindec * W_SINDEC(wcs) + cosdec * W_COSDEC(wcs) * cosdra
+	sincos = sindec * W_COSDEC(wcs) - cosdec * W_SINDEC(wcs) * cosdra
+
+	xi_r = 2.d0 * cosdec * sindra / (1.d0 + cosdist)
+	eta_r = 2.d0 * sincos / (1.d0 + cosdist)
+end
+
+
+# xt_pw_stg -- convert to ra & dec using stereographic projection
+
+procedure xt_pw_stg (wcs, xi_r, eta_r, dra_r, dec_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	xi_r		# i: standard coordinate (radians)
+double	eta_r		# i: standard coordinate (radians)
+double	dra_r		# o: RA of object - RA at reference pixel (radians)
+double	dec_r		# o: declination of object (radians)
+#--
+double	rho2		# square of distance from reference pixel
+double	scale		# factor to reduce xi, eta to y, z
+double	x, y, z		# unit vector toward object
+double	temp
+
+begin
+	rho2 = xi_r * xi_r + eta_r * eta_r
+
+	x = (4.d0 - rho2) / (4.d0 + rho2)
+	scale = (x + 1.d0) / 2.d0
+
+	y = xi_r * scale
+	z = eta_r * scale
+
+	temp = x * W_COSDEC(wcs) - z * W_SINDEC(wcs)
+	z    = x * W_SINDEC(wcs) + z * W_COSDEC(wcs)
+	x = temp
+
+	if (x == 0.d0 && y == 0.d0)
+	    dra_r = 0.d0
+	else
+	    dra_r = atan2 (y, x)
+	dec_r = atan2 (z, sqrt (x*x + y*y))
+end
+
+
+# xt_wp_ait -- convert from ra & dec using Aitoff projection
+#
+# Note that the declination at the reference pixel is ignored and is
+# assumed to be zero.  The algorithms given in the AIPS reference do
+# allow for a non-zero declination at the reference pixel.
+#
+# Reference:  AIPS Memo No. 46 by Eric W. Greisen
+
+procedure xt_wp_ait (wcs, dra_r, dec_r, xi_r, eta_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	dra_r		# i: RA of object - RA at reference pixel (radians)
+double	dec_r		# i: declination of object (radians)
+double	xi_r		# o: standard coordinate (radians)
+double	eta_r		# o: standard coordinate (radians)
+#--
+double	z		# temp variable
+double	cosdec		# cosine of declination
+
+begin
+	cosdec = cos (dec_r)
+	z = sqrt ((1.d0 + cosdec * cos (dra_r/2.d0)) / 2.d0)
+
+	xi_r = 2.d0 * cosdec * sin (dra_r/2.d0) / z
+	eta_r = sin (dec_r) / z
+end
+
+
+# xt_pw_ait -- convert to ra & dec using Aitoff projection
+#
+# Note that the declination at the reference pixel is ignored and is
+# assumed to be zero.  The algorithms given in the AIPS reference do
+# allow for a non-zero declination at the reference pixel.
+#
+# Reference:  AIPS Memo No. 46 by Eric W. Greisen
+
+procedure xt_pw_ait (wcs, xi_r, eta_r, dra_r, dec_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	xi_r		# i: standard coordinate (radians)
+double	eta_r		# i: standard coordinate (radians)
+double	dra_r		# o: RA of object - RA at reference pixel (radians)
+double	dec_r		# o: declination of object (radians)
+#--
+double	z		# temp variable
+double	cosdec		# cosine of declination
+
+begin
+	z = sqrt (1.d0 - xi_r*xi_r/16.d0 - eta_r*eta_r/4.d0)
+
+	dec_r = asin (eta_r * z)
+	cosdec = cos (dec_r)
+
+	if (cosdec > 0.d0) {
+	    dra_r = 2.d0 * asin (xi_r * z / (2.d0 * cosdec))
+	} else {
+	    dra_r = 0.d0
+	}
+end
+
+
+# xt_wp_mer -- convert from ra & dec using Mercator projection
+#
+# Note that the declination at the reference pixel is ignored and is
+# assumed to be zero.  The algorithms given in the AIPS reference do
+# allow for a non-zero declination at the reference pixel.
+#
+# Reference:  AIPS Memo No. 46 by Eric W. Greisen
+
+procedure xt_wp_mer (wcs, dra_r, dec_r, xi_r, eta_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	dra_r		# i: RA of object - RA at reference pixel (radians)
+double	dec_r		# i: declination of object (radians)
+double	xi_r		# o: standard coordinate (radians)
+double	eta_r		# o: standard coordinate (radians)
+#--
+double	temp
+
+begin
+	xi_r = dra_r
+	temp = (dec_r + HALFPI) / 2.d0
+	if (temp >= HALFPI || temp <= 0.d0)
+	    call error (1, "invalid declination for Mercator projection")
+	eta_r = log (tan (temp))
+end
+
+
+# xt_pw_mer -- convert to ra & dec using Mercator projection
+#
+# Reference:  AIPS Memo No. 46 by Eric W. Greisen
+
+procedure xt_pw_mer (wcs, xi_r, eta_r, dra_r, dec_r)
+
+pointer wcs		# i: pointer to world coord system struct
+double	xi_r		# i: standard coordinate (radians)
+double	eta_r		# i: standard coordinate (radians)
+double	dra_r		# o: RA of object - RA at reference pixel (radians)
+double	dec_r		# o: declination of object (radians)
+#--
+
+begin
+	dra_r = xi_r
+	dec_r = 2.d0 * atan (exp (eta_r)) - HALFPI
+end