1 /*** File wcslib/tnxpos.c
2 *** September 17, 2008
3 *** By Jessica Mink, jmink@cfa.harvard.edu
4 *** Harvard-Smithsonian Center for Astrophysics
5 *** After IRAF mwcs/wftnx.x and mwcs/wfgsurfit.x
6 *** Copyright (C) 1998-2008
7 *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA
8
9 This library is free software; you can redistribute it and/or
10 modify it under the terms of the GNU Lesser General Public
11 License as published by the Free Software Foundation; either
12 version 2 of the License, or (at your option) any later version.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public
20 License along with this library; if not, write to the Free Software
21 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22
23 Correspondence concerning WCSTools should be addressed as follows:
24 Internet email: jmink@cfa.harvard.edu
25 Postal address: Jessica Mink
26 Smithsonian Astrophysical Observatory
27 60 Garden St.
28 Cambridge, MA 02138 USA
29 */
30
31 #include <stdio.h>
32 #include <stdlib.h>
33 #include <math.h>
34 #include "wcs.h"
35
36 #define SPHTOL 0.00001
37 #define BADCVAL 0.0
38 #define MAX(a,b) (((a) > (b)) ? (a) : (b))
39 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
40
41 /* wftnx -- wcs function driver for the gnomonic projection with correction.
42 * tnxinit (header, wcs)
43 * tnxclose (wcs)
44 * tnxfwd (xpix, ypix, wcs, xpos, ypos) Pixels to WCS
45 * tnxrev (xpos, ypos, wcs, xpix, ypix) WCS to pixels
46 */
47
48 #define max_niter 500
49 #define SZ_ATSTRING 2000
50 static void wf_gsclose();
51 static void wf_gsb1pol();
52 static void wf_gsb1leg();
53 static void wf_gsb1cheb();
54
55 /* tnxinit -- initialize the gnomonic forward or inverse transform.
56 * initialization for this transformation consists of, determining which
57 * axis is ra / lon and which is dec / lat, computing the celestial longitude
58 * and colatitude of the native pole, reading in the the native longitude
59 * of the pole of the celestial coordinate system longpole from the attribute
60 * list, precomputing euler angles and various intermediaries derived from the
61 * coordinate reference values, and reading in the projection parameter ro
62 * from the attribute list. if longpole is undefined then a value of 180.0
63 * degrees is assumed. if ro is undefined a value of 180.0 / pi is assumed.
64 * the tan projection is equivalent to the azp projection with mu set to 0.0.
65 * in order to determine the axis order, the parameter "axtype={ra|dec}
66 * {xlon|glat}{xlon|elat}" must have been set in the attribute list for the
67 * function. the longpole and ro parameters may be set in either or both of
68 * the axes attribute lists, but the value in the ra axis attribute list takes
69 * precedence.
70 */
71
72 int
tnxinit(header,wcs)73 tnxinit (header, wcs)
74
75 const char *header; /* FITS header */
76 struct WorldCoor *wcs; /* pointer to WCS structure */
77 {
78 struct IRAFsurface *wf_gsopen();
79 char *str1, *str2, *lngstr, *latstr;
80 extern void wcsrotset();
81
82 /* allocate space for the attribute strings */
83 str1 = malloc (SZ_ATSTRING);
84 str2 = malloc (SZ_ATSTRING);
85 hgetm (header, "WAT1", SZ_ATSTRING, str1);
86 hgetm (header, "WAT2", SZ_ATSTRING, str2);
87
88 lngstr = malloc (SZ_ATSTRING);
89 latstr = malloc (SZ_ATSTRING);
90
91 /* determine the native longitude of the pole of the celestial
92 coordinate system corresponding to the FITS keyword longpole.
93 this number has no default and should normally be set to 180
94 degrees. search both axes for this quantity. */
95
96 if (wcs->longpole > 360.0) {
97 if (!igetr8 (str1, "longpole", &wcs->longpole)) {
98 if (!igetr8 (str2, "longpole", &wcs->longpole))
99 wcs->longpole = 180.0;
100 }
101 }
102
103 /* Fetch the ro projection parameter which is the radius of the
104 generating sphere for the projection. if ro is absent which
105 is the usual case set it to 180 / pi. search both axes for
106 this quantity. */
107
108 if (!igetr8 (str1, "ro", &wcs->rodeg)) {
109 if (!igetr8 (str2, "ro", &wcs->rodeg))
110 wcs->rodeg = 180.0 / PI;
111 }
112
113 /* Fetch the longitude correction surface. note that the attribute
114 string may be of any length so the length of atvalue may have
115 to be adjusted. */
116
117 if (!igets (str1, "lngcor", SZ_ATSTRING, lngstr)) {
118 if (!igets (str2, "lngcor", SZ_ATSTRING, lngstr))
119 wcs->lngcor = NULL;
120 else
121 wcs->lngcor = wf_gsopen (lngstr);
122 }
123 else
124 wcs->lngcor = wf_gsopen (lngstr);
125
126 /* Fetch the latitude correction surface. note that the attribute
127 string may be of any length so the length of atvalue may have
128 to be adjusted. */
129
130 if (!igets (str2, "latcor", SZ_ATSTRING, latstr)) {
131 if (!igets (str1, "latcor", SZ_ATSTRING, latstr))
132 wcs->latcor = NULL;
133 else
134 wcs->latcor = wf_gsopen (latstr);
135 }
136 else
137 wcs->latcor = wf_gsopen (latstr);
138
139 /* Compute image rotation */
140 wcsrotset (wcs);
141
142 /* free working space. */
143 free (str1);
144 free (str2);
145 free (lngstr);
146 free (latstr);
147
148 /* Return 1 if there are no correction coefficients */
149 if (wcs->latcor == NULL && wcs->lngcor == NULL)
150 return (1);
151 else
152 return (0);
153 }
154
155
156 /* tnxpos -- forward transform (physical to world) gnomonic projection. */
157
158 int
tnxpos(xpix,ypix,wcs,xpos,ypos)159 tnxpos (xpix, ypix, wcs, xpos, ypos)
160
161 double xpix, ypix; /*i physical coordinates (x, y) */
162 struct WorldCoor *wcs; /*i pointer to WCS descriptor */
163 double *xpos, *ypos; /*o world coordinates (ra, dec) */
164 {
165 int ira, idec;
166 double x, y, r, phi, theta, costhe, sinthe, dphi, cosphi, sinphi, dlng, z;
167 double colatp, coslatp, sinlatp, longp;
168 double xs, ys, ra, dec, xp, yp;
169 double wf_gseval();
170
171 /* Convert from pixels to image coordinates */
172 xpix = xpix - wcs->crpix[0];
173 ypix = ypix - wcs->crpix[1];
174
175 /* Scale and rotate using CD matrix */
176 if (wcs->rotmat) {
177 x = xpix * wcs->cd[0] + ypix * wcs->cd[1];
178 y = xpix * wcs->cd[2] + ypix * wcs->cd[3];
179 }
180
181 else {
182
183 /* Check axis increments - bail out if either 0 */
184 if (wcs->cdelt[0] == 0.0 || wcs->cdelt[1] == 0.0) {
185 *xpos = 0.0;
186 *ypos = 0.0;
187 return 2;
188 }
189
190 /* Scale using CDELT */
191 xs = xpix * wcs->cdelt[0];
192 ys = ypix * wcs->cdelt[1];
193
194 /* Take out rotation from CROTA */
195 if (wcs->rot != 0.0) {
196 double cosr = cos (degrad (wcs->rot));
197 double sinr = sin (degrad (wcs->rot));
198 x = xs * cosr - ys * sinr;
199 y = xs * sinr + ys * cosr;
200 }
201 else {
202 x = xs;
203 y = ys;
204 }
205 }
206
207 /* get the axis numbers */
208 if (wcs->coorflip) {
209 ira = 1;
210 idec = 0;
211 }
212 else {
213 ira = 0;
214 idec = 1;
215 }
216 colatp = degrad (90.0 - wcs->crval[idec]);
217 coslatp = cos(colatp);
218 sinlatp = sin(colatp);
219 longp = degrad(wcs->longpole);
220
221 /* Compute native spherical coordinates phi and theta in degrees from the
222 projected coordinates. this is the projection part of the computation */
223 if (wcs->lngcor != NULL)
224 xp = x + wf_gseval (wcs->lngcor, x, y);
225 else
226 xp = x;
227 if (wcs->latcor != NULL)
228 yp = y + wf_gseval (wcs->latcor, x, y);
229 else
230 yp = y;
231 x = xp;
232 y = yp;
233 r = sqrt (x * x + y * y);
234
235 /* Compute phi */
236 if (r == 0.0)
237 phi = 0.0;
238 else
239 phi = atan2 (x, -y);
240
241 /* Compute theta */
242 theta = atan2 (wcs->rodeg, r);
243
244 /* Compute the celestial coordinates ra and dec from the native
245 coordinates phi and theta. this is the spherical geometry part
246 of the computation */
247
248 costhe = cos (theta);
249 sinthe = sin (theta);
250 dphi = phi - longp;
251 cosphi = cos (dphi);
252 sinphi = sin (dphi);
253
254 /* Compute the ra */
255 x = sinthe * sinlatp - costhe * coslatp * cosphi;
256 if (fabs (x) < SPHTOL)
257 x = -cos (theta + colatp) + costhe * coslatp * (1.0 - cosphi);
258 y = -costhe * sinphi;
259 if (x != 0.0 || y != 0.0)
260 dlng = atan2 (y, x);
261 else
262 dlng = dphi + PI ;
263 ra = wcs->crval[ira] + raddeg(dlng);
264
265 /* normalize ra */
266 if (wcs->crval[ira] >= 0.0) {
267 if (ra < 0.0)
268 ra = ra + 360.0;
269 }
270 else {
271 if (ra > 0.0)
272 ra = ra - 360.0;
273 }
274 if (ra > 360.0)
275 ra = ra - 360.0;
276 else if (ra < -360.0)
277 ra = ra + 360.0;
278
279 /* compute the dec */
280 if (fmod (dphi, PI) == 0.0) {
281 dec = raddeg(theta + cosphi * colatp);
282 if (dec > 90.0)
283 dec = 180.0 - dec;
284 if (dec < -90.0)
285 dec = -180.0 - dec;
286 }
287 else {
288 z = sinthe * coslatp + costhe * sinlatp * cosphi;
289 if (fabs(z) > 0.99) {
290 if (z >= 0.0)
291 dec = raddeg(acos (sqrt(x * x + y * y)));
292 else
293 dec = raddeg(-acos (sqrt(x * x + y * y)));
294 }
295 else
296 dec = raddeg(asin (z));
297 }
298
299 /* store the results */
300 *xpos = ra;
301 *ypos = dec;
302 return (0);
303 }
304
305
306 /* tnxpix -- inverse transform (world to physical) gnomonic projection */
307
308 int
tnxpix(xpos,ypos,wcs,xpix,ypix)309 tnxpix (xpos, ypos, wcs, xpix, ypix)
310
311 double xpos, ypos; /*i world coordinates (ra, dec) */
312 struct WorldCoor *wcs; /*i pointer to WCS descriptor */
313 double *xpix, *ypix; /*o physical coordinates (x, y) */
314 {
315 int ira, idec, niter;
316 double ra, dec, cosdec, sindec, cosra, sinra, x, y, phi, theta;
317 double s, r, dphi, z, dpi, dhalfpi, twopi, tx;
318 double xm, ym, f, fx, fy, g, gx, gy, denom, dx, dy;
319 double colatp, coslatp, sinlatp, longp, sphtol;
320 double wf_gseval(), wf_gsder();
321
322 /* get the axis numbers */
323 if (wcs->coorflip) {
324 ira = 1;
325 idec = 0;
326 }
327 else {
328 ira = 0;
329 idec = 1;
330 }
331
332 /* Compute the transformation from celestial coordinates ra and
333 dec to native coordinates phi and theta. this is the spherical
334 geometry part of the transformation */
335
336 ra = degrad (xpos - wcs->crval[ira]);
337 dec = degrad (ypos);
338 cosra = cos (ra);
339 sinra = sin (ra);
340 cosdec = cos (dec);
341 sindec = sin (dec);
342 colatp = degrad (90.0 - wcs->crval[idec]);
343 coslatp = cos (colatp);
344 sinlatp = sin (colatp);
345 if (wcs->longpole == 999.0)
346 longp = degrad (180.0);
347 else
348 longp = degrad(wcs->longpole);
349 dpi = PI;
350 dhalfpi = dpi * 0.5;
351 twopi = PI + PI;
352 sphtol = SPHTOL;
353
354 /* Compute phi */
355 x = sindec * sinlatp - cosdec * coslatp * cosra;
356 if (fabs(x) < sphtol)
357 x = -cos (dec + colatp) + cosdec * coslatp * (1.0 - cosra);
358 y = -cosdec * sinra;
359 if (x != 0.0 || y != 0.0)
360 dphi = atan2 (y, x);
361 else
362 dphi = ra - dpi;
363 phi = longp + dphi;
364 if (phi > dpi)
365 phi = phi - twopi;
366 else if (phi < -dpi)
367 phi = phi + twopi;
368
369 /* Compute theta */
370 if (fmod (ra, dpi) == 0.0) {
371 theta = dec + cosra * colatp;
372 if (theta > dhalfpi)
373 theta = dpi - theta;
374 if (theta < -dhalfpi)
375 theta = -dpi - theta;
376 }
377 else {
378 z = sindec * coslatp + cosdec * sinlatp * cosra;
379 if (fabs (z) > 0.99) {
380 if (z >= 0.0)
381 theta = acos (sqrt(x * x + y * y));
382 else
383 theta = -acos (sqrt(x * x + y * y));
384 }
385 else
386 theta = asin (z);
387 }
388
389 /* Compute the transformation from native coordinates phi and theta
390 to projected coordinates x and y */
391
392 s = sin (theta);
393 if (s == 0.0) {
394 x = BADCVAL;
395 y = BADCVAL;
396 }
397 else {
398 r = wcs->rodeg * cos (theta) / s;
399 if (wcs->lngcor == NULL && wcs->latcor == NULL) {
400 if (wcs->coorflip) {
401 y = r * sin (phi);
402 x = -r * cos (phi);
403 }
404 else {
405 x = r * sin (phi);
406 y = -r * cos (phi);
407 }
408 }
409 else {
410 xm = r * sin (phi);
411 ym = -r * cos (phi);
412 x = xm;
413 y = ym;
414 niter = 0;
415 while (niter < max_niter) {
416 if (wcs->lngcor != NULL) {
417 f = x + wf_gseval (wcs->lngcor, x, y) - xm;
418 fx = wf_gsder (wcs->lngcor, x, y, 1, 0);
419 fx = 1.0 + fx;
420 fy = wf_gsder (wcs->lngcor, x, y, 0, 1);
421 }
422 else {
423 f = x - xm;
424 fx = 1.0 ;
425 fy = 0.0;
426 }
427 if (wcs->latcor != NULL) {
428 g = y + wf_gseval (wcs->latcor, x, y) - ym;
429 gx = wf_gsder (wcs->latcor, x, y, 1, 0);
430 gy = wf_gsder (wcs->latcor, x, y, 0, 1);
431 gy = 1.0 + gy;
432 }
433 else {
434 g = y - ym;
435 gx = 0.0 ;
436 gy = 1.0;
437 }
438
439 denom = fx * gy - fy * gx;
440 if (denom == 0.0)
441 break;
442 dx = (-f * gy + g * fy) / denom;
443 dy = (-g * fx + f * gx) / denom;
444 x = x + dx;
445 y = y + dy;
446 if (MAX(MAX(fabs(dx),fabs(dy)),MAX(fabs(f),fabs(g))) < 2.80e-8)
447 break;
448
449 niter = niter + 1;
450 }
451
452 /* Reverse x and y if axes flipped */
453 if (wcs->coorflip) {
454 tx = x;
455 x = y;
456 y = tx;
457 }
458 }
459 }
460
461 /* Scale and rotate using CD matrix */
462 if (wcs->rotmat) {
463 *xpix = x * wcs->dc[0] + y * wcs->dc[1];
464 *ypix = x * wcs->dc[2] + y * wcs->dc[3];
465 }
466
467 else {
468
469 /* Correct for rotation */
470 if (wcs->rot!=0.0) {
471 double cosr = cos (degrad (wcs->rot));
472 double sinr = sin (degrad (wcs->rot));
473 *xpix = x * cosr + y * sinr;
474 *ypix = y * cosr - x * sinr;
475 }
476 else {
477 *xpix = x;
478 *ypix = y;
479 }
480
481 /* Scale using CDELT */
482 if (wcs->xinc != 0.)
483 *xpix = *xpix / wcs->xinc;
484 if (wcs->yinc != 0.)
485 *ypix = *ypix / wcs->yinc;
486 }
487
488 /* Convert to pixels */
489 *xpix = *xpix + wcs->xrefpix;
490 *ypix = *ypix + wcs->yrefpix;
491
492 return (0);
493 }
494
495
496 /* TNXCLOSE -- free up the distortion surface pointers */
497
498 void
tnxclose(wcs)499 tnxclose (wcs)
500
501 struct WorldCoor *wcs; /* pointer to the WCS descriptor */
502
503 {
504 if (wcs->lngcor != NULL)
505 wf_gsclose (wcs->lngcor);
506 if (wcs->latcor != NULL)
507 wf_gsclose (wcs->latcor);
508 return;
509 }
510
511 /* copyright(c) 1986 association of universities for research in astronomy inc.
512 * wfgsurfit.x -- surface fitting package used by wcs function drivers.
513 * Translated to C from SPP by Jessica Mink, SAO, May 26, 1998
514 *
515 * the following routines are used by the experimental function drivers tnx
516 * and zpx to decode polynomial fits stored in the image header in the form
517 * of a list of parameters and coefficients into surface descriptors in
518 * ra / dec or longitude latitude. the polynomial surfaces so encoded consist
519 * of corrections to function drivers tan and zpn. the package routines are
520 * modelled after the equivalent gsurfit routines and are consistent with them.
521 * the routines are:
522 *
523 * sf = wf_gsopen (wattstr)
524 * wf_gsclose (sf)
525 *
526 * z = wf_gseval (sf, x, y)
527 * ncoeff = wf_gscoeff (sf, coeff)
528 * zder = wf_gsder (sf, x, y, nxder, nyder)
529 *
530 * wf_gsopen is used to open a surface fit encoded in a wcs attribute, returning
531 * the sf surface fitting descriptor. wf_gsclose should be called later to free
532 * the descriptor. wf_gseval is called to evaluate the surface at a point.
533 */
534
535
536 #define SZ_GSCOEFFBUF 20
537
538 /* define the structure elements for the wf_gsrestore task */
539 #define TNX_SAVETYPE 0
540 #define TNX_SAVEXORDER 1
541 #define TNX_SAVEYORDER 2
542 #define TNX_SAVEXTERMS 3
543 #define TNX_SAVEXMIN 4
544 #define TNX_SAVEXMAX 5
545 #define TNX_SAVEYMIN 6
546 #define TNX_SAVEYMAX 7
547 #define TNX_SAVECOEFF 8
548
549
550 /* wf_gsopen -- decode the longitude / latitude or ra / dec mwcs attribute
551 * and return a gsurfit compatible surface descriptor.
552 */
553
554 struct IRAFsurface *
wf_gsopen(astr)555 wf_gsopen (astr)
556
557 char *astr; /* the input mwcs attribute string */
558
559 {
560 double dval;
561 char *estr;
562 int npar, szcoeff;
563 double *coeff;
564 struct IRAFsurface *gs;
565 struct IRAFsurface *wf_gsrestore();
566
567 if (astr[1] == 0)
568 return (NULL);
569
570 gs = NULL;
571 npar = 0;
572 szcoeff = SZ_GSCOEFFBUF;
573 coeff = (double *) malloc (szcoeff * sizeof (double));
574
575 estr = astr;
576 while (*estr != (char) 0) {
577 dval = strtod (astr, &estr);
578 if (*estr == '.')
579 estr++;
580 if (*estr != (char) 0) {
581 npar++;
582 if (npar >= szcoeff) {
583 szcoeff = szcoeff + SZ_GSCOEFFBUF;
584 coeff = (double *) realloc (coeff, (szcoeff * sizeof (double)));
585 }
586 coeff[npar-1] = dval;
587 astr = estr;
588 while (*astr == ' ') astr++;
589 }
590 }
591
592 gs = wf_gsrestore (coeff);
593
594 free (coeff);
595
596 if (npar == 0)
597 return (NULL);
598 else
599 return (gs);
600 }
601
602
603 /* wf_gsclose -- procedure to free the surface descriptor */
604
605 static void
wf_gsclose(sf)606 wf_gsclose (sf)
607
608 struct IRAFsurface *sf; /* the surface descriptor */
609
610 {
611 if (sf != NULL) {
612 if (sf->xbasis != NULL)
613 free (sf->xbasis);
614 if (sf->ybasis != NULL)
615 free (sf->ybasis);
616 if (sf->coeff != NULL)
617 free (sf->coeff);
618 free (sf);
619 }
620 return;
621 }
622
623
624 /* wf_gseval -- procedure to evaluate the fitted surface at a single point.
625 * the wf->ncoeff coefficients are stored in the vector pointed to by sf->coeff.
626 */
627
628 double
wf_gseval(sf,x,y)629 wf_gseval (sf, x, y)
630
631 struct IRAFsurface *sf; /* pointer to surface descriptor structure */
632 double x; /* x value */
633 double y; /* y value */
634 {
635 double sum, accum;
636 int i, ii, k, maxorder, xorder;
637
638 /* Calculate the basis functions */
639 switch (sf->type) {
640 case TNX_CHEBYSHEV:
641 wf_gsb1cheb (x, sf->xorder, sf->xmaxmin, sf->xrange, sf->xbasis);
642 wf_gsb1cheb (y, sf->yorder, sf->ymaxmin, sf->yrange, sf->ybasis);
643 break;
644 case TNX_LEGENDRE:
645 wf_gsb1leg (x, sf->xorder, sf->xmaxmin, sf->xrange, sf->xbasis);
646 wf_gsb1leg (y, sf->yorder, sf->ymaxmin, sf->yrange, sf->ybasis);
647 break;
648 case TNX_POLYNOMIAL:
649 wf_gsb1pol (x, sf->xorder, sf->xbasis);
650 wf_gsb1pol (y, sf->yorder, sf->ybasis);
651 break;
652 default:
653 fprintf (stderr,"TNX_GSEVAL: unknown surface type\n");
654 return (0.0);
655 }
656
657 /* Initialize accumulator basis functions */
658 sum = 0.0;
659
660 /* Loop over y basis functions */
661 if (sf->xorder > sf->yorder)
662 maxorder = sf->xorder + 1;
663 else
664 maxorder = sf->yorder + 1;
665 xorder = sf->xorder;
666 ii = 0;
667
668 for (i = 0; i < sf->yorder; i++) {
669
670 /* Loop over the x basis functions */
671 accum = 0.0;
672 for (k = 0; k < xorder; k++) {
673 accum = accum + sf->coeff[ii] * sf->xbasis[k];
674 ii = ii + 1;
675 }
676 accum = accum * sf->ybasis[i];
677 sum = sum + accum;
678
679 /* Elements of the coefficient vector where neither k = 1 or i = 1
680 are not calculated if sf->xterms = no. */
681 if (sf->xterms == TNX_XNONE)
682 xorder = 1;
683 else if (sf->xterms == TNX_XHALF) {
684 if ((i + 1 + sf->xorder + 1) > maxorder)
685 xorder = xorder - 1;
686 }
687 }
688
689 return (sum);
690 }
691
692
693 /* TNX_GSCOEFF -- procedure to fetch the number and magnitude of the coefficients
694 * if the sf->xterms = wf_xbi (yes) then the number of coefficients will be
695 * (sf->xorder * sf->yorder); if wf_xterms is wf_xtri then the number
696 * of coefficients will be (sf->xorder * sf->yorder - order *
697 * (order - 1) / 2) where order is the minimum of the x and yorders; if
698 * sf->xterms = TNX_XNONE then the number of coefficients will be
699 * (sf->xorder + sf->yorder - 1).
700 */
701
702 int
wf_gscoeff(sf,coeff)703 wf_gscoeff (sf, coeff)
704
705 struct IRAFsurface *sf; /* pointer to the surface fitting descriptor */
706 double *coeff; /* the coefficients of the fit */
707
708 {
709 int ncoeff; /* the number of coefficients */
710 int i;
711
712 /* Exctract coefficients from data structure and calculate their number */
713 ncoeff = sf->ncoeff;
714 for (i = 0; i < ncoeff; i++)
715 coeff[i] = sf->coeff[i];
716 return (ncoeff);
717 }
718
719
720 static double *coeff = NULL;
721 static int nbcoeff = 0;
722
723 /* wf_gsder -- procedure to calculate a new surface which is a derivative of
724 * the input surface.
725 */
726
727 double
wf_gsder(sf1,x,y,nxd,nyd)728 wf_gsder (sf1, x, y, nxd, nyd)
729
730 struct IRAFsurface *sf1; /* pointer to the previous surface */
731 double x; /* x values */
732 double y; /* y values */
733 int nxd, nyd; /* order of the derivatives in x and y */
734 {
735 int nxder, nyder, i, j, k, nbytes;
736 int order, maxorder1, maxorder2, nmove1, nmove2;
737 struct IRAFsurface *sf2 = 0;
738 double *ptr1, *ptr2;
739 double zfit, norm;
740 double wf_gseval();
741
742 if (sf1 == NULL)
743 return (0.0);
744
745 if (nxd < 0 || nyd < 0) {
746 fprintf (stderr, "TNX_GSDER: order of derivatives cannot be < 0\n");
747 return (0.0);
748 }
749
750 if (nxd == 0 && nyd == 0) {
751 zfit = wf_gseval (sf1, x, y);
752 return (zfit);
753 }
754
755 /* Allocate space for new surface */
756 sf2 = (struct IRAFsurface *) malloc (sizeof (struct IRAFsurface));
757
758 /* Check the order of the derivatives */
759 nxder = MIN (nxd, sf1->xorder - 1);
760 nyder = MIN (nyd, sf1->yorder - 1);
761
762 /* Set up new surface */
763 sf2->type = sf1->type;
764
765 /* Set the derivative surface parameters */
766 if (sf2->type == TNX_LEGENDRE ||
767 sf2->type == TNX_CHEBYSHEV ||
768 sf2->type == TNX_POLYNOMIAL) {
769
770 sf2->xterms = sf1->xterms;
771
772 /* Find the order of the new surface */
773 switch (sf2->xterms) {
774 case TNX_XNONE:
775 if (nxder > 0 && nyder > 0) {
776 sf2->xorder = 1;
777 sf2->yorder = 1;
778 sf2->ncoeff = 1;
779 }
780 else if (nxder > 0) {
781 sf2->xorder = MAX (1, sf1->xorder - nxder);
782 sf2->yorder = 1;
783 sf2->ncoeff = sf2->xorder;
784 }
785 else if (nyder > 0) {
786 sf2->xorder = 1;
787 sf2->yorder = MAX (1, sf1->yorder - nyder);
788 sf2->ncoeff = sf2->yorder;
789 }
790 break;
791
792 case TNX_XHALF:
793 maxorder1 = MAX (sf1->xorder+1, sf1->yorder+1);
794 order = MAX(1, MIN(maxorder1-1-nyder-nxder,sf1->xorder-nxder));
795 sf2->xorder = order;
796 order = MAX(1, MIN(maxorder1-1-nyder-nxder,sf1->yorder-nyder));
797 sf2->yorder = order;
798 order = MIN (sf2->xorder, sf2->yorder);
799 sf2->ncoeff = sf2->xorder * sf2->yorder - (order*(order-1)/2);
800 break;
801
802 default:
803 sf2->xorder = MAX (1, sf1->xorder - nxder);
804 sf2->yorder = MAX (1, sf1->yorder - nyder);
805 sf2->ncoeff = sf2->xorder * sf2->yorder;
806 }
807
808 /* define the data limits */
809 sf2->xrange = sf1->xrange;
810 sf2->xmaxmin = sf1->xmaxmin;
811 sf2->yrange = sf1->yrange;
812 sf2->ymaxmin = sf1->ymaxmin;
813 }
814
815 else {
816 fprintf (stderr, "TNX_GSDER: unknown surface type %d\n", sf2->type);
817 return (0.0);
818 }
819
820 /* Allocate space for coefficients and basis functions */
821 nbytes = sf2->ncoeff * sizeof(double);
822 sf2->coeff = (double *) malloc (nbytes);
823 nbytes = sf2->xorder * sizeof(double);
824 sf2->xbasis = (double *) malloc (nbytes);
825 nbytes = sf2->yorder * sizeof(double);
826 sf2->ybasis = (double *) malloc (nbytes);
827
828 /* Get coefficients */
829 nbytes = sf1->ncoeff * sizeof(double);
830 if (nbytes > nbcoeff) {
831 if (nbcoeff > 0)
832 coeff = (double *) realloc (coeff, nbytes);
833 else
834 coeff = (double *) malloc (nbytes);
835 nbcoeff = nbytes;
836 }
837 (void) wf_gscoeff (sf1, coeff);
838
839 /* Compute the new coefficients */
840 switch (sf2->xterms) {
841 case TNX_XFULL:
842 ptr2 = sf2->coeff + (sf2->yorder - 1) * sf2->xorder;
843 ptr1 = coeff + (sf1->yorder - 1) * sf1->xorder;
844 for (i = sf1->yorder - 1; i >= nyder; i--) {
845 for (j = i; j >= i-nyder+1; j--) {
846 for (k = 0; k < sf2->xorder; k++)
847 ptr1[nxder+k] = ptr1[nxder+k] * (double)(j);
848 }
849 for (j = sf1->xorder; j >= nxder+1; j--) {
850 for (k = j; k >= j-nxder+1; k--)
851 ptr1[j-1] = ptr1[j-1] * (double)(k - 1);
852 }
853 for (j = 0; j < sf2->xorder; j++)
854 ptr2[j] = ptr1[nxder+j];
855 ptr2 = ptr2 - sf2->xorder;
856 ptr1 = ptr1 - sf1->xorder;
857 }
858 break;
859
860 case TNX_XHALF:
861 maxorder1 = MAX (sf1->xorder + 1, sf1->yorder + 1);
862 maxorder2 = MAX (sf2->xorder + 1, sf2->yorder + 1);
863 ptr2 = sf2->coeff + sf2->ncoeff;
864 ptr1 = coeff + sf1->ncoeff;
865 for (i = sf1->yorder; i >= nyder+1; i--) {
866 nmove1 = MAX (0, MIN (maxorder1 - i, sf1->xorder));
867 nmove2 = MAX (0, MIN (maxorder2 - i + nyder, sf2->xorder));
868 ptr1 = ptr1 - nmove1;
869 ptr2 = ptr2 - nmove2;
870 for (j = i; j > i - nyder + 1; j--) {
871 for (k = 0; k < nmove2; k++)
872 ptr1[nxder+k] = ptr1[nxder+k] * (double)(j-1);
873 }
874 for (j = nmove1; j >= nxder+1; j--) {
875 for (k = j; k >= j-nxder+1; k--)
876 ptr1[j-1] = ptr1[j-1] * (double)(k - 1);
877 }
878 for (j = 0; j < nmove2; j++)
879 ptr2[j] = ptr1[nxder+j];
880 }
881 break;
882
883 default:
884 if (nxder > 0 && nyder > 0)
885 sf2->coeff[0] = 0.0;
886
887 else if (nxder > 0) {
888 ptr1 = coeff;
889 ptr2 = sf2->coeff + sf2->ncoeff - 1;
890 for (j = sf1->xorder; j >= nxder+1; j--) {
891 for (k = j; k >= j - nxder + 1; k--)
892 ptr1[j-1] = ptr1[j-1] * (double)(k - 1);
893 ptr2[0] = ptr1[j-1];
894 ptr2 = ptr2 - 1;
895 }
896 }
897
898 else if (nyder > 0) {
899 ptr1 = coeff + sf1->ncoeff - 1;
900 ptr2 = sf2->coeff;
901 for (i = sf1->yorder; i >= nyder + 1; i--) {
902 for (j = i; j >= i - nyder + 1; j--)
903 *ptr1 = *ptr1 * (double)(j - 1);
904 ptr1 = ptr1 - 1;
905 }
906 for (i = 0; i < sf2->ncoeff; i++)
907 ptr2[i] = ptr1[i+1];
908 }
909 }
910
911 /* evaluate the derivatives */
912 zfit = wf_gseval (sf2, x, y);
913
914 /* normalize */
915 if (sf2->type != TNX_POLYNOMIAL) {
916 norm = pow (sf2->xrange, (double)nxder) *
917 pow (sf2->yrange, (double)nyder);
918 zfit = norm * zfit;
919 }
920
921 /* free the space */
922 wf_gsclose (sf2);
923
924 return (zfit);
925 }
926
927
928 /* wf_gsrestore -- procedure to restore the surface fit encoded in the
929 image header as a list of double precision parameters and coefficients
930 to the surface descriptor for use by the evaluating routines. the
931 surface parameters, surface type, xorder (or number of polynomial
932 terms in x), yorder (or number of polynomial terms in y), xterms,
933 xmin, xmax and ymin and ymax, are stored in the first eight elements
934 of the double array fit, followed by the wf->ncoeff surface coefficients.
935 */
936
937 struct IRAFsurface *
wf_gsrestore(fit)938 wf_gsrestore (fit)
939
940 double *fit; /* array containing the surface parameters
941 and coefficients */
942 {
943 struct IRAFsurface *sf; /* surface descriptor */
944 int surface_type, xorder, yorder, order, i;
945 double xmin, xmax, ymin, ymax;
946
947 xorder = (int) (fit[TNX_SAVEXORDER] + 0.5);
948 if (xorder < 1) {
949 fprintf (stderr, "wf_gsrestore: illegal x order %d\n", xorder);
950 return (NULL);
951 }
952
953 yorder = (int) (fit[TNX_SAVEYORDER] + 0.5);
954 if (yorder < 1) {
955 fprintf (stderr, "wf_gsrestore: illegal y order %d\n", yorder);
956 return (NULL);
957 }
958
959 xmin = fit[TNX_SAVEXMIN];
960 xmax = fit[TNX_SAVEXMAX];
961 if (xmax <= xmin) {
962 fprintf (stderr, "wf_gsrestore: illegal x range %f-%f\n",xmin,xmax);
963 return (NULL);
964 }
965 ymin = fit[TNX_SAVEYMIN];
966 ymax = fit[TNX_SAVEYMAX];
967 if (ymax <= ymin) {
968 fprintf (stderr, "wf_gsrestore: illegal y range %f-%f\n",ymin,ymax);
969 return (NULL);
970 }
971
972 /* Set surface type dependent surface descriptor parameters */
973 surface_type = (int) (fit[TNX_SAVETYPE] + 0.5);
974
975 if (surface_type == TNX_LEGENDRE ||
976 surface_type == TNX_CHEBYSHEV ||
977 surface_type == TNX_POLYNOMIAL) {
978
979 /* allocate space for the surface descriptor */
980 sf = (struct IRAFsurface *) malloc (sizeof (struct IRAFsurface));
981 sf->xorder = xorder;
982 sf->xrange = 2.0 / (xmax - xmin);
983 sf->xmaxmin = - (xmax + xmin) / 2.0;
984 sf->yorder = yorder;
985 sf->yrange = 2.0 / (ymax - ymin);
986 sf->ymaxmin = - (ymax + ymin) / 2.0;
987 sf->xterms = fit[TNX_SAVEXTERMS];
988 switch (sf->xterms) {
989 case TNX_XNONE:
990 sf->ncoeff = sf->xorder + sf->yorder - 1;
991 break;
992 case TNX_XHALF:
993 order = MIN (xorder, yorder);
994 sf->ncoeff = sf->xorder * sf->yorder - order * (order-1) / 2;
995 break;
996 case TNX_XFULL:
997 sf->ncoeff = sf->xorder * sf->yorder;
998 break;
999 }
1000 }
1001 else {
1002 fprintf (stderr, "wf_gsrestore: unknown surface type %d\n", surface_type);
1003 return (NULL);
1004 }
1005
1006 /* Set remaining curve parameters */
1007 sf->type = surface_type;
1008
1009 /* Restore coefficient array */
1010 sf->coeff = (double *) malloc (sf->ncoeff*sizeof (double));
1011 for (i = 0; i < sf->ncoeff; i++)
1012 sf->coeff[i] = fit[TNX_SAVECOEFF+i];
1013
1014 /* Allocate space for basis vectors */
1015 sf->xbasis = (double *) malloc (sf->xorder*sizeof (double));
1016 sf->ybasis = (double *) malloc (sf->yorder*sizeof (double));
1017
1018 return (sf);
1019 }
1020
1021
1022 /* wf_gsb1pol -- procedure to evaluate all the non-zero polynomial functions
1023 for a single point and given order. */
1024
1025 static void
wf_gsb1pol(x,order,basis)1026 wf_gsb1pol (x, order, basis)
1027
1028 double x; /*i data point */
1029 int order; /*i order of polynomial, order = 1, constant */
1030 double *basis; /*o basis functions */
1031 {
1032 int i;
1033
1034 basis[0] = 1.0;
1035 if (order == 1)
1036 return;
1037
1038 basis[1] = x;
1039 if (order == 2)
1040 return;
1041
1042 for (i = 2; i < order; i++)
1043 basis[i] = x * basis[i-1];
1044
1045 return;
1046 }
1047
1048
1049 /* wf_gsb1leg -- procedure to evaluate all the non-zero legendre functions for
1050 a single point and given order. */
1051
1052 static void
wf_gsb1leg(x,order,k1,k2,basis)1053 wf_gsb1leg (x, order, k1, k2, basis)
1054
1055 double x; /*i data point */
1056 int order; /*i order of polynomial, order = 1, constant */
1057 double k1, k2; /*i normalizing constants */
1058 double *basis; /*o basis functions */
1059 {
1060 int i;
1061 double ri, xnorm;
1062
1063 basis[0] = 1.0;
1064 if (order == 1)
1065 return;
1066
1067 xnorm = (x + k1) * k2 ;
1068 basis[1] = xnorm;
1069 if (order == 2)
1070 return;
1071
1072 for (i = 2; i < order; i++) {
1073 ri = i;
1074 basis[i] = ((2.0 * ri - 1.0) * xnorm * basis[i-1] -
1075 (ri - 1.0) * basis[i-2]) / ri;
1076 }
1077
1078 return;
1079 }
1080
1081
1082 /* wf_gsb1cheb -- procedure to evaluate all the non-zero chebyshev function
1083 coefficients for a given x and order. */
1084
1085 static void
wf_gsb1cheb(x,order,k1,k2,basis)1086 wf_gsb1cheb (x, order, k1, k2, basis)
1087
1088 double x; /*i number of data points */
1089 int order; /*i order of polynomial, 1 is a constant */
1090 double k1, k2; /*i normalizing constants */
1091 double *basis; /*o array of basis functions */
1092 {
1093 int i;
1094 double xnorm;
1095
1096 basis[0] = 1.0;
1097 if (order == 1)
1098 return;
1099
1100 xnorm = (x + k1) * k2;
1101 basis[1] = xnorm;
1102 if (order == 2)
1103 return;
1104
1105 for (i = 2; i < order; i++)
1106 basis[i] = 2. * xnorm * basis[i-1] - basis[i-2];
1107
1108 return;
1109 }
1110
1111 /* Set surface polynomial from arguments */
1112
1113 int
tnxpset(wcs,xorder,yorder,xterms,coeff)1114 tnxpset (wcs, xorder, yorder, xterms, coeff)
1115
1116 struct WorldCoor *wcs; /* World coordinate system structure */
1117 int xorder; /* Number of x coefficients (same for x and y) */
1118 int yorder; /* Number of y coefficients (same for x and y) */
1119 int xterms; /* Number of xy coefficients (same for x and y) */
1120 double *coeff; /* Plate fit coefficients */
1121
1122 {
1123 double *ycoeff;
1124 struct IRAFsurface *wf_gspset ();
1125
1126 wcs->prjcode = WCS_TNX;
1127
1128 wcs->lngcor = wf_gspset (xorder, yorder, xterms, coeff);
1129 ycoeff = coeff + wcs->lngcor->ncoeff;
1130 wcs->latcor = wf_gspset (xorder, yorder, xterms, ycoeff);
1131
1132 return 0;
1133 }
1134
1135
1136 /* wf_gspset -- procedure to set the surface descriptor for use by the
1137 evaluating routines. from arguments. The surface parameters are
1138 surface type, xorder (number of polynomial terms in x), yorder (number
1139 of polynomial terms in y), xterms, and the surface coefficients.
1140 */
1141
1142 struct IRAFsurface *
wf_gspset(xorder,yorder,xterms,coeff)1143 wf_gspset (xorder, yorder, xterms, coeff)
1144
1145 int xorder;
1146 int yorder;
1147 int xterms;
1148 double *coeff;
1149 {
1150 struct IRAFsurface *sf; /* surface descriptor */
1151 int surface_type, order, i;
1152 double xmin, xmax;
1153 double ymin, ymax;
1154
1155 surface_type = TNX_POLYNOMIAL;
1156 xmin = 0.0;
1157 xmax = 0.0;
1158 ymin = 0.0;
1159 ymax = 0.0;
1160
1161 if (surface_type == TNX_LEGENDRE ||
1162 surface_type == TNX_CHEBYSHEV ||
1163 surface_type == TNX_POLYNOMIAL) {
1164
1165 /* allocate space for the surface descriptor */
1166 sf = (struct IRAFsurface *) malloc (sizeof (struct IRAFsurface));
1167 sf->xorder = xorder;
1168 sf->xrange = 2.0 / (xmax - xmin);
1169 sf->xmaxmin = -(xmax + xmin) / 2.0;
1170 sf->yorder = yorder;
1171 sf->yrange = 2.0 / (ymax - ymin);
1172 sf->ymaxmin = - (ymax + ymin) / 2.0;
1173 sf->xterms = xterms;
1174 switch (sf->xterms) {
1175 case TNX_XNONE:
1176 sf->ncoeff = sf->xorder + sf->yorder - 1;
1177 break;
1178 case TNX_XHALF:
1179 order = MIN (xorder, yorder);
1180 sf->ncoeff = sf->xorder * sf->yorder - order * (order-1) / 2;
1181 break;
1182 case TNX_XFULL:
1183 sf->ncoeff = sf->xorder * sf->yorder;
1184 break;
1185 }
1186 }
1187 else {
1188 fprintf (stderr, "TNX_GSSET: unknown surface type %d\n", surface_type);
1189 return (NULL);
1190 }
1191
1192 /* Set remaining curve parameters */
1193 sf->type = surface_type;
1194
1195 /* Restore coefficient array */
1196 sf->coeff = (double *) malloc (sf->ncoeff*sizeof (double));
1197 for (i = 0; i < sf->ncoeff; i++)
1198 sf->coeff[i] = coeff[i];
1199
1200 /* Allocate space for basis vectors */
1201 sf->xbasis = (double *) malloc (sf->xorder*sizeof (double));
1202 sf->ybasis = (double *) malloc (sf->yorder*sizeof (double));
1203
1204 return (sf);
1205 }
1206
1207 /* Mar 26 1998 New subroutines, translated from SPP
1208 * Apr 28 1998 Change all local flags to TNX_* and projection flag to WCS_TNX
1209 * May 11 1998 Fix use of pole longitude default
1210 * Sep 4 1998 Fix missed assignment in tnxpos from Allen Harris, SAO
1211 * Sep 10 1998 Fix bugs in tnxpix()
1212 * Sep 10 1998 Fix missed assignment in tnxpix from Allen Harris, SAO
1213 *
1214 * Oct 22 1999 Drop unused variables, fix case statements after lint
1215 * Dec 10 1999 Fix bug in gsder() which failed to allocate enough memory
1216 * Dec 10 1999 Compute wcs->rot using wcsrotset() in tnxinit()
1217 *
1218 * Feb 14 2001 Fixed off-by-one bug in legendre evaluation (Mike Jarvis)
1219 *
1220 * Apr 11 2002 Fix bug when .-terminated substring in wf_gsopen()
1221 * Apr 29 2002 Clean up code
1222 * Jun 26 2002 Increase size of WAT strings from 500 to 2000
1223 *
1224 * Jun 27 2005 Drop unused arguments k1 and k2 from wf_gsb1pol()
1225 *
1226 * Jan 8 2007 Drop unused variable ncoeff in wf_gsder()
1227 * Jan 9 2007 Declare header const char in tnxinit()
1228 * Apr 3 2007 Fix offsets to hit last cooefficient in wf_gsder()
1229 *
1230 * Sep 5 2008 Fix wf_gseval() call in tnxpos() so unmodified x and y are used
1231 * Sep 9 2008 Fix loop in TNX_XFULL section of wf_gsder()
1232 * (last two bugs found by Ed Los)
1233 * Sep 17 2008 Fix tnxpos for null correction case (fix by Ed Los)
1234 */
1235