rm funtools for update

1 files changed, 0 insertions, 1234 deletions

diff --git a/funtools/wcs/tnxpos.c b/funtools/wcs/tnxpos.c
deleted file mode 100644
index e13d78e..0000000
--- a/funtools/wcs/tnxpos.c
+++ /dev/null
@@ -1,1234 +0,0 @@
-/*** File wcslib/tnxpos.c
- *** September 17, 2008
- *** By Jessica Mink, jmink@cfa.harvard.edu
- *** Harvard-Smithsonian Center for Astrophysics
- *** After IRAF mwcs/wftnx.x and mwcs/wfgsurfit.x
- *** Copyright (C) 1998-2008
- *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA
-
-    This library is free software; you can redistribute it and/or
-    modify it under the terms of the GNU Lesser General Public
-    License as published by the Free Software Foundation; either
-    version 2 of the License, or (at your option) any later version.
-
-    This library is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    Lesser General Public License for more details.
-    
-    You should have received a copy of the GNU Lesser General Public
-    License along with this library; if not, write to the Free Software
-    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
-
-    Correspondence concerning WCSTools should be addressed as follows:
-           Internet email: jmink@cfa.harvard.edu
-           Postal address: Jessica Mink
-                           Smithsonian Astrophysical Observatory
-                           60 Garden St.
-                           Cambridge, MA 02138 USA
- */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include "wcs.h"
-
-#define SPHTOL 0.00001
-#define BADCVAL 0.0
-#define MAX(a,b) (((a) > (b)) ? (a) : (b))
-#define MIN(a,b) (((a) < (b)) ? (a) : (b))
-
-/* wftnx -- wcs function driver for the gnomonic projection with correction.
- *    tnxinit (header, wcs)
- *    tnxclose (wcs)
- *    tnxfwd (xpix, ypix, wcs, xpos, ypos)	Pixels to WCS
- *    tnxrev (xpos, ypos, wcs, xpix, ypix)	WCS to pixels
- */
-
-#define	max_niter	500
-#define	SZ_ATSTRING	2000
-static void wf_gsclose();
-static void wf_gsb1pol();
-static void wf_gsb1leg();
-static void wf_gsb1cheb();
-
-/* tnxinit -- initialize the gnomonic forward or inverse transform.
- * initialization for this transformation consists of, determining which
- * axis is ra / lon and which is dec / lat, computing the celestial longitude
- * and colatitude of the native pole, reading in the the native longitude
- * of the pole of the celestial coordinate system longpole from the attribute
- * list, precomputing euler angles and various intermediaries derived from the
- * coordinate reference values, and reading in the projection parameter ro
- * from the attribute list. if longpole is undefined then a value of 180.0
- * degrees is assumed. if ro is undefined a value of 180.0 / pi is assumed.
- * the tan projection is equivalent to the azp projection with mu set to 0.0.
- * in order to determine the axis order, the parameter "axtype={ra|dec}
- * {xlon|glat}{xlon|elat}" must have been set in the attribute list for the
- * function. the longpole and ro parameters may be set in either or both of
- * the axes attribute lists, but the value in the ra axis attribute list takes
- * precedence. 
- */
-
-int
-tnxinit (header, wcs)
-
-const char *header;	/* FITS header */
-struct WorldCoor *wcs;	/* pointer to WCS structure */
-{
-    struct IRAFsurface *wf_gsopen();
-    char *str1, *str2, *lngstr, *latstr;
-    extern void wcsrotset();
-
-    /* allocate space for the attribute strings */
-    str1 = malloc (SZ_ATSTRING);
-    str2 = malloc (SZ_ATSTRING);
-    hgetm (header, "WAT1", SZ_ATSTRING, str1);
-    hgetm (header, "WAT2", SZ_ATSTRING, str2);
-
-    lngstr = malloc (SZ_ATSTRING);
-    latstr = malloc (SZ_ATSTRING);
-
-    /* determine the native longitude of the pole of the celestial
-	coordinate system corresponding to the FITS keyword longpole.
-	this number has no default and should normally be set to 180
-	degrees. search both axes for this quantity. */
-
-    if (wcs->longpole > 360.0) {
-	if (!igetr8 (str1, "longpole", &wcs->longpole)) {
-	    if (!igetr8 (str2, "longpole", &wcs->longpole))
-		wcs->longpole = 180.0;
-	    }
-	}
-
-    /*  Fetch the ro projection parameter which is the radius of the
-	generating sphere for the projection. if ro is absent which
-	is the usual case set it to 180 / pi. search both axes for
-	this quantity. */
-
-    if (!igetr8 (str1, "ro", &wcs->rodeg)) {
-	if (!igetr8 (str2, "ro", &wcs->rodeg))
-	    wcs->rodeg = 180.0 / PI;
-	}
-
-    /*  Fetch the longitude correction surface. note that the attribute
-	string may be of any length so the length of atvalue may have
-	to be adjusted. */
-
-    if (!igets (str1, "lngcor", SZ_ATSTRING, lngstr)) {
-	if (!igets (str2, "lngcor", SZ_ATSTRING, lngstr))
-	    wcs->lngcor = NULL;
-	else
-	    wcs->lngcor = wf_gsopen (lngstr);
-	}
-    else
-	wcs->lngcor = wf_gsopen (lngstr);
-
-    /*  Fetch the latitude correction surface. note that the attribute
-	string may be of any length so the length of atvalue may have
-	to be adjusted. */
-
-    if (!igets (str2, "latcor", SZ_ATSTRING, latstr)) {
-	if (!igets (str1, "latcor", SZ_ATSTRING, latstr))
-	    wcs->latcor = NULL;
-	else
-	    wcs->latcor = wf_gsopen (latstr);
-	}
-    else
-	wcs->latcor = wf_gsopen (latstr);
-
-    /* Compute image rotation */
-    wcsrotset (wcs);
-
-    /* free working space. */
-    free (str1);
-    free (str2);
-    free (lngstr);
-    free (latstr);
-
-    /* Return 1 if there are no correction coefficients */
-    if (wcs->latcor == NULL && wcs->lngcor == NULL)
-	return (1);
-    else
-	return (0);
-}
-
-
-/* tnxpos -- forward transform (physical to world) gnomonic projection. */
-
-int
-tnxpos (xpix, ypix, wcs, xpos, ypos)
-
-double	xpix, ypix;	/*i physical coordinates (x, y) */
-struct WorldCoor *wcs;	/*i pointer to WCS descriptor */
-double	*xpos, *ypos;	/*o world coordinates (ra, dec) */
-{
-    int	ira, idec;
-    double x, y, r, phi, theta, costhe, sinthe, dphi, cosphi, sinphi, dlng, z;
-    double colatp, coslatp, sinlatp, longp;
-    double xs, ys, ra, dec, xp, yp;
-    double wf_gseval();
-
-    /* Convert from pixels to image coordinates */
-    xpix = xpix - wcs->crpix[0];
-    ypix = ypix - wcs->crpix[1];
-
-    /* Scale and rotate using CD matrix */
-    if (wcs->rotmat) {
-	x = xpix * wcs->cd[0] + ypix * wcs->cd[1];
-	y = xpix * wcs->cd[2] + ypix * wcs->cd[3];
-	}
-
-    else {
-
-	/* Check axis increments - bail out if either 0 */
-	if (wcs->cdelt[0] == 0.0 || wcs->cdelt[1] == 0.0) {
-	    *xpos = 0.0;
-	    *ypos = 0.0;
-	    return 2;
-	    }
-
-	/* Scale using CDELT */
-	xs = xpix * wcs->cdelt[0];
-	ys = ypix * wcs->cdelt[1];
-
-	/* Take out rotation from CROTA */
-	if (wcs->rot != 0.0) {
-	    double cosr = cos (degrad (wcs->rot));
-	    double sinr = sin (degrad (wcs->rot));
-	    x = xs * cosr - ys * sinr;
-	    y = xs * sinr + ys * cosr;
-    	    }
-	else {
-	    x = xs;
-	    y = ys;
-	    }
-	}
-
-    /* get the axis numbers */
-    if (wcs->coorflip) {
-	ira = 1;
-	idec = 0;
-	}
-    else {
-	ira = 0;
-	idec = 1;
-	}
-    colatp = degrad (90.0 - wcs->crval[idec]);
-    coslatp = cos(colatp);
-    sinlatp = sin(colatp);
-    longp = degrad(wcs->longpole);
-
-    /*  Compute native spherical coordinates phi and theta in degrees from the
-	projected coordinates. this is the projection part of the computation */
-    if (wcs->lngcor != NULL)
-	xp = x + wf_gseval (wcs->lngcor, x, y);
-    else
-	xp = x;
-    if (wcs->latcor != NULL)
-	yp = y + wf_gseval (wcs->latcor, x, y);
-    else
-	yp = y;
-    x = xp;
-    y = yp;
-    r = sqrt (x * x + y * y);
-
-    /* Compute phi */
-    if (r == 0.0)
-	phi = 0.0;
-    else
-	phi = atan2 (x, -y);
-
-    /* Compute theta */
-    theta = atan2 (wcs->rodeg, r);
-
-    /*  Compute the celestial coordinates ra and dec from the native
-	coordinates phi and theta. this is the spherical geometry part
-	of the computation */
-
-    costhe = cos (theta);
-    sinthe = sin (theta);
-    dphi = phi - longp;
-    cosphi = cos (dphi);
-    sinphi = sin (dphi);
-
-    /* Compute the ra */
-    x = sinthe * sinlatp - costhe * coslatp * cosphi;
-    if (fabs (x) < SPHTOL)
-	x = -cos (theta + colatp) + costhe * coslatp * (1.0 - cosphi);
-    y = -costhe * sinphi;
-    if (x != 0.0 || y != 0.0)
-	dlng = atan2 (y, x);
-    else
-	dlng = dphi + PI ;
-    ra =  wcs->crval[ira] + raddeg(dlng);
-
-    /* normalize ra */
-    if (wcs->crval[ira] >= 0.0) {
-	if (ra < 0.0)
-	    ra = ra + 360.0;
-	}
-    else {
-	if (ra > 0.0)
-	    ra = ra - 360.0;
-	}
-    if (ra > 360.0)
-	ra = ra - 360.0;
-    else if (ra < -360.0)
-	ra = ra + 360.0;
-
-    /* compute the dec */
-    if (fmod (dphi, PI) == 0.0) {
-	dec = raddeg(theta + cosphi * colatp);
-	if (dec > 90.0)
-	    dec = 180.0 - dec;
-	if (dec < -90.0)
-	    dec = -180.0 - dec;
-	}
-    else {
-	z = sinthe * coslatp + costhe * sinlatp * cosphi;
-	if (fabs(z) > 0.99) {
-	    if (z >= 0.0)
-		dec = raddeg(acos (sqrt(x * x + y * y)));
-	    else
-		dec = raddeg(-acos (sqrt(x * x + y * y)));
-	    }
-	else
-		dec = raddeg(asin (z));
-	}
-
-    /* store the results */
-    *xpos  = ra;
-    *ypos = dec;
-    return (0);
-}
-
-
-/* tnxpix -- inverse transform (world to physical) gnomonic projection */
-
-int
-tnxpix (xpos, ypos, wcs, xpix, ypix)
-
-double	xpos, ypos;	/*i world coordinates (ra, dec) */
-struct WorldCoor *wcs;	/*i pointer to WCS descriptor */
-double	*xpix, *ypix;	/*o physical coordinates (x, y) */
-{
-    int	ira, idec, niter;
-    double ra, dec, cosdec, sindec, cosra, sinra, x, y, phi, theta;
-    double s, r, dphi, z, dpi, dhalfpi, twopi, tx;
-    double xm, ym, f, fx, fy, g, gx, gy, denom, dx, dy;
-    double colatp, coslatp, sinlatp, longp, sphtol;
-    double wf_gseval(), wf_gsder();
-
-    /* get the axis numbers */
-    if (wcs->coorflip) {
-	ira = 1;
-	idec = 0;
-	}
-    else {
-	ira = 0;
-	idec = 1;
-	}
-
-    /*  Compute the transformation from celestial coordinates ra and
-	dec to native coordinates phi and theta. this is the spherical
-	geometry part of the transformation */
-
-    ra  = degrad (xpos - wcs->crval[ira]);
-    dec = degrad (ypos);
-    cosra = cos (ra);
-    sinra = sin (ra);
-    cosdec = cos (dec);
-    sindec = sin (dec);
-    colatp = degrad (90.0 - wcs->crval[idec]);
-    coslatp = cos (colatp);
-    sinlatp = sin (colatp);
-    if (wcs->longpole == 999.0)
-	longp = degrad (180.0);
-    else
-	longp = degrad(wcs->longpole);
-    dpi = PI;
-    dhalfpi = dpi * 0.5;
-    twopi = PI + PI;
-    sphtol = SPHTOL;
-
-    /* Compute phi */
-    x = sindec * sinlatp - cosdec * coslatp * cosra;
-    if (fabs(x) < sphtol)
-	x = -cos (dec + colatp) + cosdec * coslatp * (1.0 - cosra);
-    y = -cosdec * sinra;
-    if (x != 0.0 || y != 0.0)
-	dphi = atan2 (y, x);
-    else
-	dphi = ra - dpi;
-    phi = longp + dphi;
-    if (phi > dpi)
-	phi = phi - twopi;
-    else if (phi < -dpi)
-	phi = phi + twopi;
-
-    /* Compute theta */
-    if (fmod (ra, dpi) == 0.0) {
-	theta = dec + cosra * colatp;
-	if (theta > dhalfpi)
-	    theta = dpi - theta;
-	if (theta < -dhalfpi)
-	    theta = -dpi - theta;
-	}
-    else {
-	z = sindec * coslatp + cosdec * sinlatp * cosra;
-	if (fabs (z) > 0.99) {
-	    if (z >= 0.0)
-		theta = acos (sqrt(x * x + y * y));
-	    else
-		theta = -acos (sqrt(x * x + y * y));
-	    }
-	else
-	    theta = asin (z);
-	}
-
-    /*  Compute the transformation from native coordinates phi and theta
-	to projected coordinates x and y */
-
-    s = sin (theta);
-    if (s == 0.0) {
-	x = BADCVAL;
-	y = BADCVAL;
-	}
-    else {
-	r = wcs->rodeg * cos (theta) / s;
-	if (wcs->lngcor == NULL && wcs->latcor == NULL) {
-	    if (wcs->coorflip) {
-		y  = r * sin (phi);
-		x = -r * cos (phi);
-		}
-	    else {
-		x  = r * sin (phi);
-		y = -r * cos (phi);
-		}
-	    }
-	else {
-	    xm  = r * sin (phi);
-	    ym = -r * cos (phi);
-	    x = xm;
-	    y = ym;
-	    niter = 0;
-	    while (niter < max_niter) {
-		if (wcs->lngcor != NULL) {
-		    f = x + wf_gseval (wcs->lngcor, x, y) - xm;
-		    fx = wf_gsder (wcs->lngcor, x, y, 1, 0);
-		    fx = 1.0 + fx;
-		    fy = wf_gsder (wcs->lngcor, x, y, 0, 1);
-		    }
-		else {
-		    f = x - xm;
-		    fx = 1.0 ;
-		    fy = 0.0;
-		    }
-		if (wcs->latcor != NULL) {
-		    g = y + wf_gseval (wcs->latcor, x, y) - ym;
-		    gx = wf_gsder (wcs->latcor, x, y, 1, 0);
-		    gy = wf_gsder (wcs->latcor, x, y, 0, 1);
-		    gy = 1.0 + gy;
-		    }
-		else {
-		    g = y - ym;
-		    gx = 0.0 ;
-		    gy = 1.0;
-		    }
-
-		denom = fx * gy - fy * gx;
-		if (denom == 0.0)
-		    break;
-		dx = (-f * gy + g * fy) / denom;
-		dy = (-g * fx + f * gx) / denom;
-		x = x + dx;
-		y = y + dy;
-		if (MAX(MAX(fabs(dx),fabs(dy)),MAX(fabs(f),fabs(g))) < 2.80e-8)
-		    break;
-
-		niter = niter + 1;
-		}
-
-	    /* Reverse x and y if axes flipped */
-	    if (wcs->coorflip) {
-		tx = x;
-		x = y;
-		y = tx;
-		}
-	    }
-	}
-
-    /* Scale and rotate using CD matrix */
-    if (wcs->rotmat) {
-	*xpix = x * wcs->dc[0] + y * wcs->dc[1];
-	*ypix = x * wcs->dc[2] + y * wcs->dc[3];
-	}
-
-    else {
-
-	/* Correct for rotation */
-	if (wcs->rot!=0.0) {
-	    double cosr = cos (degrad (wcs->rot));
-	    double sinr = sin (degrad (wcs->rot));
-	    *xpix = x * cosr + y * sinr;
-	    *ypix = y * cosr - x * sinr;
-	    }
-	else {
-	    *xpix = x;
-	    *ypix = y;
-	    }
-
-	/* Scale using CDELT */
-	if (wcs->xinc != 0.)
-	    *xpix = *xpix / wcs->xinc;
-	if (wcs->yinc != 0.)
-	    *ypix = *ypix / wcs->yinc;
-	}
-
-    /* Convert to pixels  */
-    *xpix = *xpix + wcs->xrefpix;
-    *ypix = *ypix + wcs->yrefpix;
-
-    return (0);
-}
-
-
-/* TNXCLOSE -- free up the distortion surface pointers */
-
-void
-tnxclose (wcs)
-
-struct WorldCoor *wcs;		/* pointer to the WCS descriptor */
-
-{
-    if (wcs->lngcor != NULL)
-	wf_gsclose (wcs->lngcor);
-    if (wcs->latcor != NULL)
-	wf_gsclose (wcs->latcor);
-    return;
-}
-
-/* copyright(c) 1986 association of universities for research in astronomy inc.
- * wfgsurfit.x -- surface fitting package used by wcs function drivers.
- * Translated to C from SPP by Jessica Mink, SAO, May 26, 1998
- *
- * the following routines are used by the experimental function drivers tnx
- * and zpx to decode polynomial fits stored in the image header in the form
- * of a list of parameters and coefficients into surface descriptors in
- * ra / dec or longitude latitude. the polynomial surfaces so encoded consist
- * of corrections to function drivers tan and zpn. the package routines are
- * modelled after the equivalent gsurfit routines and are consistent with them.
- * the routines are:
- *
- *                 sf = wf_gsopen (wattstr)
- *                     wf_gsclose (sf)
- *
- *                  z = wf_gseval (sf, x, y)
- *             ncoeff = wf_gscoeff (sf, coeff)
- *               zder = wf_gsder (sf, x, y, nxder, nyder)
- *
- * wf_gsopen is used to open a surface fit encoded in a wcs attribute, returning
- * the sf surface fitting descriptor.  wf_gsclose should be called later to free
- * the descriptor.  wf_gseval is called to evaluate the surface at a point.
- */
-
-
-#define  SZ_GSCOEFFBUF     20
-
-/* define the structure elements for the wf_gsrestore task */
-#define  TNX_SAVETYPE     0
-#define  TNX_SAVEXORDER   1
-#define  TNX_SAVEYORDER   2
-#define  TNX_SAVEXTERMS   3
-#define  TNX_SAVEXMIN     4
-#define  TNX_SAVEXMAX     5
-#define  TNX_SAVEYMIN     6
-#define  TNX_SAVEYMAX     7
-#define  TNX_SAVECOEFF    8
-
-
-/* wf_gsopen -- decode the longitude / latitude or ra / dec mwcs attribute
- * and return a gsurfit compatible surface descriptor.
- */
-
-struct IRAFsurface *
-wf_gsopen (astr)
-
-char    *astr;		/* the input mwcs attribute string */
-
-{
-    double dval;
-    char *estr;
-    int npar, szcoeff;
-    double *coeff;
-    struct IRAFsurface *gs;
-    struct IRAFsurface *wf_gsrestore();
-
-    if (astr[1] == 0)
-	return (NULL);
-
-    gs = NULL;
-    npar = 0;
-    szcoeff = SZ_GSCOEFFBUF;
-    coeff = (double *) malloc (szcoeff * sizeof (double));
-
-    estr = astr;
-    while (*estr != (char) 0) {
-	dval = strtod (astr, &estr);
-	if (*estr == '.')
-	    estr++;
-	if (*estr != (char) 0) {
-	    npar++;
-	    if (npar >= szcoeff) {
-		szcoeff = szcoeff + SZ_GSCOEFFBUF;
-		coeff = (double *) realloc (coeff, (szcoeff * sizeof (double)));
-		}
-	    coeff[npar-1] = dval;
-	    astr = estr;
-	    while (*astr == ' ') astr++;
-	    }
-        }
-
-    gs = wf_gsrestore (coeff);
-
-    free (coeff);
-
-    if (npar == 0)
-	return (NULL);
-    else
-	return (gs);
-}
-
-
-/* wf_gsclose -- procedure to free the surface descriptor */
-
-static void
-wf_gsclose (sf)
-
-struct IRAFsurface *sf;	/* the surface descriptor */
-
-{
-    if (sf != NULL) {
-	if (sf->xbasis != NULL)
-	    free (sf->xbasis);
-	if (sf->ybasis != NULL)
-	    free (sf->ybasis);
-	if (sf->coeff != NULL)
-	    free (sf->coeff);
-	free (sf);
-	}
-    return;
-}
-
-
-/* wf_gseval -- procedure to evaluate the fitted surface at a single point.
- * the wf->ncoeff coefficients are stored in the vector pointed to by sf->coeff.
- */
-
-double
-wf_gseval (sf, x, y)
-
-struct IRAFsurface *sf;	/* pointer to surface descriptor structure */
-double  x;		/* x value */
-double  y;		/* y value */
-{
-    double sum, accum;
-    int i, ii, k, maxorder, xorder;
-
-    /* Calculate the basis functions */
-    switch (sf->type) {
-        case TNX_CHEBYSHEV:
-            wf_gsb1cheb (x, sf->xorder, sf->xmaxmin, sf->xrange, sf->xbasis);
-            wf_gsb1cheb (y, sf->yorder, sf->ymaxmin, sf->yrange, sf->ybasis);
-	    break;
-        case TNX_LEGENDRE:
-            wf_gsb1leg (x, sf->xorder, sf->xmaxmin, sf->xrange, sf->xbasis);
-            wf_gsb1leg (y, sf->yorder, sf->ymaxmin, sf->yrange, sf->ybasis);
-	    break;
-        case TNX_POLYNOMIAL:
-            wf_gsb1pol (x, sf->xorder, sf->xbasis);
-            wf_gsb1pol (y, sf->yorder, sf->ybasis);
-	    break;
-        default:
-            fprintf (stderr,"TNX_GSEVAL: unknown surface type\n");
-	    return (0.0);
-        }
-
-    /* Initialize accumulator basis functions */
-    sum = 0.0;
-
-    /* Loop over y basis functions */
-    if (sf->xorder > sf->yorder)
-	maxorder = sf->xorder + 1;
-    else
-	maxorder = sf->yorder + 1;
-    xorder = sf->xorder;
-    ii = 0;
-
-    for (i = 0; i < sf->yorder; i++) {
-
-	/* Loop over the x basis functions */
-	accum = 0.0;
-	for (k = 0; k < xorder; k++) {
-	    accum = accum + sf->coeff[ii] * sf->xbasis[k];
-	    ii = ii + 1;
-	    }
-	accum = accum * sf->ybasis[i];
-	sum = sum + accum;
-
-        /* Elements of the coefficient vector where neither k = 1 or i = 1
-           are not calculated if sf->xterms = no. */
-        if (sf->xterms == TNX_XNONE)
-            xorder = 1;
-        else if (sf->xterms == TNX_XHALF) {
-            if ((i + 1 + sf->xorder + 1) > maxorder)
-                xorder = xorder - 1;
-	    }
-        }
-
-    return (sum);
-}
-
-
-/* TNX_GSCOEFF -- procedure to fetch the number and magnitude of the coefficients
- * if the sf->xterms = wf_xbi (yes) then the number of coefficients will be
- * (sf->xorder * sf->yorder); if wf_xterms is wf_xtri then the number
- * of coefficients will be (sf->xorder *  sf->yorder - order *
- * (order - 1) / 2) where order is the minimum of the x and yorders;  if
- * sf->xterms = TNX_XNONE then the number of coefficients will be
- * (sf->xorder + sf->yorder - 1).
- */
-
-int
-wf_gscoeff (sf, coeff)
-
-struct IRAFsurface *sf;	/* pointer to the surface fitting descriptor */
-double	*coeff;		/* the coefficients of the fit */
-
-{
-    int ncoeff;		/* the number of coefficients */
-    int i;
-
-    /* Exctract coefficients from data structure and calculate their number */
-    ncoeff = sf->ncoeff;
-    for (i = 0; i < ncoeff; i++)
-	coeff[i] = sf->coeff[i];
-    return (ncoeff);
-}
-
-
-static double *coeff = NULL;
-static int nbcoeff = 0;
-
-/* wf_gsder -- procedure to calculate a new surface which is a derivative of
- * the input surface.
- */
-
-double
-wf_gsder (sf1, x, y, nxd, nyd)
-
-struct IRAFsurface *sf1; /* pointer to the previous surface */
-double	x;		/* x values */
-double	y;		/* y values */
-int	nxd, nyd;	/* order of the derivatives in x and y */
-{
-    int nxder, nyder, i, j, k, nbytes;
-    int order, maxorder1, maxorder2, nmove1, nmove2;
-    struct IRAFsurface *sf2 = 0;
-    double *ptr1, *ptr2;
-    double zfit, norm;
-    double wf_gseval();
-
-    if (sf1 == NULL)
-	return (0.0);
-
-    if (nxd < 0 || nyd < 0) {
-	fprintf (stderr, "TNX_GSDER: order of derivatives cannot be < 0\n");
-	return (0.0);
-	}
-
-    if (nxd == 0 && nyd == 0) {
-	zfit = wf_gseval (sf1, x, y);
-	return (zfit);
-	}
-
-    /* Allocate space for new surface */
-    sf2 = (struct IRAFsurface *) malloc (sizeof (struct IRAFsurface));
-
-    /* Check the order of the derivatives */
-    nxder = MIN (nxd, sf1->xorder - 1);
-    nyder = MIN (nyd, sf1->yorder - 1);
-
-    /* Set up new surface */
-    sf2->type = sf1->type;
-
-    /* Set the derivative surface parameters */
-    if (sf2->type == TNX_LEGENDRE ||
-	sf2->type == TNX_CHEBYSHEV ||
-	sf2->type == TNX_POLYNOMIAL) {
-
-	sf2->xterms = sf1->xterms;
-
-	/* Find the order of the new surface */
-	switch (sf2->xterms) {
-	    case TNX_XNONE: 
-		if (nxder > 0 && nyder > 0) {
-		    sf2->xorder = 1;
-		    sf2->yorder = 1;
-		    sf2->ncoeff = 1;
-		    }
-		else if (nxder > 0) {
-		    sf2->xorder = MAX (1, sf1->xorder - nxder);
-		    sf2->yorder = 1;
-		    sf2->ncoeff = sf2->xorder;
-		    }
-		else if (nyder > 0) {
-		    sf2->xorder = 1;
-		    sf2->yorder = MAX (1, sf1->yorder - nyder);
-		    sf2->ncoeff = sf2->yorder;
-		    }
-		break;
-
-	    case TNX_XHALF:
-		maxorder1 = MAX (sf1->xorder+1, sf1->yorder+1);
-		order = MAX(1, MIN(maxorder1-1-nyder-nxder,sf1->xorder-nxder));
-		sf2->xorder = order;
-		order = MAX(1, MIN(maxorder1-1-nyder-nxder,sf1->yorder-nyder));
-		sf2->yorder = order;
-		order = MIN (sf2->xorder, sf2->yorder);
-		sf2->ncoeff = sf2->xorder * sf2->yorder - (order*(order-1)/2);
-		break;
-
-	    default:
-		sf2->xorder = MAX (1, sf1->xorder - nxder);
-		sf2->yorder = MAX (1, sf1->yorder - nyder);
-		sf2->ncoeff = sf2->xorder * sf2->yorder;
-	    }
-
-	/* define the data limits */
-	sf2->xrange = sf1->xrange;
-	sf2->xmaxmin = sf1->xmaxmin;
-	sf2->yrange = sf1->yrange;
-	sf2->ymaxmin = sf1->ymaxmin;
-	}
-
-    else {
-	fprintf (stderr, "TNX_GSDER: unknown surface type %d\n", sf2->type);
-	return (0.0);
-	}
-
-    /* Allocate space for coefficients and basis functions */
-    nbytes = sf2->ncoeff * sizeof(double);
-    sf2->coeff = (double *) malloc (nbytes);
-    nbytes = sf2->xorder * sizeof(double);
-    sf2->xbasis = (double *) malloc (nbytes);
-    nbytes = sf2->yorder * sizeof(double);
-    sf2->ybasis = (double *) malloc (nbytes);
-
-    /* Get coefficients */
-    nbytes = sf1->ncoeff * sizeof(double);
-    if (nbytes > nbcoeff) {
-	if (nbcoeff > 0)
-	    coeff = (double *) realloc (coeff, nbytes);
-	else
-	    coeff = (double *) malloc (nbytes);
-	nbcoeff = nbytes;
-	}
-    (void) wf_gscoeff (sf1, coeff);
-
-    /* Compute the new coefficients */
-    switch (sf2->xterms) {
-	case TNX_XFULL:
-	    ptr2 = sf2->coeff + (sf2->yorder - 1) * sf2->xorder;
-	    ptr1 = coeff + (sf1->yorder - 1) * sf1->xorder;
-	    for (i = sf1->yorder - 1; i >= nyder; i--) {
-		for (j = i; j >= i-nyder+1; j--) {
-		    for (k = 0; k < sf2->xorder; k++)
-			ptr1[nxder+k] = ptr1[nxder+k] * (double)(j);
-		    }
-		for (j = sf1->xorder; j >= nxder+1; j--) {
-		    for (k = j; k >= j-nxder+1; k--)
-			ptr1[j-1] = ptr1[j-1] * (double)(k - 1);
-		    }
-		for (j = 0; j < sf2->xorder; j++)
-		    ptr2[j] = ptr1[nxder+j];
-		ptr2 = ptr2 - sf2->xorder;
-		ptr1 = ptr1 - sf1->xorder;
-		}
-	    break;
-
-	case TNX_XHALF:
-	    maxorder1 = MAX (sf1->xorder + 1, sf1->yorder + 1);
-	    maxorder2 = MAX (sf2->xorder + 1, sf2->yorder + 1);
-	    ptr2 = sf2->coeff + sf2->ncoeff;
-	    ptr1 = coeff + sf1->ncoeff;
-	    for (i = sf1->yorder; i >= nyder+1; i--) {
-		nmove1 = MAX (0, MIN (maxorder1 - i, sf1->xorder));
-		nmove2 = MAX (0, MIN (maxorder2 - i + nyder, sf2->xorder));
-		ptr1 = ptr1 - nmove1;
-		ptr2 = ptr2 - nmove2;
-		for (j = i; j > i - nyder + 1; j--) {
-		    for (k = 0; k < nmove2; k++)
-			ptr1[nxder+k] = ptr1[nxder+k] * (double)(j-1);
-		    }
-		for (j = nmove1; j >= nxder+1; j--) {
-		    for (k = j;  k >= j-nxder+1; k--)
-			ptr1[j-1] = ptr1[j-1] * (double)(k - 1);
-		    }
-		for (j = 0; j < nmove2; j++)
-		    ptr2[j] = ptr1[nxder+j];
-		}
-	    break;
-
-	default:
-	    if (nxder > 0 && nyder > 0)
-		sf2->coeff[0] = 0.0;
-
-	    else if (nxder > 0) { 
-		ptr1 = coeff;
-		ptr2 = sf2->coeff + sf2->ncoeff - 1;
-		for (j = sf1->xorder; j >= nxder+1; j--) {
-		    for (k = j; k >= j - nxder + 1; k--)
-			ptr1[j-1] = ptr1[j-1] * (double)(k - 1);
-		    ptr2[0] = ptr1[j-1];
-		    ptr2 = ptr2 - 1;
-		    }
-		}
-
-	    else if (nyder > 0) {
-		ptr1 = coeff + sf1->ncoeff - 1;
-		ptr2 = sf2->coeff;
-		for (i = sf1->yorder; i >= nyder + 1; i--) {
-		    for (j = i; j >= i - nyder + 1; j--)
-			*ptr1 = *ptr1 * (double)(j - 1);
-		    ptr1 = ptr1 - 1;
-		    }
-		for (i = 0; i < sf2->ncoeff; i++)
-		    ptr2[i] = ptr1[i+1];
-		}
-	}
-
-    /* evaluate the derivatives */
-    zfit = wf_gseval (sf2, x, y);
-
-    /* normalize */
-    if (sf2->type != TNX_POLYNOMIAL) { 
-	norm = pow (sf2->xrange, (double)nxder) *
-	       pow (sf2->yrange, (double)nyder);
-	zfit = norm * zfit;
-	}
-
-    /* free the space */
-    wf_gsclose (sf2);
-
-    return (zfit);
-}
-
-
-/* wf_gsrestore -- procedure to restore the surface fit encoded in the
-   image header as a list of double precision parameters and coefficients
-   to the surface descriptor for use by the evaluating routines. the
-   surface parameters, surface type, xorder (or number of polynomial
-   terms in x), yorder (or number of polynomial terms in y), xterms,
-   xmin, xmax and ymin and ymax, are stored in the first eight elements
-   of the double array fit, followed by the wf->ncoeff surface coefficients.
- */
-
-struct IRAFsurface *
-wf_gsrestore (fit)
-
-double	*fit;			/* array containing the surface parameters
-				   and coefficients */
-{
-    struct IRAFsurface	*sf;	/* surface descriptor */
-    int	surface_type, xorder, yorder, order, i;
-    double xmin, xmax, ymin, ymax;
-
-    xorder = (int) (fit[TNX_SAVEXORDER] + 0.5);
-    if (xorder < 1) {
-	fprintf (stderr, "wf_gsrestore: illegal x order %d\n", xorder);
-	return (NULL);
-	}
-
-    yorder = (int) (fit[TNX_SAVEYORDER] + 0.5);
-    if (yorder < 1) {
-	fprintf (stderr, "wf_gsrestore: illegal y order %d\n", yorder);
-	return (NULL);
-	}
-
-    xmin = fit[TNX_SAVEXMIN];
-    xmax = fit[TNX_SAVEXMAX];
-    if (xmax <= xmin) {
-	fprintf (stderr, "wf_gsrestore: illegal x range %f-%f\n",xmin,xmax);
-	return (NULL);
-	}
-    ymin = fit[TNX_SAVEYMIN];
-    ymax = fit[TNX_SAVEYMAX];
-    if (ymax <= ymin) {
-	fprintf (stderr, "wf_gsrestore: illegal y range %f-%f\n",ymin,ymax);
-	return (NULL);
-	}
-
-    /* Set surface type dependent surface descriptor parameters */
-    surface_type = (int) (fit[TNX_SAVETYPE] + 0.5);
-
-    if (surface_type == TNX_LEGENDRE ||
-	surface_type == TNX_CHEBYSHEV ||
-	surface_type == TNX_POLYNOMIAL) {
-
-	/* allocate space for the surface descriptor */
-	sf = (struct IRAFsurface *) malloc (sizeof (struct IRAFsurface));
-	sf->xorder = xorder;
-	sf->xrange = 2.0 / (xmax - xmin);
-	sf->xmaxmin =  - (xmax + xmin) / 2.0;
-	sf->yorder = yorder;
-	sf->yrange = 2.0 / (ymax - ymin);
-	sf->ymaxmin =  - (ymax + ymin) / 2.0;
-	sf->xterms = fit[TNX_SAVEXTERMS];
-	switch (sf->xterms) {
-	    case TNX_XNONE:
-		sf->ncoeff = sf->xorder + sf->yorder - 1;
-		break;
-	    case TNX_XHALF:
-		order = MIN (xorder, yorder);
-		sf->ncoeff = sf->xorder * sf->yorder - order * (order-1) / 2;
-		break;
-	    case TNX_XFULL:
-		sf->ncoeff = sf->xorder * sf->yorder;
-		break;
-	    }
-	}
-    else {
-	fprintf (stderr, "wf_gsrestore: unknown surface type %d\n", surface_type);
-	return (NULL);
-	}
-
-    /* Set remaining curve parameters */
-    sf->type = surface_type;
-
-    /* Restore coefficient array */
-    sf->coeff = (double *) malloc (sf->ncoeff*sizeof (double));
-    for (i = 0; i < sf->ncoeff; i++)
-	sf->coeff[i] = fit[TNX_SAVECOEFF+i];
-
-    /* Allocate space for basis vectors */
-    sf->xbasis = (double *) malloc (sf->xorder*sizeof (double));
-    sf->ybasis = (double *) malloc (sf->yorder*sizeof (double));
-
-    return (sf);
-}
-
-
-/* wf_gsb1pol -- procedure to evaluate all the non-zero polynomial functions
-   for a single point and given order. */
-
-static void
-wf_gsb1pol (x, order, basis)
-
-double  x;		/*i data point */
-int     order;		/*i order of polynomial, order = 1, constant */
-double  *basis;		/*o basis functions */
-{
-    int     i;
-
-    basis[0] = 1.0;
-    if (order == 1)
-	return;
-
-    basis[1] = x;
-    if (order == 2)
-	return;
-
-    for (i = 2; i < order; i++)
-	basis[i] = x * basis[i-1];
-
-    return;
-}
-
-
-/* wf_gsb1leg -- procedure to evaluate all the non-zero legendre functions for
-   a single point and given order. */
-
-static void
-wf_gsb1leg (x, order, k1, k2, basis)
-
-double  x;		/*i data point */
-int     order;		/*i order of polynomial, order = 1, constant */
-double  k1, k2;		/*i normalizing constants */
-double	*basis;		/*o basis functions */
-{
-    int i;
-    double ri, xnorm;
-
-    basis[0] = 1.0;
-    if (order == 1)
-	return;
-
-    xnorm = (x + k1) * k2 ;
-    basis[1] = xnorm;
-    if (order == 2)
-        return;
-
-    for (i = 2; i < order; i++) {
-	ri = i;
-        basis[i] = ((2.0 * ri - 1.0) * xnorm * basis[i-1] -
-                       (ri - 1.0) * basis[i-2]) / ri;
-        }
-
-    return;
-}
-
-
-/* wf_gsb1cheb -- procedure to evaluate all the non-zero chebyshev function
-   coefficients for a given x and order. */
-
-static void
-wf_gsb1cheb (x, order, k1, k2, basis)
-
-double	x;		/*i number of data points */
-int	order;		/*i order of polynomial, 1 is a constant */
-double	k1, k2;		/*i normalizing constants */
-double	*basis;		/*o array of basis functions */
-{
-    int i;
-    double xnorm;
-
-    basis[0] = 1.0;
-    if (order == 1)
-	return;
-
-    xnorm = (x + k1) * k2;
-    basis[1] = xnorm;
-    if (order == 2)
-	return;
-
-    for (i = 2; i < order; i++)
-	basis[i] = 2. * xnorm * basis[i-1] - basis[i-2];
-
-    return;
-}
-
-/* Set surface polynomial from arguments */
-
-int
-tnxpset (wcs, xorder, yorder, xterms, coeff)
-
-struct WorldCoor *wcs;  /* World coordinate system structure */
-int	xorder;		/* Number of x coefficients (same for x and y) */
-int	yorder;		/* Number of y coefficients (same for x and y) */
-int	xterms;		/* Number of xy coefficients (same for x and y) */
-double	*coeff;		/* Plate fit coefficients */
-
-{
-    double *ycoeff;
-    struct IRAFsurface *wf_gspset ();
-
-    wcs->prjcode = WCS_TNX;
-
-    wcs->lngcor = wf_gspset (xorder, yorder, xterms, coeff);
-    ycoeff = coeff + wcs->lngcor->ncoeff;
-    wcs->latcor = wf_gspset (xorder, yorder, xterms, ycoeff);
-
-    return 0;
-}
-
-
-/* wf_gspset -- procedure to set the surface descriptor for use by the
-   evaluating routines.  from arguments.  The surface parameters are
-   surface type, xorder (number of polynomial terms in x), yorder (number
-   of polynomial terms in y), xterms, and the surface coefficients.
- */
-
-struct IRAFsurface *
-wf_gspset (xorder, yorder, xterms, coeff)
-
-int	xorder;
-int	yorder;
-int	xterms;
-double	*coeff;
-{
-    struct IRAFsurface	*sf;	/* surface descriptor */
-    int	surface_type, order, i;
-    double xmin, xmax;
-    double ymin, ymax;
-
-    surface_type = TNX_POLYNOMIAL;
-    xmin = 0.0;
-    xmax = 0.0;
-    ymin = 0.0;
-    ymax = 0.0;
-
-    if (surface_type == TNX_LEGENDRE ||
-	surface_type == TNX_CHEBYSHEV ||
-	surface_type == TNX_POLYNOMIAL) {
-
-	/* allocate space for the surface descriptor */
-	sf = (struct IRAFsurface *) malloc (sizeof (struct IRAFsurface));
-	sf->xorder = xorder;
-	sf->xrange = 2.0 / (xmax - xmin);
-	sf->xmaxmin =  -(xmax + xmin) / 2.0;
-	sf->yorder = yorder;
-	sf->yrange = 2.0 / (ymax - ymin);
-	sf->ymaxmin =  - (ymax + ymin) / 2.0;
-	sf->xterms = xterms;
-	switch (sf->xterms) {
-	    case TNX_XNONE:
-		sf->ncoeff = sf->xorder + sf->yorder - 1;
-		break;
-	    case TNX_XHALF:
-		order = MIN (xorder, yorder);
-		sf->ncoeff = sf->xorder * sf->yorder - order * (order-1) / 2;
-		break;
-	    case TNX_XFULL:
-		sf->ncoeff = sf->xorder * sf->yorder;
-		break;
-	    }
-	}
-    else {
-	fprintf (stderr, "TNX_GSSET: unknown surface type %d\n", surface_type);
-	return (NULL);
-	}
-
-    /* Set remaining curve parameters */
-    sf->type = surface_type;
-
-    /* Restore coefficient array */
-    sf->coeff = (double *) malloc (sf->ncoeff*sizeof (double));
-    for (i = 0; i < sf->ncoeff; i++)
-	sf->coeff[i] = coeff[i];
-
-    /* Allocate space for basis vectors */
-    sf->xbasis = (double *) malloc (sf->xorder*sizeof (double));
-    sf->ybasis = (double *) malloc (sf->yorder*sizeof (double));
-
-    return (sf);
-}
-
-/* Mar 26 1998	New subroutines, translated from SPP
- * Apr 28 1998  Change all local flags to TNX_* and projection flag to WCS_TNX
- * May 11 1998	Fix use of pole longitude default
- * Sep  4 1998	Fix missed assignment in tnxpos from Allen Harris, SAO
- * Sep 10 1998	Fix bugs in tnxpix()
- * Sep 10 1998	Fix missed assignment in tnxpix from Allen Harris, SAO
- *
- * Oct 22 1999	Drop unused variables, fix case statements after lint
- * Dec 10 1999	Fix bug in gsder() which failed to allocate enough memory
- * Dec 10 1999	Compute wcs->rot using wcsrotset() in tnxinit()
- *
- * Feb 14 2001	Fixed off-by-one bug in legendre evaluation (Mike Jarvis)
- *
- * Apr 11 2002	Fix bug when .-terminated substring in wf_gsopen()
- * Apr 29 2002	Clean up code
- * Jun 26 2002	Increase size of WAT strings from 500 to 2000
- *
- * Jun 27 2005	Drop unused arguments k1 and k2 from wf_gsb1pol()
- *
- * Jan  8 2007	Drop unused variable ncoeff in wf_gsder()
- * Jan  9 2007	Declare header const char in tnxinit()
- * Apr  3 2007	Fix offsets to hit last cooefficient in wf_gsder()
- *
- * Sep  5 2008	Fix wf_gseval() call in tnxpos() so unmodified x and y are used
- * Sep  9 2008	Fix loop in TNX_XFULL section of wf_gsder()
- * 		(last two bugs found by Ed Los)
- * Sep 17 2008	Fix tnxpos for null correction case (fix by Ed Los)
- */