/*** File wcslib/tnxpos.c *** February 14, 2001 *** By Doug Mink, dmink@cfa.harvard.edu *** Harvard-Smithsonian Center for Astrophysics *** After IRAF mwcs/wftnx.x and mwcs/wfgsurfit.x */ #include #include #include #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 (wcs, dir) * tnxdestroy (wcs) * tnxfwd (wcs, v1, v2) Pixels to WCS * tnxrev (wcs, v1, v2) WCS to pixels */ #define max_niter 500 #define SZ_ATSTRING 500 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) 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; 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) x = x + wf_gseval (wcs->lngcor, x, y); if (wcs->latcor != NULL) y = y + wf_gseval (wcs->latcor, x, y); 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); } /* tnxdestroy -- free up the distortion surface pointers */ void tnxdestroy (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. * * 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 != (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->xmaxmin, sf->xrange, sf->xbasis); wf_gsb1pol (y, sf->yorder, sf->ymaxmin, sf->yrange, 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 ncoeff, 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; } ncoeff = 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+1; i--) { for (j = i; j >= i-nyder+1; j--) { for (k = 0; k < sf2->xorder; k++) ptr1[nxder+k] = ptr1[nxder+k] * (double)(j-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); } 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; 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; 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, k1, k2, basis) double x; /*i data point */ int order; /*i order of polynomial, order = 1, constant */ double k1, k2; /*i nomalizing constants, dummy in this case */ 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) */