/*** 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)
 */