1 /*============================================================================
2 *
3 * WCSLIB - an implementation of the FITS WCS proposal.
4 * Copyright (C) 1995-1999, Mark Calabretta
5 *
6 * This library is free software; you can redistribute it and/or modify it
7 * under the terms of the GNU Library General Public License as published
8 * by the Free Software Foundation; either version 2 of the License, or (at
9 * your option) any later version.
10 *
11 * This library is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library
14 * General Public License for more details.
15 *
16 * You should have received a copy of the GNU Library General Public License
17 * along with this library; if not, write to the Free Software Foundation,
18 * Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
19 *
20 * Correspondence concerning WCSLIB may be directed to:
21 * Internet email: mcalabre@atnf.csiro.au
22 * Postal address: Dr. Mark Calabretta,
23 * Australia Telescope National Facility,
24 * P.O. Box 76,
25 * Epping, NSW, 2121,
26 * AUSTRALIA
27 *
28 *=============================================================================
29 *
30 * C implementation of the spherical map projections recognized by the FITS
31 * "World Coordinate System" (WCS) convention.
32 *
33 * Summary of routines
34 * -------------------
35 * Each projection is implemented via separate functions for the forward,
36 * *fwd(), and reverse, *rev(), transformation.
37 *
38 * Initialization routines, *set(), compute intermediate values from the
39 * projection parameters but need not be called explicitly - see the
40 * explanation of prj.flag below.
41 *
42 * azpset azpfwd azprev AZP: zenithal/azimuthal perspective
43 * tanset tanfwd tanrev TAN: gnomonic
44 * sinset sinfwd sinrev SIN: orthographic/synthesis
45 * stgset stgfwd stgrev STG: stereographic
46 * arcset arcfwd arcrev ARC: zenithal/azimuthal equidistant
47 * zpnset zpnfwd zpnrev ZPN: zenithal/azimuthal polynomial
48 * zeaset zeafwd zearev ZEA: zenithal/azimuthal equal area
49 * airset airfwd airrev AIR: Airy
50 * cypset cypfwd cyprev CYP: cylindrical perspective
51 * carset carfwd carrev CAR: Cartesian
52 * merset merfwd merrev MER: Mercator
53 * ceaset ceafwd cearev CEA: cylindrical equal area
54 * copset copfwd coprev COP: conic perspective
55 * codset codfwd codrev COD: conic equidistant
56 * coeset coefwd coerev COE: conic equal area
57 * cooset coofwd coorev COO: conic orthomorphic
58 * bonset bonfwd bonrev BON: Bonne
59 * pcoset pcofwd pcorev PCO: polyconic
60 * glsset glsfwd glsrev GLS: Sanson-Flamsteed (global sinusoidal)
61 * parset parfwd parrev PAR: parabolic
62 * aitset aitfwd aitrev AIT: Hammer-Aitoff
63 * molset molfwd molrev MOL: Mollweide
64 * cscset cscfwd cscrev CSC: COBE quadrilateralized spherical cube
65 * qscset qscfwd qscrev QSC: quadrilateralized spherical cube
66 * tscset tscfwd tscrev TSC: tangential spherical cube
67 * tnxset tnxfwd tnxrev TNX: IRAF's gnomonic
68 *
69 *
70 * Initialization routine; *set()
71 * ------------------------------
72 * Initializes members of a prjprm data structure which hold intermediate
73 * values. Note that this routine need not be called directly; it will be
74 * invoked by prjfwd() and prjrev() if the "flag" structure member is
75 * anything other than a predefined magic value.
76 *
77 * Given and/or returned:
78 * prj prjprm* Projection parameters (see below).
79 *
80 * Function return value:
81 * int Error status
82 * 0: Success.
83 * 1: Invalid projection parameters.
84 *
85 * Forward transformation; *fwd()
86 * -----------------------------
87 * Compute (x,y) coordinates in the plane of projection from native spherical
88 * coordinates (phi,theta).
89 *
90 * Given:
91 * phi, const double
92 * theta Longitude and latitude of the projected point in
93 * native spherical coordinates, in degrees.
94 *
95 * Given and returned:
96 * prj prjprm* Projection parameters (see below).
97 *
98 * Returned:
99 * x,y double* Projected coordinates.
100 *
101 * Function return value:
102 * int Error status
103 * 0: Success.
104 * 1: Invalid projection parameters.
105 * 2: Invalid value of (phi,theta).
106 *
107 * Reverse transformation; *rev()
108 * -----------------------------
109 * Compute native spherical coordinates (phi,theta) from (x,y) coordinates in
110 * the plane of projection.
111 *
112 * Given:
113 * x,y const double
114 * Projected coordinates.
115 *
116 * Given and returned:
117 * prj prjprm* Projection parameters (see below).
118 *
119 * Returned:
120 * phi, double* Longitude and latitude of the projected point in
121 * theta native spherical coordinates, in degrees.
122 *
123 * Function return value:
124 * int Error status
125 * 0: Success.
126 * 1: Invalid projection parameters.
127 * 2: Invalid value of (x,y).
128 *
129 * Projection parameters
130 * ---------------------
131 * The prjprm struct consists of the following:
132 *
133 * int flag
134 * This flag must be set to zero whenever any of p[10] or r0 are set
135 * or changed. This signals the initialization routine to recompute
136 * intermediaries. flag may also be set to -1 to disable strict bounds
137 * checking for the AZP, TAN, SIN, ZPN, and COP projections.
138 * double r0
139 * r0; The radius of the generating sphere for the projection, a linear
140 * scaling parameter. If this is zero, it will be reset to the default
141 * value of 180/pi (the value for FITS WCS).
142 * double p[10]
143 * The first 10 elements contain projection parameters which correspond
144 * to the PROJPn keywords in FITS, so p[0] is PROJP0, and p[9] is
145 * PROJP9. Many projections use p[1] (PROJP1) and some also use p[2]
146 * (PROJP2). ZPN is the only projection which uses any of the others.
147 *
148 * The remaining members of the prjprm struct are maintained by the
149 * initialization routines and should not be modified. This is done for the
150 * sake of efficiency and to allow an arbitrary number of contexts to be
151 * maintained simultaneously.
152 *
153 * int n
154 * double w[10]
155 * Intermediate values derived from the projection parameters.
156 *
157 * Usage of the p[] array as it applies to each projection is described in
158 * the prologue to each trio of projection routines.
159 *
160 * Argument checking
161 * -----------------
162 * Forward routines:
163 *
164 * The values of phi and theta (the native longitude and latitude)
165 * normally lie in the range [-180,180] for phi, and [-90,90] for theta.
166 * However, all forward projections will accept any value of phi and will
167 * not normalize it.
168 *
169 * The forward projection routines do not explicitly check that theta lies
170 * within the range [-90,90]. They do check for any value of theta which
171 * produces an invalid argument to the projection equations (e.g. leading
172 * to division by zero). The forward routines for AZP, TAN, SIN, ZPN, and
173 * COP also return error 2 if (phi,theta) corresponds to the overlapped
174 * (far) side of the projection but also return the corresponding value of
175 * (x,y). This strict bounds checking may be relaxed by setting prj->flag
176 * to -1 (rather than 0) when these projections are initialized.
177 *
178 * Reverse routines:
179 *
180 * Error checking on the projected coordinates (x,y) is limited to that
181 * required to ascertain whether a solution exists. Where a solution does
182 * exist no check is made that the value of phi and theta obtained lie
183 * within the ranges [-180,180] for phi, and [-90,90] for theta.
184 *
185 * Accuracy
186 * --------
187 * Closure to a precision of at least 1E-10 degree of longitude and latitude
188 * has been verified for typical projection parameters on the 1 degree grid
189 * of native longitude and latitude (to within 5 degrees of any latitude
190 * where the projection may diverge).
191 *
192 * Author: Mark Calabretta, Australia Telescope National Facility
193 * IRAF's TNX added by E.Bertin 2000/08/23
194 * $Id: proj.c,v 1.1.1.1 2004/01/04 21:33:26 bertin Exp $
195 *===========================================================================*/
196
197 #ifdef HAVE_CONFIG_H
198 #include "config.h"
199 #endif
200
201 #ifdef HAVE_MATHIMF_H
202 #include <mathimf.h>
203 #else
204 #include <math.h>
205 #endif
206 #include <stdio.h>
207 #include <stdlib.h>
208 #include "poly.h"
209 #include "proj.h"
210 #include "tnx.h"
211 #include "wcsmath.h"
212 #include "wcstrig.h"
213
214 /* Map error number to error message for each function. */
215 const char *prjset_errmsg[] = {
216 0,
217 "Invalid projection parameters"};
218
219 const char *prjfwd_errmsg[] = {
220 0,
221 "Invalid projection parameters",
222 "Invalid value of (phi,theta)"};
223
224 const char *prjrev_errmsg[] = {
225 0,
226 "Invalid projection parameters",
227 "Invalid value of (x,y)"};
228
229 #define wcs_copysign(X, Y) ((Y) < 0.0 ? -fabs(X) : fabs(X))
230
231 /*============================================================================
232 * AZP: zenithal/azimuthal perspective projection.
233 *
234 * Given:
235 * prj->p[1] AZP distance parameter, mu in units of r0.
236 *
237 * Given and/or returned:
238 * prj->r0 r0; reset to 180/pi if 0.
239 * prj->w[0] r0*(mu+1)
240 * prj->w[1] 1/prj->w[0]
241 * prj->w[2] Boundary parameter, -mu for |mu| <= 1,
242 * -1/mu for |mu| >= 1.
243 *===========================================================================*/
244
azpset(prj)245 int azpset(prj)
246
247 struct prjprm *prj;
248
249 {
250 if (prj->r0 == 0.0) prj->r0 = R2D;
251
252 prj->w[0] = prj->r0*(prj->p[1] + 1.0);
253 if (prj->w[0] == 0.0) {
254 return 1;
255 }
256
257 prj->w[1] = 1.0/prj->w[0];
258 if (fabs(prj->p[1]) <= 1.0) {
259 prj->w[2] = -prj->p[1];
260 } else {
261 prj->w[2] = -1.0/prj->p[1];
262 }
263
264 if (prj->flag == -1) {
265 prj->flag = -PRJSET;
266 } else {
267 prj->flag = PRJSET;
268 }
269
270 return 0;
271 }
272
273 /*--------------------------------------------------------------------------*/
274
azpfwd(phi,theta,prj,x,y)275 int azpfwd(phi, theta, prj, x, y)
276
277 const double phi, theta;
278 struct prjprm *prj;
279 double *x, *y;
280
281 {
282 double r, s, sthe;
283
284 if (abs(prj->flag) != PRJSET) {
285 if (azpset(prj)) return 1;
286 }
287
288 sthe = wcs_sind(theta);
289
290 s = prj->p[1] + sthe;
291 if (s == 0.0) {
292 return 2;
293 }
294
295 r = prj->w[0]*wcs_cosd(theta)/s;
296 *x = r*wcs_sind(phi);
297 *y = -r*wcs_cosd(phi);
298
299 if (prj->flag == PRJSET && sthe < prj->w[2]) {
300 return 2;
301 }
302
303 return 0;
304 }
305
306 /*--------------------------------------------------------------------------*/
307
azprev(x,y,prj,phi,theta)308 int azprev(x, y, prj, phi, theta)
309
310 const double x, y;
311 struct prjprm *prj;
312 double *phi, *theta;
313
314 {
315 double r, rho, s;
316 const double tol = 1.0e-13;
317
318 if (abs(prj->flag) != PRJSET) {
319 if (azpset(prj)) return 1;
320 }
321
322 r = sqrt(x*x + y*y);
323 if (r == 0.0) {
324 *phi = 0.0;
325 } else {
326 *phi = wcs_atan2d(x, -y);
327 }
328
329 rho = r*prj->w[1];
330 s = rho*prj->p[1]/sqrt(rho*rho+1.0);
331 if (fabs(s) > 1.0) {
332 if (fabs(s) > 1.0+tol) {
333 return 2;
334 }
335 *theta = wcs_atan2d(1.0,rho) - wcs_copysign(90.0,s);
336 } else {
337 *theta = wcs_atan2d(1.0,rho) - wcs_asind(s);
338 }
339
340 return 0;
341 }
342
343 /*============================================================================
344 * TAN: gnomonic projection.
345 *
346 * Given and/or returned:
347 * prj->r0 r0; reset to 180/pi if 0.
348 *===========================================================================*/
349
tanset(prj)350 int tanset(prj)
351
352 struct prjprm *prj;
353
354 {
355 int k;
356
357 if (prj->r0 == 0.0) prj->r0 = R2D;
358
359 if (prj->flag == -1) {
360 prj->flag = -PRJSET;
361 } else {
362 prj->flag = PRJSET;
363 }
364
365 for (k = 99; k >= 0 && prj->p[k] == 0.0 && prj->p[k+100] == 0.0; k--);
366 if (k < 0)
367 k = 0;
368
369 prj->n = k;
370
371 return 0;
372 }
373
374 /*--------------------------------------------------------------------------*/
375
tanfwd(phi,theta,prj,x,y)376 int tanfwd(phi, theta, prj, x, y)
377
378 const double phi, theta;
379 struct prjprm *prj;
380 double *x, *y;
381
382 {
383 double r, s, xp[2];
384
385 if (abs(prj->flag) != PRJSET) {
386 if(tanset(prj)) return 1;
387 }
388
389 s = wcs_sind(theta);
390 if (s == 0.0) return 2;
391
392 r = prj->r0*wcs_cosd(theta)/s;
393 xp[0] = r*wcs_sind(phi);
394 xp[1] = -r*wcs_cosd(phi);
395 *x = prj->inv_x? poly_func(prj->inv_x, xp) : xp[0];
396 *y = prj->inv_y? poly_func(prj->inv_y, xp) : xp[1];
397
398 if (prj->flag == PRJSET && s < 0.0) {
399 return 2;
400 }
401
402 return 0;
403 }
404
405 /*--------------------------------------------------------------------------*/
406
tanrev(x,y,prj,phi,theta)407 int tanrev(x, y, prj, phi, theta)
408
409 const double x, y;
410 struct prjprm *prj;
411 double *phi, *theta;
412
413 {
414 double xp,yp, rp;
415
416 if (abs(prj->flag) != PRJSET) {
417 if (tanset(prj)) return 1;
418 }
419
420 if (prj->n)
421 raw_to_pv(prj, x,y, &xp, &yp);
422 else
423 {
424 xp = x;
425 yp = y;
426 }
427 rp = sqrt(xp*xp+yp*yp);
428 if (rp == 0.0) {
429 *phi = 0.0;
430 } else {
431 *phi = wcs_atan2d(xp, -yp);
432 }
433 *theta = wcs_atan2d(prj->r0, rp);
434
435 return 0;
436 }
437
438 /*============================================================================
439 * SIN: orthographic/synthesis projection.
440 *
441 * Given:
442 * prj->p[1:2] SIN obliqueness parameters, alpha and beta.
443 *
444 * Given and/or returned:
445 * prj->r0 r0; reset to 180/pi if 0.
446 * prj->w[0] 1/r0
447 * prj->w[1] alpha**2 + beta**2
448 * prj->w[2] 2*(alpha**2 + beta**2)
449 * prj->w[3] 2*(alpha**2 + beta**2 + 1)
450 * prj->w[4] alpha**2 + beta**2 - 1
451 *===========================================================================*/
452
sinset(prj)453 int sinset(prj)
454
455 struct prjprm *prj;
456
457 {
458 if (prj->r0 == 0.0) prj->r0 = R2D;
459
460 prj->w[0] = 1.0/prj->r0;
461 prj->w[1] = prj->p[1]*prj->p[1] + prj->p[2]*prj->p[2];
462 prj->w[2] = 2.0*prj->w[1];
463 prj->w[3] = prj->w[2] + 2.0;
464 prj->w[4] = prj->w[1] - 1.0;
465
466 if (prj->flag == -1) {
467 prj->flag = -PRJSET;
468 } else {
469 prj->flag = PRJSET;
470 }
471
472 return 0;
473 }
474
475 /*--------------------------------------------------------------------------*/
476
sinfwd(phi,theta,prj,x,y)477 int sinfwd(phi, theta, prj, x, y)
478
479 const double phi, theta;
480 struct prjprm *prj;
481 double *x, *y;
482
483 {
484 double cphi, cthe, sphi, t, z;
485
486 if (abs(prj->flag) != PRJSET) {
487 if (sinset(prj)) return 1;
488 }
489
490 t = (90.0 - fabs(theta))*D2R;
491 if (t < 1.0e-5) {
492 if (theta > 0.0) {
493 z = -t*t/2.0;
494 } else {
495 z = -2.0 + t*t/2.0;
496 }
497 cthe = t;
498 } else {
499 z = wcs_sind(theta) - 1.0;
500 cthe = wcs_cosd(theta);
501 }
502
503 cphi = wcs_cosd(phi);
504 sphi = wcs_sind(phi);
505 *x = prj->r0*(cthe*sphi + prj->p[1]*z);
506 *y = -prj->r0*(cthe*cphi + prj->p[2]*z);
507 /* Validate this solution. */
508 if (prj->flag == PRJSET) {
509 if (prj->w[1] == 0.0) {
510 /* Orthographic projection. */
511 if (theta < 0.0) {
512 return 2;
513 }
514 } else {
515 /* "Synthesis" projection. */
516 t = wcs_atand(prj->p[1]*sphi + prj->p[2]*cphi);
517 if (theta < t) {
518 return 2;
519 }
520 }
521 }
522
523 return 0;
524 }
525
526 /*--------------------------------------------------------------------------*/
527
sinrev(x,y,prj,phi,theta)528 int sinrev (x, y, prj, phi, theta)
529
530 const double x, y;
531 struct prjprm *prj;
532 double *phi, *theta;
533
534 {
535 const double tol = 1.0e-13;
536 double a, b, c, d, r2, sth, sth1, sth2, sxy, x0, xp, y0, yp, z;
537
538 if (abs(prj->flag) != PRJSET) {
539 if (sinset(prj)) return 1;
540 }
541
542 /* Compute intermediaries. */
543 x0 = x*prj->w[0];
544 y0 = y*prj->w[0];
545 r2 = x0*x0 + y0*y0;
546
547 if (prj->w[1] == 0.0) {
548 /* Orthographic projection. */
549 if (r2 != 0.0) {
550 *phi = wcs_atan2d(x0, -y0);
551 } else {
552 *phi = 0.0;
553 }
554
555 if (r2 < 0.5) {
556 *theta = wcs_acosd(sqrt(r2));
557 } else if (r2 <= 1.0) {
558 *theta = wcs_asind(sqrt(1.0 - r2));
559 } else {
560 return 2;
561 }
562
563 } else {
564 /* "Synthesis" projection. */
565 if (r2 < 1.0e-10) {
566 /* Use small angle formula. */
567 z = -r2/2.0;
568 *theta = 90.0 - R2D*sqrt(r2/(1.0 - x0*prj->p[1] + y0*prj->p[2]));
569
570 } else {
571 sxy = 2.0*(prj->p[1]*x0 - prj->p[2]*y0);
572
573 a = prj->w[3];
574 b = -(sxy + prj->w[2]);
575 c = r2 + sxy + prj->w[4];
576 d = b*b - 2.0*a*c;
577
578 /* Check for a solution. */
579 if (d < 0.0) {
580 return 2;
581 }
582 d = sqrt(d);
583
584 /* Choose solution closest to pole. */
585 sth1 = (-b + d)/a;
586 sth2 = (-b - d)/a;
587 sth = (sth1>sth2) ? sth1 : sth2;
588 if (sth > 1.0) {
589 if (sth-1.0 < tol) {
590 sth = 1.0;
591 } else {
592 sth = (sth1<sth2) ? sth1 : sth2;
593 }
594 }
595 if (sth > 1.0 || sth < -1.0) {
596 return 2;
597 }
598
599 *theta = wcs_asind(sth);
600 z = sth - 1.0;
601 }
602
603 xp = -y0 - prj->p[2]*z;
604 yp = x0 - prj->p[1]*z;
605 if (xp == 0.0 && yp == 0.0) {
606 *phi = 0.0;
607 } else {
608 *phi = wcs_atan2d(yp,xp);
609 }
610 }
611
612 return 0;
613 }
614
615 /*============================================================================
616 * STG: stereographic projection.
617 *
618 * Given and/or returned:
619 * prj->r0 r0; reset to 180/pi if 0.
620 * prj->w[0] 2*r0
621 * prj->w[1] 1/(2*r0)
622 *===========================================================================*/
623
stgset(prj)624 int stgset(prj)
625
626 struct prjprm *prj;
627
628 {
629 if (prj->r0 == 0.0) {
630 prj->r0 = R2D;
631 prj->w[0] = 360.0/PI;
632 prj->w[1] = PI/360.0;
633 } else {
634 prj->w[0] = 2.0*prj->r0;
635 prj->w[1] = 1.0/prj->w[0];
636 }
637
638 prj->flag = PRJSET;
639 return 0;
640 }
641
642 /*--------------------------------------------------------------------------*/
643
stgfwd(phi,theta,prj,x,y)644 int stgfwd(phi, theta, prj, x, y)
645
646 const double phi, theta;
647 struct prjprm *prj;
648 double *x, *y;
649
650 {
651 double r, s;
652
653 if (prj->flag != PRJSET) {
654 if (stgset(prj)) return 1;
655 }
656
657 s = 1.0 + wcs_sind(theta);
658 if (s == 0.0) {
659 return 2;
660 }
661
662 r = prj->w[0]*wcs_cosd(theta)/s;
663 *x = r*wcs_sind(phi);
664 *y = -r*wcs_cosd(phi);
665
666 return 0;
667 }
668
669 /*--------------------------------------------------------------------------*/
670
stgrev(x,y,prj,phi,theta)671 int stgrev(x, y, prj, phi, theta)
672
673 const double x, y;
674 struct prjprm *prj;
675 double *phi, *theta;
676
677 {
678 double r;
679
680 if (prj->flag != PRJSET) {
681 if (stgset(prj)) return 1;
682 }
683
684 r = sqrt(x*x + y*y);
685 if (r == 0.0) {
686 *phi = 0.0;
687 } else {
688 *phi = wcs_atan2d(x, -y);
689 }
690 *theta = 90.0 - 2.0*wcs_atand(r*prj->w[1]);
691
692 return 0;
693 }
694
695 /*============================================================================
696 * ARC: zenithal/azimuthal equidistant projection.
697 *
698 * Given and/or returned:
699 * prj->r0 r0; reset to 180/pi if 0.
700 * prj->w[0] r0*(pi/180)
701 * prj->w[1] (180/pi)/r0
702 *===========================================================================*/
703
arcset(prj)704 int arcset(prj)
705
706 struct prjprm *prj;
707
708 {
709 if (prj->r0 == 0.0) {
710 prj->r0 = R2D;
711 prj->w[0] = 1.0;
712 prj->w[1] = 1.0;
713 } else {
714 prj->w[0] = prj->r0*D2R;
715 prj->w[1] = 1.0/prj->w[0];
716 }
717
718 prj->flag = PRJSET;
719
720 return 0;
721 }
722
723 /*--------------------------------------------------------------------------*/
724
arcfwd(phi,theta,prj,x,y)725 int arcfwd(phi, theta, prj, x, y)
726
727 const double phi, theta;
728 struct prjprm *prj;
729 double *x, *y;
730
731 {
732 double r;
733
734 if (prj->flag != PRJSET) {
735 if (arcset(prj)) return 1;
736 }
737
738 r = prj->w[0]*(90.0 - theta);
739 *x = r*wcs_sind(phi);
740 *y = -r*wcs_cosd(phi);
741
742 return 0;
743 }
744
745 /*--------------------------------------------------------------------------*/
746
arcrev(x,y,prj,phi,theta)747 int arcrev(x, y, prj, phi, theta)
748
749 const double x, y;
750 struct prjprm *prj;
751 double *phi, *theta;
752
753 {
754 double r;
755
756 if (prj->flag != PRJSET) {
757 if (arcset(prj)) return 1;
758 }
759
760 r = sqrt(x*x + y*y);
761 if (r == 0.0) {
762 *phi = 0.0;
763 } else {
764 *phi = wcs_atan2d(x, -y);
765 }
766 *theta = 90.0 - r*prj->w[1];
767
768 return 0;
769 }
770
771 /*============================================================================
772 * ZPN: zenithal/azimuthal polynomial projection.
773 *
774 * Given:
775 * prj->p[0:99] Polynomial coefficients.
776 *
777 * Given and/or returned:
778 * prj->r0 r0; reset to 180/pi if 0.
779 * prj->n Degree of the polynomial, N.
780 * prj->w[0] Co-latitude of the first point of inflection (N > 2).
781 * prj->w[1] Radius of the first point of inflection (N > 2).
782 *===========================================================================*/
783
zpnset(prj)784 int zpnset(prj)
785
786 struct prjprm *prj;
787
788 {
789 int i, j, k;
790 double d, d1, d2, r, zd, zd1, zd2;
791 const double tol = 1.0e-13;
792
793 if (prj->r0 == 0.0) prj->r0 = R2D;
794
795 /* Find the highest non-zero coefficient. */
796 for (k = 99; k >= 0 && prj->p[k] == 0.0; k--);
797 if (k < 0) return 1;
798
799 prj->n = k;
800
801 if (k >= 3) {
802 /* Find the point of inflection closest to the pole. */
803 zd1 = 0.0;
804 d1 = prj->p[1];
805 if (d1 <= 0.0) {
806 return 1;
807 }
808
809 /* Find the point where the derivative first goes negative. */
810 for (i = 0; i < 180; i++) {
811 zd2 = i*D2R;
812 d2 = 0.0;
813 for (j = k; j > 0; j--) {
814 d2 = d2*zd2 + j*prj->p[j];
815 }
816
817 if (d2 <= 0.0) break;
818 zd1 = zd2;
819 d1 = d2;
820 }
821
822 if (i == 180) {
823 /* No negative derivative -> no point of inflection. */
824 zd = PI;
825 } else {
826 /* Find where the derivative is zero. */
827 for (i = 1; i <= 10; i++) {
828 zd = zd1 - d1*(zd2-zd1)/(d2-d1);
829
830 d = 0.0;
831 for (j = k; j > 0; j--) {
832 d = d*zd + j*prj->p[j];
833 }
834
835 if (fabs(d) < tol) break;
836
837 if (d < 0.0) {
838 zd2 = zd;
839 d2 = d;
840 } else {
841 zd1 = zd;
842 d1 = d;
843 }
844 }
845 }
846
847 r = 0.0;
848 for (j = k; j >= 0; j--) {
849 r = r*zd + prj->p[j];
850 }
851 prj->w[0] = zd;
852 prj->w[1] = r;
853 }
854
855 if (prj->flag == -1) {
856 prj->flag = -PRJSET;
857 } else {
858 prj->flag = PRJSET;
859 }
860
861 return 0;
862 }
863
864 /*--------------------------------------------------------------------------*/
865
zpnfwd(phi,theta,prj,x,y)866 int zpnfwd(phi, theta, prj, x, y)
867
868 const double phi, theta;
869 struct prjprm *prj;
870 double *x, *y;
871
872 {
873 int j;
874 double r, s;
875
876 if (abs(prj->flag) != PRJSET) {
877 if (zpnset(prj)) return 1;
878 }
879
880 s = (90.0 - theta)*D2R;
881 r = 0.0;
882 for (j = prj->n; j >= 0; j--) {
883 r = r*s + prj->p[j];
884 }
885 r = prj->r0*r;
886
887 *x = r*wcs_sind(phi);
888 *y = -r*wcs_cosd(phi);
889
890 if (prj->flag == PRJSET && s > prj->w[0]) {
891 return 2;
892 }
893
894 return 0;
895 }
896
897 /*--------------------------------------------------------------------------*/
898
zpnrev(x,y,prj,phi,theta)899 int zpnrev(x, y, prj, phi, theta)
900
901 const double x, y;
902 struct prjprm *prj;
903 double *phi, *theta;
904
905 {
906 int i, j, k;
907 double a, b, c, d, lambda, r, r1, r2, rt, zd, zd1, zd2;
908 const double tol = 1.0e-13;
909
910 if (abs(prj->flag) != PRJSET) {
911 if (zpnset(prj)) return 1;
912 }
913
914 k = prj->n;
915
916 r = sqrt(x*x + y*y)/prj->r0;
917
918 if (k < 1) {
919 /* Constant - no solution. */
920 return 1;
921 } else if (k == 1) {
922 /* Linear. */
923 zd = (r - prj->p[0])/prj->p[1];
924 } else if (k == 2) {
925 /* Quadratic. */
926 a = prj->p[2];
927 b = prj->p[1];
928 c = prj->p[0] - r;
929
930 d = b*b - 4.0*a*c;
931 if (d < 0.0) {
932 return 2;
933 }
934 d = sqrt(d);
935
936 /* Choose solution closest to pole. */
937 zd1 = (-b + d)/(2.0*a);
938 zd2 = (-b - d)/(2.0*a);
939 zd = (zd1<zd2) ? zd1 : zd2;
940 if (zd < -tol) zd = (zd1>zd2) ? zd1 : zd2;
941 if (zd < 0.0) {
942 if (zd < -tol) {
943 return 2;
944 }
945 zd = 0.0;
946 } else if (zd > PI) {
947 if (zd > PI+tol) {
948 return 2;
949 }
950 zd = PI;
951 }
952 } else {
953 /* Higher order - solve iteratively. */
954 zd1 = 0.0;
955 r1 = prj->p[0];
956 zd2 = prj->w[0];
957 r2 = prj->w[1];
958
959 if (r < r1) {
960 if (r < r1-tol) {
961 return 2;
962 }
963 zd = zd1;
964 } else if (r > r2) {
965 if (r > r2+tol) {
966 return 2;
967 }
968 zd = zd2;
969 } else {
970 /* Disect the interval. */
971 for (j = 0; j < 100; j++) {
972 lambda = (r2 - r)/(r2 - r1);
973 if (lambda < 0.1) {
974 lambda = 0.1;
975 } else if (lambda > 0.9) {
976 lambda = 0.9;
977 }
978
979 zd = zd2 - lambda*(zd2 - zd1);
980
981 rt = 0.0;
982 for (i = k; i >= 0; i--) {
983 rt = (rt * zd) + prj->p[i];
984 }
985
986 if (rt < r) {
987 if (r-rt < tol) break;
988 r1 = rt;
989 zd1 = zd;
990 } else {
991 if (rt-r < tol) break;
992 r2 = rt;
993 zd2 = zd;
994 }
995
996 if (fabs(zd2-zd1) < tol) break;
997 }
998 }
999 }
1000
1001 if (r == 0.0) {
1002 *phi = 0.0;
1003 } else {
1004 *phi = wcs_atan2d(x, -y);
1005 }
1006 *theta = 90.0 - zd*R2D;
1007
1008 return 0;
1009 }
1010
1011 /*============================================================================
1012 * ZEA: zenithal/azimuthal equal area projection.
1013 *
1014 * Given and/or returned:
1015 * prj->r0 r0; reset to 180/pi if 0.
1016 * prj->w[0] 2*r0
1017 * prj->w[1] 1/(2*r0)
1018 *===========================================================================*/
1019
zeaset(prj)1020 int zeaset(prj)
1021
1022 struct prjprm *prj;
1023
1024 {
1025 if (prj->r0 == 0.0) {
1026 prj->r0 = R2D;
1027 prj->w[0] = 360.0/PI;
1028 prj->w[1] = PI/360.0;
1029 } else {
1030 prj->w[0] = 2.0*prj->r0;
1031 prj->w[1] = 1.0/prj->w[0];
1032 }
1033
1034 prj->flag = PRJSET;
1035
1036 return 0;
1037 }
1038
1039 /*--------------------------------------------------------------------------*/
1040
zeafwd(phi,theta,prj,x,y)1041 int zeafwd(phi, theta, prj, x, y)
1042
1043 const double phi, theta;
1044 struct prjprm *prj;
1045 double *x, *y;
1046
1047 {
1048 double r;
1049
1050 if (prj->flag != PRJSET) {
1051 if (zeaset(prj)) return 1;
1052 }
1053
1054 r = prj->w[0]*wcs_sind((90.0 - theta)/2.0);
1055 *x = r*wcs_sind(phi);
1056 *y = -r*wcs_cosd(phi);
1057
1058 return 0;
1059 }
1060
1061 /*--------------------------------------------------------------------------*/
1062
zearev(x,y,prj,phi,theta)1063 int zearev(x, y, prj, phi, theta)
1064
1065 const double x, y;
1066 struct prjprm *prj;
1067 double *phi, *theta;
1068
1069 {
1070 double r, s;
1071 const double tol = 1.0e-12;
1072
1073 if (prj->flag != PRJSET) {
1074 if (zeaset(prj)) return 1;
1075 }
1076
1077 r = sqrt(x*x + y*y);
1078 if (r == 0.0) {
1079 *phi = 0.0;
1080 } else {
1081 *phi = wcs_atan2d(x, -y);
1082 }
1083
1084 s = r*prj->w[1];
1085 if (fabs(s) > 1.0) {
1086 if (fabs(r - prj->w[0]) < tol) {
1087 *theta = -90.0;
1088 } else {
1089 return 2;
1090 }
1091 } else {
1092 *theta = 90.0 - 2.0*wcs_asind(s);
1093 }
1094
1095 return 0;
1096 }
1097
1098 /*============================================================================
1099 * AIR: Airy's projection.
1100 *
1101 * Given:
1102 * prj->p[1] Latitude theta_b within which the error is minimized,
1103 * in degrees.
1104 *
1105 * Given and/or returned:
1106 * prj->r0 r0; reset to 180/pi if 0.
1107 * prj->w[0] 2*r0
1108 * prj->w[1] ln(cos(xi_b))/tan(xi_b)**2, where xi_b = (90-theta_b)/2
1109 * prj->w[2] 1/2 - prj->w[1]
1110 * prj->w[3] 2*r0*prj->w[2]
1111 * prj->w[4] tol, cutoff for using small angle approximation, in
1112 * radians.
1113 * prj->w[5] prj->w[2]*tol
1114 * prj->w[6] (180/pi)/prj->w[2]
1115 *===========================================================================*/
1116
airset(prj)1117 int airset(prj)
1118
1119 struct prjprm *prj;
1120
1121 {
1122 const double tol = 1.0e-4;
1123 double cxi;
1124
1125 if (prj->r0 == 0.0) prj->r0 = R2D;
1126
1127 prj->w[0] = 2.0*prj->r0;
1128 if (prj->p[1] == 90.0) {
1129 prj->w[1] = -0.5;
1130 prj->w[2] = 1.0;
1131 } else if (prj->p[1] > -90.0) {
1132 cxi = wcs_cosd((90.0 - prj->p[1])/2.0);
1133 prj->w[1] = log(cxi)*(cxi*cxi)/(1.0-cxi*cxi);
1134 prj->w[2] = 0.5 - prj->w[1];
1135 } else {
1136 return 1;
1137 }
1138
1139 prj->w[3] = prj->w[0] * prj->w[2];
1140 prj->w[4] = tol;
1141 prj->w[5] = prj->w[2]*tol;
1142 prj->w[6] = R2D/prj->w[2];
1143
1144 prj->flag = PRJSET;
1145
1146 return 0;
1147 }
1148
1149 /*--------------------------------------------------------------------------*/
1150
airfwd(phi,theta,prj,x,y)1151 int airfwd(phi, theta, prj, x, y)
1152
1153 const double phi, theta;
1154 struct prjprm *prj;
1155 double *x, *y;
1156
1157 {
1158 double cxi, r, txi, xi;
1159
1160 if (prj->flag != PRJSET) {
1161 if (airset(prj)) return 1;
1162 }
1163
1164 if (theta == 90.0) {
1165 r = 0.0;
1166 } else if (theta > -90.0) {
1167 xi = D2R*(90.0 - theta)/2.0;
1168 if (xi < prj->w[4]) {
1169 r = xi*prj->w[3];
1170 } else {
1171 cxi = wcs_cosd((90.0 - theta)/2.0);
1172 txi = sqrt(1.0-cxi*cxi)/cxi;
1173 r = -prj->w[0]*(log(cxi)/txi + prj->w[1]*txi);
1174 }
1175 } else {
1176 return 2;
1177 }
1178
1179 *x = r*wcs_sind(phi);
1180 *y = -r*wcs_cosd(phi);
1181
1182 return 0;
1183 }
1184
1185 /*--------------------------------------------------------------------------*/
1186
airrev(x,y,prj,phi,theta)1187 int airrev(x, y, prj, phi, theta)
1188
1189 const double x, y;
1190 struct prjprm *prj;
1191 double *phi, *theta;
1192
1193 {
1194 int j;
1195 double cxi, lambda, r, r1, r2, rt, txi, x1, x2, xi;
1196 const double tol = 1.0e-12;
1197
1198 if (prj->flag != PRJSET) {
1199 if (airset(prj)) return 1;
1200 }
1201
1202 r = sqrt(x*x + y*y)/prj->w[0];
1203
1204 if (r == 0.0) {
1205 xi = 0.0;
1206 } else if (r < prj->w[5]) {
1207 xi = r*prj->w[6];
1208 } else {
1209 /* Find a solution interval. */
1210 x1 = 1.0;
1211 r1 = 0.0;
1212 for (j = 0; j < 30; j++) {
1213 x2 = x1/2.0;
1214 txi = sqrt(1.0-x2*x2)/x2;
1215 r2 = -(log(x2)/txi + prj->w[1]*txi);
1216
1217 if (r2 >= r) break;
1218 x1 = x2;
1219 r1 = r2;
1220 }
1221 if (j == 30) return 2;
1222
1223 for (j = 0; j < 100; j++) {
1224 /* Weighted division of the interval. */
1225 lambda = (r2-r)/(r2-r1);
1226 if (lambda < 0.1) {
1227 lambda = 0.1;
1228 } else if (lambda > 0.9) {
1229 lambda = 0.9;
1230 }
1231 cxi = x2 - lambda*(x2-x1);
1232
1233 txi = sqrt(1.0-cxi*cxi)/cxi;
1234 rt = -(log(cxi)/txi + prj->w[1]*txi);
1235
1236 if (rt < r) {
1237 if (r-rt < tol) break;
1238 r1 = rt;
1239 x1 = cxi;
1240 } else {
1241 if (rt-r < tol) break;
1242 r2 = rt;
1243 x2 = cxi;
1244 }
1245 }
1246 if (j == 100) return 2;
1247
1248 xi = wcs_acosd(cxi);
1249 }
1250
1251 if (r == 0.0) {
1252 *phi = 0.0;
1253 } else {
1254 *phi = wcs_atan2d(x, -y);
1255 }
1256 *theta = 90.0 - 2.0*xi;
1257
1258 return 0;
1259 }
1260
1261 /*============================================================================
1262 * CYP: cylindrical perspective projection.
1263 *
1264 * Given:
1265 * prj->p[1] Distance of point of projection from the centre of the
1266 * generating sphere, mu, in units of r0.
1267 * prj->p[2] Radius of the cylinder of projection, lambda, in units
1268 * of r0.
1269 *
1270 * Given and/or returned:
1271 * prj->r0 r0; reset to 180/pi if 0.
1272 * prj->w[0] r0*lambda*(pi/180)
1273 * prj->w[1] (180/pi)/(r0*lambda)
1274 * prj->w[2] r0*(mu + lambda)
1275 * prj->w[3] 1/(r0*(mu + lambda))
1276 *===========================================================================*/
1277
cypset(prj)1278 int cypset(prj)
1279
1280 struct prjprm *prj;
1281
1282 {
1283 if (prj->r0 == 0.0) {
1284 prj->r0 = R2D;
1285
1286 prj->w[0] = prj->p[2];
1287 if (prj->w[0] == 0.0) {
1288 return 1;
1289 }
1290
1291 prj->w[1] = 1.0/prj->w[0];
1292
1293 prj->w[2] = R2D*(prj->p[1] + prj->p[2]);
1294 if (prj->w[2] == 0.0) {
1295 return 1;
1296 }
1297
1298 prj->w[3] = 1.0/prj->w[2];
1299 } else {
1300 prj->w[0] = prj->r0*prj->p[2]*D2R;
1301 if (prj->w[0] == 0.0) {
1302 return 1;
1303 }
1304
1305 prj->w[1] = 1.0/prj->w[0];
1306
1307 prj->w[2] = prj->r0*(prj->p[1] + prj->p[2]);
1308 if (prj->w[2] == 0.0) {
1309 return 1;
1310 }
1311
1312 prj->w[3] = 1.0/prj->w[2];
1313 }
1314
1315 prj->flag = PRJSET;
1316
1317 return 0;
1318 }
1319
1320 /*--------------------------------------------------------------------------*/
1321
cypfwd(phi,theta,prj,x,y)1322 int cypfwd(phi, theta, prj, x, y)
1323
1324 const double phi, theta;
1325 struct prjprm *prj;
1326 double *x, *y;
1327
1328 {
1329 double s;
1330
1331 if (prj->flag != PRJSET) {
1332 if (cypset(prj)) return 1;
1333 }
1334
1335 s = prj->p[1] + wcs_cosd(theta);
1336 if (s == 0.0) {
1337 return 2;
1338 }
1339
1340 *x = prj->w[0]*phi;
1341 *y = prj->w[2]*wcs_sind(theta)/s;
1342
1343 return 0;
1344 }
1345
1346 /*--------------------------------------------------------------------------*/
1347
cyprev(x,y,prj,phi,theta)1348 int cyprev(x, y, prj, phi, theta)
1349
1350 const double x, y;
1351 struct prjprm *prj;
1352 double *phi, *theta;
1353
1354 {
1355 double eta;
1356
1357 if (prj->flag != PRJSET) {
1358 if (cypset(prj)) return 1;
1359 }
1360
1361 *phi = x*prj->w[1];
1362 eta = y*prj->w[3];
1363 *theta = wcs_atan2d(eta,1.0) + wcs_asind(eta*prj->p[1]/sqrt(eta*eta+1.0));
1364
1365 return 0;
1366 }
1367
1368 /*============================================================================
1369 * CAR: Cartesian projection.
1370 *
1371 * Given and/or returned:
1372 * prj->r0 r0; reset to 180/pi if 0.
1373 * prj->w[0] r0*(pi/180)
1374 * prj->w[1] (180/pi)/r0
1375 *===========================================================================*/
1376
carset(prj)1377 int carset(prj)
1378
1379 struct prjprm *prj;
1380
1381 {
1382
1383 if (prj->r0 == 0.0) {
1384 prj->r0 = R2D;
1385 prj->w[0] = 1.0;
1386 prj->w[1] = 1.0;
1387 } else {
1388 prj->w[0] = prj->r0*D2R;
1389 prj->w[1] = 1.0/prj->w[0];
1390 }
1391
1392 prj->flag = PRJSET;
1393
1394 return 0;
1395 }
1396
1397 /*--------------------------------------------------------------------------*/
1398
carfwd(phi,theta,prj,x,y)1399 int carfwd(phi, theta, prj, x, y)
1400
1401 const double phi, theta;
1402 struct prjprm *prj;
1403 double *x, *y;
1404
1405 {
1406 if (prj->flag != PRJSET) {
1407 if (carset(prj)) return 1;
1408 }
1409
1410 *x = prj->w[0]*phi;
1411 *y = prj->w[0]*theta;
1412
1413 return 0;
1414 }
1415
1416 /*--------------------------------------------------------------------------*/
1417
carrev(x,y,prj,phi,theta)1418 int carrev(x, y, prj, phi, theta)
1419
1420 const double x, y;
1421 struct prjprm *prj;
1422 double *phi, *theta;
1423
1424 {
1425 if (prj->flag != PRJSET) {
1426 if (carset(prj)) return 1;
1427 }
1428
1429 *phi = prj->w[1]*x;
1430 *theta = prj->w[1]*y;
1431
1432 return 0;
1433 }
1434
1435 /*============================================================================
1436 * MER: Mercator's projection.
1437 *
1438 * Given and/or returned:
1439 * prj->r0 r0; reset to 180/pi if 0.
1440 * prj->w[0] r0*(pi/180)
1441 * prj->w[1] (180/pi)/r0
1442 *===========================================================================*/
1443
merset(prj)1444 int merset(prj)
1445
1446 struct prjprm *prj;
1447
1448 {
1449 if (prj->r0 == 0.0) {
1450 prj->r0 = R2D;
1451 prj->w[0] = 1.0;
1452 prj->w[1] = 1.0;
1453 } else {
1454 prj->w[0] = prj->r0*D2R;
1455 prj->w[1] = 1.0/prj->w[0];
1456 }
1457
1458 prj->flag = PRJSET;
1459
1460 return 0;
1461 }
1462
1463 /*--------------------------------------------------------------------------*/
1464
merfwd(phi,theta,prj,x,y)1465 int merfwd(phi, theta, prj, x, y)
1466
1467 const double phi, theta;
1468 struct prjprm *prj;
1469 double *x, *y;
1470
1471 {
1472 if (prj->flag != PRJSET) {
1473 if (merset(prj)) return 1;
1474 }
1475
1476 if (theta <= -90.0 || theta >= 90.0) {
1477 return 2;
1478 }
1479
1480 *x = prj->w[0]*phi;
1481 *y = prj->r0*log(wcs_tand((90.0+theta)/2.0));
1482
1483 return 0;
1484 }
1485
1486 /*--------------------------------------------------------------------------*/
1487
merrev(x,y,prj,phi,theta)1488 int merrev(x, y, prj, phi, theta)
1489
1490 const double x, y;
1491 struct prjprm *prj;
1492 double *phi, *theta;
1493
1494 {
1495 if (prj->flag != PRJSET) {
1496 if (merset(prj)) return 1;
1497 }
1498
1499 *phi = x*prj->w[1];
1500 *theta = 2.0*wcs_atand(exp(y/prj->r0)) - 90.0;
1501
1502 return 0;
1503 }
1504
1505 /*============================================================================
1506 * CEA: cylindrical equal area projection.
1507 *
1508 * Given:
1509 * prj->p[1] Square of the cosine of the latitude at which the
1510 * projection is conformal, lambda.
1511 *
1512 * Given and/or returned:
1513 * prj->r0 r0; reset to 180/pi if 0.
1514 * prj->w[0] r0*(pi/180)
1515 * prj->w[1] (180/pi)/r0
1516 * prj->w[2] r0/lambda
1517 * prj->w[3] lambda/r0
1518 *===========================================================================*/
1519
ceaset(prj)1520 int ceaset(prj)
1521
1522 struct prjprm *prj;
1523
1524 {
1525 if (prj->r0 == 0.0) {
1526 prj->r0 = R2D;
1527 prj->w[0] = 1.0;
1528 prj->w[1] = 1.0;
1529 if (prj->p[1] <= 0.0 || prj->p[1] > 1.0) {
1530 return 1;
1531 }
1532 prj->w[2] = prj->r0/prj->p[1];
1533 prj->w[3] = prj->p[1]/prj->r0;
1534 } else {
1535 prj->w[0] = prj->r0*D2R;
1536 prj->w[1] = R2D/prj->r0;
1537 if (prj->p[1] <= 0.0 || prj->p[1] > 1.0) {
1538 return 1;
1539 }
1540 prj->w[2] = prj->r0/prj->p[1];
1541 prj->w[3] = prj->p[1]/prj->r0;
1542 }
1543
1544 prj->flag = PRJSET;
1545
1546 return 0;
1547 }
1548
1549 /*--------------------------------------------------------------------------*/
1550
ceafwd(phi,theta,prj,x,y)1551 int ceafwd(phi, theta, prj, x, y)
1552
1553 const double phi, theta;
1554 struct prjprm *prj;
1555 double *x, *y;
1556
1557 {
1558 if (prj->flag != PRJSET) {
1559 if (ceaset(prj)) return 1;
1560 }
1561
1562 *x = prj->w[0]*phi;
1563 *y = prj->w[2]*wcs_sind(theta);
1564
1565 return 0;
1566 }
1567
1568 /*--------------------------------------------------------------------------*/
1569
cearev(x,y,prj,phi,theta)1570 int cearev(x, y, prj, phi, theta)
1571
1572 const double x, y;
1573 struct prjprm *prj;
1574 double *phi, *theta;
1575
1576 {
1577 double s;
1578 const double tol = 1.0e-13;
1579
1580 if (prj->flag != PRJSET) {
1581 if (ceaset(prj)) return 1;
1582 }
1583
1584 s = y*prj->w[3];
1585 if (fabs(s) > 1.0) {
1586 if (fabs(s) > 1.0+tol) {
1587 return 2;
1588 }
1589 s = copysign(1.0,s);
1590 }
1591
1592 *phi = x*prj->w[1];
1593 *theta = wcs_asind(s);
1594
1595 return 0;
1596 }
1597
1598 /*============================================================================
1599 * COP: conic perspective projection.
1600 *
1601 * Given:
1602 * prj->p[1] sigma = (theta2+theta1)/2
1603 * prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the
1604 * latitudes of the standard parallels, in degrees.
1605 *
1606 * Given and/or returned:
1607 * prj->r0 r0; reset to 180/pi if 0.
1608 * prj->w[0] C = sin(sigma)
1609 * prj->w[1] 1/C
1610 * prj->w[2] Y0 = r0*cos(delta)*cot(sigma)
1611 * prj->w[3] r0*cos(delta)
1612 * prj->w[4] 1/(r0*cos(delta)
1613 * prj->w[5] cot(sigma)
1614 *===========================================================================*/
1615
copset(prj)1616 int copset(prj)
1617
1618 struct prjprm *prj;
1619
1620 {
1621 if (prj->r0 == 0.0) prj->r0 = R2D;
1622
1623 prj->w[0] = wcs_sind(prj->p[1]);
1624 if (prj->w[0] == 0.0) {
1625 return 1;
1626 }
1627
1628 prj->w[1] = 1.0/prj->w[0];
1629
1630 prj->w[3] = prj->r0*wcs_cosd(prj->p[2]);
1631 if (prj->w[3] == 0.0) {
1632 return 1;
1633 }
1634
1635 prj->w[4] = 1.0/prj->w[3];
1636 prj->w[5] = 1.0/wcs_tand(prj->p[1]);
1637
1638 prj->w[2] = prj->w[3]*prj->w[5];
1639
1640 if (prj->flag == -1) {
1641 prj->flag = -PRJSET;
1642 } else {
1643 prj->flag = PRJSET;
1644 }
1645
1646 return 0;
1647 }
1648
1649 /*--------------------------------------------------------------------------*/
1650
copfwd(phi,theta,prj,x,y)1651 int copfwd(phi, theta, prj, x, y)
1652
1653 const double phi, theta;
1654 struct prjprm *prj;
1655 double *x, *y;
1656
1657 {
1658 double a, r, s, t;
1659
1660 if (abs(prj->flag) != PRJSET) {
1661 if (copset(prj)) return 1;
1662 }
1663
1664 t = theta - prj->p[1];
1665 s = wcs_cosd(t);
1666 if (s == 0.0) {
1667 return 2;
1668 }
1669
1670 a = prj->w[0]*phi;
1671 r = prj->w[2] - prj->w[3]*wcs_sind(t)/s;
1672
1673 *x = r*wcs_sind(a);
1674 *y = prj->w[2] - r*wcs_cosd(a);
1675
1676 if (prj->flag == PRJSET && r*prj->w[0] < 0.0) {
1677 return 2;
1678 }
1679
1680 return 0;
1681 }
1682
1683 /*--------------------------------------------------------------------------*/
1684
coprev(x,y,prj,phi,theta)1685 int coprev(x, y, prj, phi, theta)
1686
1687 const double x, y;
1688 struct prjprm *prj;
1689 double *phi, *theta;
1690
1691 {
1692 double a, dy, r;
1693
1694 if (abs(prj->flag) != PRJSET) {
1695 if (copset(prj)) return 1;
1696 }
1697
1698 dy = prj->w[2] - y;
1699 r = sqrt(x*x + dy*dy);
1700 if (prj->p[1] < 0.0) r = -r;
1701
1702 if (r == 0.0) {
1703 a = 0.0;
1704 } else {
1705 a = wcs_atan2d(x/r, dy/r);
1706 }
1707
1708 *phi = a*prj->w[1];
1709 *theta = prj->p[1] + wcs_atand(prj->w[5] - r*prj->w[4]);
1710
1711 return 0;
1712 }
1713
1714 /*============================================================================
1715 * COD: conic equidistant projection.
1716 *
1717 * Given:
1718 * prj->p[1] sigma = (theta2+theta1)/2
1719 * prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the
1720 * latitudes of the standard parallels, in degrees.
1721 *
1722 * Given and/or returned:
1723 * prj->r0 r0; reset to 180/pi if 0.
1724 * prj->w[0] C = r0*sin(sigma)*sin(delta)/delta
1725 * prj->w[1] 1/C
1726 * prj->w[2] Y0 = delta*cot(delta)*cot(sigma)
1727 * prj->w[3] Y0 + sigma
1728 *===========================================================================*/
1729
codset(prj)1730 int codset(prj)
1731
1732 struct prjprm *prj;
1733
1734 {
1735 if (prj->r0 == 0.0) prj->r0 = R2D;
1736
1737 if (prj->p[2] == 0.0) {
1738 prj->w[0] = prj->r0*wcs_sind(prj->p[1])*D2R;
1739 } else {
1740 prj->w[0] = prj->r0*wcs_sind(prj->p[1])*wcs_sind(prj->p[2])/prj->p[2];
1741 }
1742
1743 if (prj->w[0] == 0.0) {
1744 return 1;
1745 }
1746
1747 prj->w[1] = 1.0/prj->w[0];
1748 prj->w[2] = prj->r0*wcs_cosd(prj->p[2])*wcs_cosd(prj->p[1])/prj->w[0];
1749 prj->w[3] = prj->w[2] + prj->p[1];
1750
1751 prj->flag = PRJSET;
1752
1753 return 0;
1754 }
1755
1756 /*--------------------------------------------------------------------------*/
1757
codfwd(phi,theta,prj,x,y)1758 int codfwd(phi, theta, prj, x, y)
1759
1760 const double phi, theta;
1761 struct prjprm *prj;
1762 double *x, *y;
1763
1764 {
1765 double a, r;
1766
1767 if (prj->flag != PRJSET) {
1768 if (codset(prj)) return 1;
1769 }
1770
1771 a = prj->w[0]*phi;
1772 r = prj->w[3] - theta;
1773
1774 *x = r*wcs_sind(a);
1775 *y = prj->w[2] - r*wcs_cosd(a);
1776
1777 return 0;
1778 }
1779
1780 /*--------------------------------------------------------------------------*/
1781
codrev(x,y,prj,phi,theta)1782 int codrev(x, y, prj, phi, theta)
1783
1784 const double x, y;
1785 struct prjprm *prj;
1786 double *phi, *theta;
1787
1788 {
1789 double a, dy, r;
1790
1791 if (prj->flag != PRJSET) {
1792 if (codset(prj)) return 1;
1793 }
1794
1795 dy = prj->w[2] - y;
1796 r = sqrt(x*x + dy*dy);
1797 if (prj->p[1] < 0.0) r = -r;
1798
1799 if (r == 0.0) {
1800 a = 0.0;
1801 } else {
1802 a = wcs_atan2d(x/r, dy/r);
1803 }
1804
1805 *phi = a*prj->w[1];
1806 *theta = prj->w[3] - r;
1807
1808 return 0;
1809 }
1810
1811 /*============================================================================
1812 * COE: conic equal area projection.
1813 *
1814 * Given:
1815 * prj->p[1] sigma = (theta2+theta1)/2
1816 * prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the
1817 * latitudes of the standard parallels, in degrees.
1818 *
1819 * Given and/or returned:
1820 * prj->r0 r0; reset to 180/pi if 0.
1821 * prj->w[0] C = (sin(theta1) + sin(theta2))/2
1822 * prj->w[1] 1/C
1823 * prj->w[2] Y0 = chi*sqrt(psi - 2C*wcs_sind(sigma))
1824 * prj->w[3] chi = r0/C
1825 * prj->w[4] psi = 1 + sin(theta1)*sin(theta2)
1826 * prj->w[5] 2C
1827 * prj->w[6] (1 + sin(theta1)*sin(theta2))*(r0/C)**2
1828 * prj->w[7] C/(2*r0**2)
1829 * prj->w[8] chi*sqrt(psi + 2C)
1830 *===========================================================================*/
1831
coeset(prj)1832 int coeset(prj)
1833
1834 struct prjprm *prj;
1835
1836 {
1837 double theta1, theta2;
1838
1839 if (prj->r0 == 0.0) prj->r0 = R2D;
1840
1841 theta1 = prj->p[1] - prj->p[2];
1842 theta2 = prj->p[1] + prj->p[2];
1843
1844 prj->w[0] = (wcs_sind(theta1) + wcs_sind(theta2))/2.0;
1845 if (prj->w[0] == 0.0) {
1846 return 1;
1847 }
1848
1849 prj->w[1] = 1.0/prj->w[0];
1850
1851 prj->w[3] = prj->r0/prj->w[0];
1852 prj->w[4] = 1.0 + wcs_sind(theta1)*wcs_sind(theta2);
1853 prj->w[5] = 2.0*prj->w[0];
1854 prj->w[6] = prj->w[3]*prj->w[3]*prj->w[4];
1855 prj->w[7] = 1.0/(2.0*prj->r0*prj->w[3]);
1856 prj->w[8] = prj->w[3]*sqrt(prj->w[4] + prj->w[5]);
1857
1858 prj->w[2] = prj->w[3]*sqrt(prj->w[4] - prj->w[5]*wcs_sind(prj->p[1]));
1859
1860 prj->flag = PRJSET;
1861
1862 return 0;
1863 }
1864
1865 /*--------------------------------------------------------------------------*/
1866
coefwd(phi,theta,prj,x,y)1867 int coefwd(phi, theta, prj, x, y)
1868
1869 const double phi, theta;
1870 struct prjprm *prj;
1871 double *x, *y;
1872
1873 {
1874 double a, r;
1875
1876 if (prj->flag != PRJSET) {
1877 if (coeset(prj)) return 1;
1878 }
1879
1880 a = phi*prj->w[0];
1881 if (theta == -90.0) {
1882 r = prj->w[8];
1883 } else {
1884 r = prj->w[3]*sqrt(prj->w[4] - prj->w[5]*wcs_sind(theta));
1885 }
1886
1887 *x = r*wcs_sind(a);
1888 *y = prj->w[2] - r*wcs_cosd(a);
1889
1890 return 0;
1891 }
1892
1893 /*--------------------------------------------------------------------------*/
1894
coerev(x,y,prj,phi,theta)1895 int coerev(x, y, prj, phi, theta)
1896
1897 const double x, y;
1898 struct prjprm *prj;
1899 double *phi, *theta;
1900
1901 {
1902 double a, dy, r, w;
1903 const double tol = 1.0e-12;
1904
1905 if (prj->flag != PRJSET) {
1906 if (coeset(prj)) return 1;
1907 }
1908
1909 dy = prj->w[2] - y;
1910 r = sqrt(x*x + dy*dy);
1911 if (prj->p[1] < 0.0) r = -r;
1912
1913 if (r == 0.0) {
1914 a = 0.0;
1915 } else {
1916 a = wcs_atan2d(x/r, dy/r);
1917 }
1918
1919 *phi = a*prj->w[1];
1920 if (fabs(r - prj->w[8]) < tol) {
1921 *theta = -90.0;
1922 } else {
1923 w = (prj->w[6] - r*r)*prj->w[7];
1924 if (fabs(w) > 1.0) {
1925 if (fabs(w-1.0) < tol) {
1926 *theta = 90.0;
1927 } else if (fabs(w+1.0) < tol) {
1928 *theta = -90.0;
1929 } else {
1930 return 2;
1931 }
1932 } else {
1933 *theta = wcs_asind(w);
1934 }
1935 }
1936
1937 return 0;
1938 }
1939
1940 /*============================================================================
1941 * COO: conic orthomorphic projection.
1942 *
1943 * Given:
1944 * prj->p[1] sigma = (theta2+theta1)/2
1945 * prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the
1946 * latitudes of the standard parallels, in degrees.
1947 *
1948 * Given and/or returned:
1949 * prj->r0 r0; reset to 180/pi if 0.
1950 * prj->w[0] C = ln(cos(theta2)/cos(theta1))/ln(tan(tau2)/tan(tau1))
1951 * where tau1 = (90 - theta1)/2
1952 * tau2 = (90 - theta2)/2
1953 * prj->w[1] 1/C
1954 * prj->w[2] Y0 = psi*tan((90-sigma)/2)**C
1955 * prj->w[3] psi = (r0*cos(theta1)/C)/tan(tau1)**C
1956 * prj->w[4] 1/psi
1957 *===========================================================================*/
1958
cooset(prj)1959 int cooset(prj)
1960
1961 struct prjprm *prj;
1962
1963 {
1964 double cos1, cos2, tan1, tan2, theta1, theta2;
1965
1966 if (prj->r0 == 0.0) prj->r0 = R2D;
1967
1968 theta1 = prj->p[1] - prj->p[2];
1969 theta2 = prj->p[1] + prj->p[2];
1970
1971 tan1 = wcs_tand((90.0 - theta1)/2.0);
1972 cos1 = wcs_cosd(theta1);
1973
1974 if (theta1 == theta2) {
1975 prj->w[0] = wcs_sind(theta1);
1976 } else {
1977 tan2 = wcs_tand((90.0 - theta2)/2.0);
1978 cos2 = wcs_cosd(theta2);
1979 prj->w[0] = log(cos2/cos1)/log(tan2/tan1);
1980 }
1981 if (prj->w[0] == 0.0) {
1982 return 1;
1983 }
1984
1985 prj->w[1] = 1.0/prj->w[0];
1986
1987 prj->w[3] = prj->r0*(cos1/prj->w[0])/pow(tan1,prj->w[0]);
1988 if (prj->w[3] == 0.0) {
1989 return 1;
1990 }
1991 prj->w[2] = prj->w[3]*pow(wcs_tand((90.0 - prj->p[1])/2.0),prj->w[0]);
1992 prj->w[4] = 1.0/prj->w[3];
1993
1994 prj->flag = PRJSET;
1995
1996 return 0;
1997 }
1998
1999 /*--------------------------------------------------------------------------*/
2000
coofwd(phi,theta,prj,x,y)2001 int coofwd(phi, theta, prj, x, y)
2002
2003 const double phi, theta;
2004 struct prjprm *prj;
2005 double *x, *y;
2006
2007 {
2008 double a, r;
2009
2010 if (prj->flag != PRJSET) {
2011 if (cooset(prj)) return 1;
2012 }
2013
2014 a = prj->w[0]*phi;
2015 if (theta == -90.0) {
2016 if (prj->w[0] < 0.0) {
2017 r = 0.0;
2018 } else {
2019 return 2;
2020 }
2021 } else {
2022 r = prj->w[3]*pow(wcs_tand((90.0 - theta)/2.0),prj->w[0]);
2023 }
2024
2025 *x = r*wcs_sind(a);
2026 *y = prj->w[2] - r*wcs_cosd(a);
2027
2028 return 0;
2029 }
2030
2031 /*--------------------------------------------------------------------------*/
2032
coorev(x,y,prj,phi,theta)2033 int coorev(x, y, prj, phi, theta)
2034
2035 const double x, y;
2036 struct prjprm *prj;
2037 double *phi, *theta;
2038
2039 {
2040 double a, dy, r;
2041
2042 if (prj->flag != PRJSET) {
2043 if (cooset(prj)) return 1;
2044 }
2045
2046 dy = prj->w[2] - y;
2047 r = sqrt(x*x + dy*dy);
2048 if (prj->p[1] < 0.0) r = -r;
2049
2050 if (r == 0.0) {
2051 a = 0.0;
2052 } else {
2053 a = wcs_atan2d(x/r, dy/r);
2054 }
2055
2056 *phi = a*prj->w[1];
2057 if (r == 0.0) {
2058 if (prj->w[0] < 0.0) {
2059 *theta = -90.0;
2060 } else {
2061 return 2;
2062 }
2063 } else {
2064 *theta = 90.0 - 2.0*wcs_atand(pow(r*prj->w[4],prj->w[1]));
2065 }
2066
2067 return 0;
2068 }
2069
2070 /*============================================================================
2071 * BON: Bonne's projection.
2072 *
2073 * Given:
2074 * prj->p[1] Bonne conformal latitude, theta1, in degrees.
2075 *
2076 * Given and/or returned:
2077 * prj->r0 r0; reset to 180/pi if 0.
2078 * prj->w[1] r0*pi/180
2079 * prj->w[2] Y0 = r0*cot(theta1) + theta1*pi/180)
2080 *===========================================================================*/
2081
bonset(prj)2082 int bonset(prj)
2083
2084 struct prjprm *prj;
2085
2086 {
2087 if (prj->r0 == 0.0) {
2088 prj->r0 = R2D;
2089 prj->w[1] = 1.0;
2090 prj->w[2] = prj->r0*wcs_cosd(prj->p[1])/wcs_sind(prj->p[1]) + prj->p[1];
2091 } else {
2092 prj->w[1] = prj->r0*D2R;
2093 prj->w[2] = prj->r0*(wcs_cosd(prj->p[1])/wcs_sind(prj->p[1]) + prj->p[1]*D2R);
2094 }
2095
2096 prj->flag = PRJSET;
2097
2098 return 0;
2099 }
2100
2101 /*--------------------------------------------------------------------------*/
2102
bonfwd(phi,theta,prj,x,y)2103 int bonfwd(phi, theta, prj, x, y)
2104
2105 const double phi, theta;
2106 struct prjprm *prj;
2107 double *x, *y;
2108
2109 {
2110 double a, r;
2111
2112 if (prj->p[1] == 0.0) {
2113 /* Sanson-Flamsteed. */
2114 return glsfwd(phi, theta, prj, x, y);
2115 }
2116
2117 if (prj->flag != PRJSET) {
2118 if (bonset(prj)) return 1;
2119 }
2120
2121 r = prj->w[2] - theta*prj->w[1];
2122 a = prj->r0*phi*wcs_cosd(theta)/r;
2123
2124 *x = r*wcs_sind(a);
2125 *y = prj->w[2] - r*wcs_cosd(a);
2126
2127 return 0;
2128 }
2129
2130 /*--------------------------------------------------------------------------*/
2131
bonrev(x,y,prj,phi,theta)2132 int bonrev(x, y, prj, phi, theta)
2133
2134 const double x, y;
2135 struct prjprm *prj;
2136 double *phi, *theta;
2137
2138 {
2139 double a, dy, costhe, r;
2140
2141 if (prj->p[1] == 0.0) {
2142 /* Sanson-Flamsteed. */
2143 return glsrev(x, y, prj, phi, theta);
2144 }
2145
2146 if (prj->flag != PRJSET) {
2147 if (bonset(prj)) return 1;
2148 }
2149
2150 dy = prj->w[2] - y;
2151 r = sqrt(x*x + dy*dy);
2152 if (prj->p[1] < 0.0) r = -r;
2153
2154 if (r == 0.0) {
2155 a = 0.0;
2156 } else {
2157 a = wcs_atan2d(x/r, dy/r);
2158 }
2159
2160 *theta = (prj->w[2] - r)/prj->w[1];
2161 costhe = wcs_cosd(*theta);
2162 if (costhe == 0.0) {
2163 *phi = 0.0;
2164 } else {
2165 *phi = a*(r/prj->r0)/wcs_cosd(*theta);
2166 }
2167
2168 return 0;
2169 }
2170
2171 /*============================================================================
2172 * PCO: polyconic projection.
2173 *
2174 * Given and/or returned:
2175 * prj->r0 r0; reset to 180/pi if 0.
2176 * prj->w[0] r0*(pi/180)
2177 * prj->w[1] 1/r0
2178 * prj->w[2] 2*r0
2179 *===========================================================================*/
2180
pcoset(prj)2181 int pcoset(prj)
2182
2183 struct prjprm *prj;
2184
2185 {
2186 if (prj->r0 == 0.0) {
2187 prj->r0 = R2D;
2188 prj->w[0] = 1.0;
2189 prj->w[1] = 1.0;
2190 prj->w[2] = 360.0/PI;
2191 } else {
2192 prj->w[0] = prj->r0*D2R;
2193 prj->w[1] = 1.0/prj->w[0];
2194 prj->w[2] = 2.0*prj->r0;
2195 }
2196
2197 prj->flag = PRJSET;
2198
2199 return 0;
2200 }
2201
2202 /*--------------------------------------------------------------------------*/
2203
pcofwd(phi,theta,prj,x,y)2204 int pcofwd(phi, theta, prj, x, y)
2205
2206 const double phi, theta;
2207 struct prjprm *prj;
2208 double *x, *y;
2209
2210 {
2211 double a, costhe, cotthe, sinthe;
2212
2213 if (prj->flag != PRJSET) {
2214 if (pcoset(prj)) return 1;
2215 }
2216
2217 costhe = wcs_cosd(theta);
2218 sinthe = wcs_sind(theta);
2219 a = phi*sinthe;
2220
2221 if (sinthe == 0.0) {
2222 *x = prj->w[0]*phi;
2223 *y = 0.0;
2224 } else {
2225 cotthe = costhe/sinthe;
2226 *x = prj->r0*cotthe*wcs_sind(a);
2227 *y = prj->r0*(cotthe*(1.0 - wcs_cosd(a)) + theta*D2R);
2228 }
2229
2230 return 0;
2231 }
2232
2233 /*--------------------------------------------------------------------------*/
2234
pcorev(x,y,prj,phi,theta)2235 int pcorev(x, y, prj, phi, theta)
2236
2237 const double x, y;
2238 struct prjprm *prj;
2239 double *phi, *theta;
2240
2241 {
2242 int j;
2243 double f, fneg, fpos, lambda, tanthe, theneg, thepos, w, xp, xx, ymthe, yp;
2244 const double tol = 1.0e-12;
2245
2246 if (prj->flag != PRJSET) {
2247 if (pcoset(prj)) return 1;
2248 }
2249
2250 w = fabs(y*prj->w[1]);
2251 if (w < tol) {
2252 *phi = x*prj->w[1];
2253 *theta = 0.0;
2254 } else if (fabs(w-90.0) < tol) {
2255 *phi = 0.0;
2256 *theta = wcs_copysign(90.0,y);
2257 } else {
2258 /* Iterative solution using weighted division of the interval. */
2259 if (y > 0.0) {
2260 thepos = 90.0;
2261 } else {
2262 thepos = -90.0;
2263 }
2264 theneg = 0.0;
2265
2266 xx = x*x;
2267 ymthe = y - prj->w[0]*thepos;
2268 fpos = xx + ymthe*ymthe;
2269 fneg = -999.0;
2270
2271 for (j = 0; j < 64; j++) {
2272 if (fneg < -100.0) {
2273 /* Equal division of the interval. */
2274 *theta = (thepos+theneg)/2.0;
2275 } else {
2276 /* Weighted division of the interval. */
2277 lambda = fpos/(fpos-fneg);
2278 if (lambda < 0.1) {
2279 lambda = 0.1;
2280 } else if (lambda > 0.9) {
2281 lambda = 0.9;
2282 }
2283 *theta = thepos - lambda*(thepos-theneg);
2284 }
2285
2286 /* Compute the residue. */
2287 ymthe = y - prj->w[0]*(*theta);
2288 tanthe = wcs_tand(*theta);
2289 f = xx + ymthe*(ymthe - prj->w[2]/tanthe);
2290
2291 /* Check for convergence. */
2292 if (fabs(f) < tol) break;
2293 if (fabs(thepos-theneg) < tol) break;
2294
2295 /* Redefine the interval. */
2296 if (f > 0.0) {
2297 thepos = *theta;
2298 fpos = f;
2299 } else {
2300 theneg = *theta;
2301 fneg = f;
2302 }
2303 }
2304
2305 xp = prj->r0 - ymthe*tanthe;
2306 yp = x*tanthe;
2307 if (xp == 0.0 && yp == 0.0) {
2308 *phi = 0.0;
2309 } else {
2310 *phi = wcs_atan2d(yp, xp)/wcs_sind(*theta);
2311 }
2312 }
2313
2314 return 0;
2315 }
2316
2317 /*============================================================================
2318 * GLS: Sanson-Flamsteed ("global sinusoid") projection.
2319 *
2320 * Given and/or returned:
2321 * prj->r0 r0; reset to 180/pi if 0.
2322 * prj->w[0] r0*(pi/180)
2323 * prj->w[1] (180/pi)/r0
2324 *===========================================================================*/
2325
glsset(prj)2326 int glsset(prj)
2327
2328 struct prjprm *prj;
2329
2330 {
2331 if (prj->r0 == 0.0) {
2332 prj->r0 = R2D;
2333 prj->w[0] = 1.0;
2334 prj->w[1] = 1.0;
2335 } else {
2336 prj->w[0] = prj->r0*D2R;
2337 prj->w[1] = 1.0/prj->w[0];
2338 }
2339
2340 prj->flag = PRJSET;
2341
2342 return 0;
2343 }
2344
2345 /*--------------------------------------------------------------------------*/
2346
glsfwd(phi,theta,prj,x,y)2347 int glsfwd(phi, theta, prj, x, y)
2348
2349 const double phi, theta;
2350 struct prjprm *prj;
2351 double *x, *y;
2352
2353 {
2354 if (prj->flag != PRJSET) {
2355 if (glsset(prj)) return 1;
2356 }
2357
2358 *x = prj->w[0]*phi*wcs_cosd(theta);
2359 *y = prj->w[0]*theta;
2360
2361 return 0;
2362 }
2363
2364 /*--------------------------------------------------------------------------*/
2365
glsrev(x,y,prj,phi,theta)2366 int glsrev(x, y, prj, phi, theta)
2367
2368 const double x, y;
2369 struct prjprm *prj;
2370 double *phi, *theta;
2371
2372 {
2373 double w;
2374
2375 if (prj->flag != PRJSET) {
2376 if (glsset(prj)) return 1;
2377 }
2378
2379 w = cos(y/prj->r0);
2380 if (w == 0.0) {
2381 *phi = 0.0;
2382 } else {
2383 *phi = x*prj->w[1]/cos(y/prj->r0);
2384 }
2385 *theta = y*prj->w[1];
2386
2387 return 0;
2388 }
2389
2390 /*============================================================================
2391 * PAR: parabolic projection.
2392 *
2393 * Given and/or returned:
2394 * prj->r0 r0; reset to 180/pi if 0.
2395 * prj->w[0] r0*(pi/180)
2396 * prj->w[1] (180/pi)/r0
2397 * prj->w[2] pi*r0
2398 * prj->w[3] 1/(pi*r0)
2399 *===========================================================================*/
2400
parset(prj)2401 int parset(prj)
2402
2403 struct prjprm *prj;
2404
2405 {
2406 if (prj->r0 == 0.0) {
2407 prj->r0 = R2D;
2408 prj->w[0] = 1.0;
2409 prj->w[1] = 1.0;
2410 prj->w[2] = 180.0;
2411 prj->w[3] = 1.0/prj->w[2];
2412 } else {
2413 prj->w[0] = prj->r0*D2R;
2414 prj->w[1] = 1.0/prj->w[0];
2415 prj->w[2] = PI*prj->r0;
2416 prj->w[3] = 1.0/prj->w[2];
2417 }
2418
2419 prj->flag = PRJSET;
2420
2421 return 0;
2422 }
2423
2424 /*--------------------------------------------------------------------------*/
2425
parfwd(phi,theta,prj,x,y)2426 int parfwd(phi, theta, prj, x, y)
2427
2428 const double phi, theta;
2429 struct prjprm *prj;
2430 double *x, *y;
2431
2432 {
2433 double s;
2434
2435 if (prj->flag != PRJSET) {
2436 if (parset(prj)) return 1;
2437 }
2438
2439 s = wcs_sind(theta/3.0);
2440 *x = prj->w[0]*phi*(1.0 - 4.0*s*s);
2441 *y = prj->w[2]*s;
2442
2443 return 0;
2444 }
2445
2446 /*--------------------------------------------------------------------------*/
2447
parrev(x,y,prj,phi,theta)2448 int parrev(x, y, prj, phi, theta)
2449
2450 const double x, y;
2451 struct prjprm *prj;
2452 double *phi, *theta;
2453
2454 {
2455 double s, t;
2456
2457 if (prj->flag != PRJSET) {
2458 if (parset(prj)) return 1;
2459 }
2460
2461 s = y*prj->w[3];
2462 if (s > 1.0 || s < -1.0) {
2463 return 2;
2464 }
2465
2466 t = 1.0 - 4.0*s*s;
2467 if (t == 0.0) {
2468 if (x == 0.0) {
2469 *phi = 0.0;
2470 } else {
2471 return 2;
2472 }
2473 } else {
2474 *phi = prj->w[1]*x/t;
2475 }
2476
2477 *theta = 3.0*wcs_asind(s);
2478
2479 return 0;
2480 }
2481
2482 /*============================================================================
2483 * AIT: Hammer-Aitoff projection.
2484 *
2485 * Given and/or returned:
2486 * prj->r0 r0; reset to 180/pi if 0.
2487 * prj->w[0] 2*r0**2
2488 * prj->w[1] 1/(2*r0)**2
2489 * prj->w[2] 1/(4*r0)**2
2490 * prj->w[3] 1/(2*r0)
2491 *===========================================================================*/
2492
aitset(prj)2493 int aitset(prj)
2494
2495 struct prjprm *prj;
2496
2497 {
2498 if (prj->r0 == 0.0) prj->r0 = R2D;
2499
2500 prj->w[0] = 2.0*prj->r0*prj->r0;
2501 prj->w[1] = 1.0/(2.0*prj->w[0]);
2502 prj->w[2] = prj->w[1]/4.0;
2503 prj->w[3] = 1.0/(2.0*prj->r0);
2504
2505 prj->flag = PRJSET;
2506
2507 return 0;
2508 }
2509
2510 /*--------------------------------------------------------------------------*/
2511
aitfwd(phi,theta,prj,x,y)2512 int aitfwd(phi, theta, prj, x, y)
2513
2514 const double phi, theta;
2515 struct prjprm *prj;
2516 double *x, *y;
2517
2518 {
2519 double costhe, w;
2520
2521 if (prj->flag != PRJSET) {
2522 if (aitset(prj)) return 1;
2523 }
2524
2525 costhe = wcs_cosd(theta);
2526 w = sqrt(prj->w[0]/(1.0 + costhe*wcs_cosd(phi/2.0)));
2527 *x = 2.0*w*costhe*wcs_sind(phi/2.0);
2528 *y = w*wcs_sind(theta);
2529
2530 return 0;
2531 }
2532
2533 /*--------------------------------------------------------------------------*/
2534
aitrev(x,y,prj,phi,theta)2535 int aitrev(x, y, prj, phi, theta)
2536
2537 const double x, y;
2538 struct prjprm *prj;
2539 double *phi, *theta;
2540
2541 {
2542 double s, u, xp, yp, z;
2543 const double tol = 1.0e-13;
2544
2545 if (prj->flag != PRJSET) {
2546 if (aitset(prj)) return 1;
2547 }
2548
2549 u = 1.0 - x*x*prj->w[2] - y*y*prj->w[1];
2550 if (u < 0.0) {
2551 if (u < -tol) {
2552 return 2;
2553 }
2554
2555 u = 0.0;
2556 }
2557
2558 z = sqrt(u);
2559 s = z*y/prj->r0;
2560 if (fabs(s) > 1.0) {
2561 if (fabs(s) > 1.0+tol) {
2562 return 2;
2563 }
2564 s = wcs_copysign(1.0,s);
2565 }
2566
2567 xp = 2.0*z*z - 1.0;
2568 yp = z*x*prj->w[3];
2569 if (xp == 0.0 && yp == 0.0) {
2570 *phi = 0.0;
2571 } else {
2572 *phi = 2.0*wcs_atan2d(yp, xp);
2573 }
2574 *theta = wcs_asind(s);
2575
2576 return 0;
2577 }
2578
2579 /*============================================================================
2580 * MOL: Mollweide's projection.
2581 *
2582 * Given and/or returned:
2583 * prj->r0 r0; reset to 180/pi if 0.
2584 * prj->w[0] sqrt(2)*r0
2585 * prj->w[1] sqrt(2)*r0/90
2586 * prj->w[2] 1/(sqrt(2)*r0)
2587 * prj->w[3] 90/r0
2588 *===========================================================================*/
2589
molset(prj)2590 int molset(prj)
2591
2592 struct prjprm *prj;
2593
2594 {
2595 if (prj->r0 == 0.0) prj->r0 = R2D;
2596
2597 prj->w[0] = SQRT2*prj->r0;
2598 prj->w[1] = prj->w[0]/90.0;
2599 prj->w[2] = 1.0/prj->w[0];
2600 prj->w[3] = 90.0/prj->r0;
2601 prj->w[4] = 2.0/PI;
2602
2603 prj->flag = PRJSET;
2604
2605 return 0;
2606 }
2607
2608 /*--------------------------------------------------------------------------*/
2609
molfwd(phi,theta,prj,x,y)2610 int molfwd(phi, theta, prj, x, y)
2611
2612 const double phi, theta;
2613 struct prjprm *prj;
2614 double *x, *y;
2615
2616 {
2617 int j;
2618 double alpha, resid, u, v, v0, v1;
2619 const double tol = 1.0e-13;
2620
2621 if (prj->flag != PRJSET) {
2622 if (molset(prj)) return 1;
2623 }
2624
2625 if (fabs(theta) == 90.0) {
2626 *x = 0.0;
2627 *y = wcs_copysign(prj->w[0],theta);
2628 } else if (theta == 0.0) {
2629 *x = prj->w[1]*phi;
2630 *y = 0.0;
2631 } else {
2632 u = PI*wcs_sind(theta);
2633 v0 = -PI;
2634 v1 = PI;
2635 v = u;
2636 for (j = 0; j < 100; j++) {
2637 resid = (v - u) + sin(v);
2638 if (resid < 0.0) {
2639 if (resid > -tol) break;
2640 v0 = v;
2641 } else {
2642 if (resid < tol) break;
2643 v1 = v;
2644 }
2645 v = (v0 + v1)/2.0;
2646 }
2647
2648 alpha = v/2.0;
2649 *x = prj->w[1]*phi*cos(alpha);
2650 *y = prj->w[0]*sin(alpha);
2651 }
2652
2653 return 0;
2654 }
2655
2656 /*--------------------------------------------------------------------------*/
2657
molrev(x,y,prj,phi,theta)2658 int molrev(x, y, prj, phi, theta)
2659
2660 const double x, y;
2661 struct prjprm *prj;
2662 double *phi, *theta;
2663
2664 {
2665 double s, y0, z;
2666 const double tol = 1.0e-12;
2667
2668 if (prj->flag != PRJSET) {
2669 if (molset(prj)) return 1;
2670 }
2671
2672 y0 = y/prj->r0;
2673 s = 2.0 - y0*y0;
2674 if (s <= tol) {
2675 if (s < -tol) {
2676 return 2;
2677 }
2678 s = 0.0;
2679
2680 if (fabs(x) > tol) {
2681 return 2;
2682 }
2683 *phi = 0.0;
2684 } else {
2685 s = sqrt(s);
2686 *phi = prj->w[3]*x/s;
2687 }
2688
2689 z = y*prj->w[2];
2690 if (fabs(z) > 1.0) {
2691 if (fabs(z) > 1.0+tol) {
2692 return 2;
2693 }
2694 z = wcs_copysign(1.0,z) + y0*s/PI;
2695 } else {
2696 z = asin(z)*prj->w[4] + y0*s/PI;
2697 }
2698
2699 if (fabs(z) > 1.0) {
2700 if (fabs(z) > 1.0+tol) {
2701 return 2;
2702 }
2703 z = wcs_copysign(1.0,z);
2704 }
2705
2706 *theta = wcs_asind(z);
2707
2708 return 0;
2709 }
2710
2711 /*============================================================================
2712 * CSC: COBE quadrilateralized spherical cube projection.
2713 *
2714 * Given and/or returned:
2715 * prj->r0 r0; reset to 180/pi if 0.
2716 * prj->w[0] r0*(pi/4)
2717 * prj->w[1] (4/pi)/r0
2718 *===========================================================================*/
2719
cscset(prj)2720 int cscset(prj)
2721
2722 struct prjprm *prj;
2723
2724 {
2725 if (prj->r0 == 0.0) {
2726 prj->r0 = R2D;
2727 prj->w[0] = 45.0;
2728 prj->w[1] = 1.0/45.0;
2729 } else {
2730 prj->w[0] = prj->r0*PI/4.0;
2731 prj->w[1] = 1.0/prj->w[0];
2732 }
2733
2734 prj->flag = PRJSET;
2735
2736 return 0;
2737 }
2738
2739 /*--------------------------------------------------------------------------*/
2740
cscfwd(phi,theta,prj,x,y)2741 int cscfwd(phi, theta, prj, x, y)
2742
2743 const double phi, theta;
2744 struct prjprm *prj;
2745 double *x, *y;
2746
2747 {
2748 int face;
2749 double costhe, eta, l, m, n, rho, xi;
2750 const float tol = 1.0e-7;
2751
2752 float a, a2, a2b2, a4, ab, b, b2, b4, ca2, cb2, x0, xf, y0, yf;
2753 const float gstar = 1.37484847732;
2754 const float mm = 0.004869491981;
2755 const float gamma = -0.13161671474;
2756 const float omega1 = -0.159596235474;
2757 const float d0 = 0.0759196200467;
2758 const float d1 = -0.0217762490699;
2759 const float c00 = 0.141189631152;
2760 const float c10 = 0.0809701286525;
2761 const float c01 = -0.281528535557;
2762 const float c11 = 0.15384112876;
2763 const float c20 = -0.178251207466;
2764 const float c02 = 0.106959469314;
2765
2766 if (prj->flag != PRJSET) {
2767 if (cscset(prj)) return 1;
2768 }
2769
2770 costhe = wcs_cosd(theta);
2771 l = costhe*wcs_cosd(phi);
2772 m = costhe*wcs_sind(phi);
2773 n = wcs_sind(theta);
2774
2775 face = 0;
2776 rho = n;
2777 if (l > rho) {
2778 face = 1;
2779 rho = l;
2780 }
2781 if (m > rho) {
2782 face = 2;
2783 rho = m;
2784 }
2785 if (-l > rho) {
2786 face = 3;
2787 rho = -l;
2788 }
2789 if (-m > rho) {
2790 face = 4;
2791 rho = -m;
2792 }
2793 if (-n > rho) {
2794 face = 5;
2795 rho = -n;
2796 }
2797
2798 if (face == 0) {
2799 xi = m;
2800 eta = -l;
2801 x0 = 0.0;
2802 y0 = 2.0;
2803 } else if (face == 1) {
2804 xi = m;
2805 eta = n;
2806 x0 = 0.0;
2807 y0 = 0.0;
2808 } else if (face == 2) {
2809 xi = -l;
2810 eta = n;
2811 x0 = 2.0;
2812 y0 = 0.0;
2813 } else if (face == 3) {
2814 xi = -m;
2815 eta = n;
2816 x0 = 4.0;
2817 y0 = 0.0;
2818 } else if (face == 4) {
2819 xi = l;
2820 eta = n;
2821 x0 = 6.0;
2822 y0 = 0.0;
2823 } else {
2824 xi = m;
2825 eta = l;
2826 x0 = 0.0;
2827 y0 = -2.0;
2828 }
2829
2830 a = xi/rho;
2831 b = eta/rho;
2832
2833 a2 = a*a;
2834 b2 = b*b;
2835 ca2 = 1.0 - a2;
2836 cb2 = 1.0 - b2;
2837
2838 /* Avoid floating underflows. */
2839 ab = fabs(a*b);
2840 a4 = (a2 > 1.0e-16) ? a2*a2 : 0.0;
2841 b4 = (b2 > 1.0e-16) ? b2*b2 : 0.0;
2842 a2b2 = (ab > 1.0e-16) ? a2*b2 : 0.0;
2843
2844 xf = a*(a2 + ca2*(gstar + b2*(gamma*ca2 + mm*a2 +
2845 cb2*(c00 + c10*a2 + c01*b2 + c11*a2b2 + c20*a4 + c02*b4)) +
2846 a2*(omega1 - ca2*(d0 + d1*a2))));
2847 yf = b*(b2 + cb2*(gstar + a2*(gamma*cb2 + mm*b2 +
2848 ca2*(c00 + c10*b2 + c01*a2 + c11*a2b2 + c20*b4 + c02*a4)) +
2849 b2*(omega1 - cb2*(d0 + d1*b2))));
2850
2851 if (fabs(xf) > 1.0) {
2852 if (fabs(xf) > 1.0+tol) {
2853 return 2;
2854 }
2855 xf = wcs_copysign(1.0,xf);
2856 }
2857 if (fabs(yf) > 1.0) {
2858 if (fabs(yf) > 1.0+tol) {
2859 return 2;
2860 }
2861 yf = wcs_copysign(1.0,yf);
2862 }
2863
2864 *x = prj->w[0]*(x0 + xf);
2865 *y = prj->w[0]*(y0 + yf);
2866
2867 return 0;
2868 }
2869
2870 /*--------------------------------------------------------------------------*/
2871
cscrev(x,y,prj,phi,theta)2872 int cscrev(x, y, prj, phi, theta)
2873
2874 const double x, y;
2875 struct prjprm *prj;
2876 double *phi, *theta;
2877
2878 {
2879 int face;
2880 double l, m, n;
2881
2882 float a, b, xf, xx, yf, yy, z0, z1, z2, z3, z4, z5, z6;
2883 const float p00 = -0.27292696;
2884 const float p10 = -0.07629969;
2885 const float p20 = -0.22797056;
2886 const float p30 = 0.54852384;
2887 const float p40 = -0.62930065;
2888 const float p50 = 0.25795794;
2889 const float p60 = 0.02584375;
2890 const float p01 = -0.02819452;
2891 const float p11 = -0.01471565;
2892 const float p21 = 0.48051509;
2893 const float p31 = -1.74114454;
2894 const float p41 = 1.71547508;
2895 const float p51 = -0.53022337;
2896 const float p02 = 0.27058160;
2897 const float p12 = -0.56800938;
2898 const float p22 = 0.30803317;
2899 const float p32 = 0.98938102;
2900 const float p42 = -0.83180469;
2901 const float p03 = -0.60441560;
2902 const float p13 = 1.50880086;
2903 const float p23 = -0.93678576;
2904 const float p33 = 0.08693841;
2905 const float p04 = 0.93412077;
2906 const float p14 = -1.41601920;
2907 const float p24 = 0.33887446;
2908 const float p05 = -0.63915306;
2909 const float p15 = 0.52032238;
2910 const float p06 = 0.14381585;
2911
2912 if (prj->flag != PRJSET) {
2913 if (cscset(prj)) return 1;
2914 }
2915
2916 xf = x*prj->w[1];
2917 yf = y*prj->w[1];
2918
2919 /* Check bounds. */
2920 if (fabs(xf) <= 1.0) {
2921 if (fabs(yf) > 3.0) return 2;
2922 } else {
2923 if (fabs(xf) > 7.0) return 2;
2924 if (fabs(yf) > 1.0) return 2;
2925 }
2926
2927 /* Map negative faces to the other side. */
2928 if (xf < -1.0) xf += 8.0;
2929
2930 /* Determine the face. */
2931 if (xf > 5.0) {
2932 face = 4;
2933 xf = xf - 6.0;
2934 } else if (xf > 3.0) {
2935 face = 3;
2936 xf = xf - 4.0;
2937 } else if (xf > 1.0) {
2938 face = 2;
2939 xf = xf - 2.0;
2940 } else if (yf > 1.0) {
2941 face = 0;
2942 yf = yf - 2.0;
2943 } else if (yf < -1.0) {
2944 face = 5;
2945 yf = yf + 2.0;
2946 } else {
2947 face = 1;
2948 }
2949
2950 xx = xf*xf;
2951 yy = yf*yf;
2952
2953 z0 = p00 + xx*(p10 + xx*(p20 + xx*(p30 + xx*(p40 + xx*(p50 + xx*(p60))))));
2954 z1 = p01 + xx*(p11 + xx*(p21 + xx*(p31 + xx*(p41 + xx*(p51)))));
2955 z2 = p02 + xx*(p12 + xx*(p22 + xx*(p32 + xx*(p42))));
2956 z3 = p03 + xx*(p13 + xx*(p23 + xx*(p33)));
2957 z4 = p04 + xx*(p14 + xx*(p24));
2958 z5 = p05 + xx*(p15);
2959 z6 = p06;
2960
2961 a = z0 + yy*(z1 + yy*(z2 + yy*(z3 + yy*(z4 + yy*(z5 + yy*z6)))));
2962 a = xf + xf*(1.0 - xx)*a;
2963
2964 z0 = p00 + yy*(p10 + yy*(p20 + yy*(p30 + yy*(p40 + yy*(p50 + yy*(p60))))));
2965 z1 = p01 + yy*(p11 + yy*(p21 + yy*(p31 + yy*(p41 + yy*(p51)))));
2966 z2 = p02 + yy*(p12 + yy*(p22 + yy*(p32 + yy*(p42))));
2967 z3 = p03 + yy*(p13 + yy*(p23 + yy*(p33)));
2968 z4 = p04 + yy*(p14 + yy*(p24));
2969 z5 = p05 + yy*(p15);
2970 z6 = p06;
2971
2972 b = z0 + xx*(z1 + xx*(z2 + xx*(z3 + xx*(z4 + xx*(z5 + xx*z6)))));
2973 b = yf + yf*(1.0 - yy)*b;
2974
2975 if (face == 0) {
2976 n = 1.0/sqrt(a*a + b*b + 1.0);
2977 l = -b*n;
2978 m = a*n;
2979 } else if (face == 1) {
2980 l = 1.0/sqrt(a*a + b*b + 1.0);
2981 m = a*l;
2982 n = b*l;
2983 } else if (face == 2) {
2984 m = 1.0/sqrt(a*a + b*b + 1.0);
2985 l = -a*m;
2986 n = b*m;
2987 } else if (face == 3) {
2988 l = -1.0/sqrt(a*a + b*b + 1.0);
2989 m = a*l;
2990 n = -b*l;
2991 } else if (face == 4) {
2992 m = -1.0/sqrt(a*a + b*b + 1.0);
2993 l = -a*m;
2994 n = -b*m;
2995 } else {
2996 n = -1.0/sqrt(a*a + b*b + 1.0);
2997 l = -b*n;
2998 m = -a*n;
2999 }
3000
3001 if (l == 0.0 && m == 0.0) {
3002 *phi = 0.0;
3003 } else {
3004 *phi = wcs_atan2d(m, l);
3005 }
3006 *theta = wcs_asind(n);
3007
3008 return 0;
3009 }
3010
3011 /*============================================================================
3012 * QSC: quadrilaterilized spherical cube projection.
3013 *
3014 * Given and/or returned:
3015 * prj->r0 r0; reset to 180/pi if 0.
3016 * prj->w[0] r0*(pi/4)
3017 * prj->w[1] (4/pi)/r0
3018 *===========================================================================*/
3019
qscset(prj)3020 int qscset(prj)
3021
3022 struct prjprm *prj;
3023
3024 {
3025 if (prj->r0 == 0.0) {
3026 prj->r0 = R2D;
3027 prj->w[0] = 45.0;
3028 prj->w[1] = 1.0/45.0;
3029 } else {
3030 prj->w[0] = prj->r0*PI/4.0;
3031 prj->w[1] = 1.0/prj->w[0];
3032 }
3033
3034 prj->flag = PRJSET;
3035
3036 return 0;
3037 }
3038
3039 /*--------------------------------------------------------------------------*/
3040
qscfwd(phi,theta,prj,x,y)3041 int qscfwd(phi, theta, prj, x, y)
3042
3043 const double phi, theta;
3044 struct prjprm *prj;
3045 double *x, *y;
3046
3047 {
3048 int face;
3049 double chi, costhe, eta, l, m, n, p, psi, rho, rhu, t, x0, xf, xi, y0, yf;
3050 const double tol = 1.0e-12;
3051
3052 if (prj->flag != PRJSET) {
3053 if (qscset(prj)) return 1;
3054 }
3055
3056 if (fabs(theta) == 90.0) {
3057 *x = 0.0;
3058 *y = wcs_copysign(2.0*prj->w[0],theta);
3059 return 0;
3060 }
3061
3062 costhe = wcs_cosd(theta);
3063 l = costhe*wcs_cosd(phi);
3064 m = costhe*wcs_sind(phi);
3065 n = wcs_sind(theta);
3066
3067 face = 0;
3068 rho = n;
3069 if (l > rho) {
3070 face = 1;
3071 rho = l;
3072 }
3073 if (m > rho) {
3074 face = 2;
3075 rho = m;
3076 }
3077 if (-l > rho) {
3078 face = 3;
3079 rho = -l;
3080 }
3081 if (-m > rho) {
3082 face = 4;
3083 rho = -m;
3084 }
3085 if (-n > rho) {
3086 face = 5;
3087 rho = -n;
3088 }
3089
3090 rhu = 1.0 - rho;
3091
3092 if (face == 0) {
3093 xi = m;
3094 eta = -l;
3095 if (rhu < 1.0e-8) {
3096 /* Small angle formula. */
3097 t = (90.0 - theta)*D2R;
3098 rhu = t*t/2.0;
3099 }
3100 x0 = 0.0;
3101 y0 = 2.0;
3102 } else if (face == 1) {
3103 xi = m;
3104 eta = n;
3105 if (rhu < 1.0e-8) {
3106 /* Small angle formula. */
3107 t = theta*D2R;
3108 p = fmod(phi,360.0);
3109 if (p < -180.0) p += 360.0;
3110 if (p > 180.0) p -= 360.0;
3111 p *= D2R;
3112 rhu = (p*p + t*t)/2.0;
3113 }
3114 x0 = 0.0;
3115 y0 = 0.0;
3116 } else if (face == 2) {
3117 xi = -l;
3118 eta = n;
3119 if (rhu < 1.0e-8) {
3120 /* Small angle formula. */
3121 t = theta*D2R;
3122 p = fmod(phi,360.0);
3123 if (p < -180.0) p += 360.0;
3124 p = (90.0 - p)*D2R;
3125 rhu = (p*p + t*t)/2.0;
3126 }
3127 x0 = 2.0;
3128 y0 = 0.0;
3129 } else if (face == 3) {
3130 xi = -m;
3131 eta = n;
3132 if (rhu < 1.0e-8) {
3133 /* Small angle formula. */
3134 t = theta*D2R;
3135 p = fmod(phi,360.0);
3136 if (p < 0.0) p += 360.0;
3137 p = (180.0 - p)*D2R;
3138 rhu = (p*p + t*t)/2.0;
3139 }
3140 x0 = 4.0;
3141 y0 = 0.0;
3142 } else if (face == 4) {
3143 xi = l;
3144 eta = n;
3145 if (rhu < 1.0e-8) {
3146 /* Small angle formula. */
3147 t = theta*D2R;
3148 p = fmod(phi,360.0);
3149 if (p > 180.0) p -= 360.0;
3150 p *= (90.0 + p)*D2R;
3151 rhu = (p*p + t*t)/2.0;
3152 }
3153 x0 = 6;
3154 y0 = 0.0;
3155 } else {
3156 xi = m;
3157 eta = l;
3158 if (rhu < 1.0e-8) {
3159 /* Small angle formula. */
3160 t = (90.0 + theta)*D2R;
3161 rhu = t*t/2.0;
3162 }
3163 x0 = 0.0;
3164 y0 = -2;
3165 }
3166
3167 if (xi == 0.0 && eta == 0.0) {
3168 xf = 0.0;
3169 yf = 0.0;
3170 } else if (-xi >= fabs(eta)) {
3171 psi = eta/xi;
3172 chi = 1.0 + psi*psi;
3173 xf = -sqrt(rhu/(1.0-1.0/sqrt(1.0+chi)));
3174 yf = (xf/15.0)*(wcs_atand(psi) - wcs_asind(psi/sqrt(chi+chi)));
3175 } else if (xi >= fabs(eta)) {
3176 psi = eta/xi;
3177 chi = 1.0 + psi*psi;
3178 xf = sqrt(rhu/(1.0-1.0/sqrt(1.0+chi)));
3179 yf = (xf/15.0)*(wcs_atand(psi) - wcs_asind(psi/sqrt(chi+chi)));
3180 } else if (-eta > fabs(xi)) {
3181 psi = xi/eta;
3182 chi = 1.0 + psi*psi;
3183 yf = -sqrt(rhu/(1.0-1.0/sqrt(1.0+chi)));
3184 xf = (yf/15.0)*(wcs_atand(psi) - wcs_asind(psi/sqrt(chi+chi)));
3185 } else {
3186 psi = xi/eta;
3187 chi = 1.0 + psi*psi;
3188 yf = sqrt(rhu/(1.0-1.0/sqrt(1.0+chi)));
3189 xf = (yf/15.0)*(wcs_atand(psi) - wcs_asind(psi/sqrt(chi+chi)));
3190 }
3191
3192 if (fabs(xf) > 1.0) {
3193 if (fabs(xf) > 1.0+tol) {
3194 return 2;
3195 }
3196 xf = wcs_copysign(1.0,xf);
3197 }
3198 if (fabs(yf) > 1.0) {
3199 if (fabs(yf) > 1.0+tol) {
3200 return 2;
3201 }
3202 yf = wcs_copysign(1.0,yf);
3203 }
3204
3205 *x = prj->w[0]*(xf + x0);
3206 *y = prj->w[0]*(yf + y0);
3207
3208
3209 return 0;
3210 }
3211
3212 /*--------------------------------------------------------------------------*/
3213
qscrev(x,y,prj,phi,theta)3214 int qscrev(x, y, prj, phi, theta)
3215
3216 const double x, y;
3217 struct prjprm *prj;
3218 double *phi, *theta;
3219
3220 {
3221 int direct, face;
3222 double chi, l, m, n, psi, rho, rhu, xf, yf, w;
3223 const double tol = 1.0e-12;
3224
3225 if (prj->flag != PRJSET) {
3226 if (qscset(prj)) return 1;
3227 }
3228
3229 xf = x*prj->w[1];
3230 yf = y*prj->w[1];
3231
3232 /* Check bounds. */
3233 if (fabs(xf) <= 1.0) {
3234 if (fabs(yf) > 3.0) return 2;
3235 } else {
3236 if (fabs(xf) > 7.0) return 2;
3237 if (fabs(yf) > 1.0) return 2;
3238 }
3239
3240 /* Map negative faces to the other side. */
3241 if (xf < -1.0) xf += 8.0;
3242
3243 /* Determine the face. */
3244 if (xf > 5.0) {
3245 face = 4;
3246 xf = xf - 6.0;
3247 } else if (xf > 3.0) {
3248 face = 3;
3249 xf = xf - 4.0;
3250 } else if (xf > 1.0) {
3251 face = 2;
3252 xf = xf - 2.0;
3253 } else if (yf > 1.0) {
3254 face = 0;
3255 yf = yf - 2.0;
3256 } else if (yf < -1.0) {
3257 face = 5;
3258 yf = yf + 2.0;
3259 } else {
3260 face = 1;
3261 }
3262
3263 direct = (fabs(xf) > fabs(yf));
3264 if (direct) {
3265 if (xf == 0.0) {
3266 psi = 0.0;
3267 chi = 1.0;
3268 rho = 1.0;
3269 rhu = 0.0;
3270 } else {
3271 w = 15.0*yf/xf;
3272 psi = wcs_sind(w)/(wcs_cosd(w) - SQRT2INV);
3273 chi = 1.0 + psi*psi;
3274 rhu = xf*xf*(1.0 - 1.0/sqrt(1.0 + chi));
3275 rho = 1.0 - rhu;
3276 }
3277 } else {
3278 if (yf == 0.0) {
3279 psi = 0.0;
3280 chi = 1.0;
3281 rho = 1.0;
3282 rhu = 0.0;
3283 } else {
3284 w = 15.0*xf/yf;
3285 psi = wcs_sind(w)/(wcs_cosd(w) - SQRT2INV);
3286 chi = 1.0 + psi*psi;
3287 rhu = yf*yf*(1.0 - 1.0/sqrt(1.0 + chi));
3288 rho = 1.0 - rhu;
3289 }
3290 }
3291
3292 if (rho < -1.0) {
3293 if (rho < -1.0-tol) {
3294 return 2;
3295 }
3296
3297 rho = -1.0;
3298 rhu = 2.0;
3299 w = 0.0;
3300 } else {
3301 w = sqrt(rhu*(2.0-rhu)/chi);
3302 }
3303
3304 if (face == 0) {
3305 n = rho;
3306 if (direct) {
3307 m = w;
3308 if (xf < 0.0) m = -m;
3309 l = -m*psi;
3310 } else {
3311 l = w;
3312 if (yf > 0.0) l = -l;
3313 m = -l*psi;
3314 }
3315 } else if (face == 1) {
3316 l = rho;
3317 if (direct) {
3318 m = w;
3319 if (xf < 0.0) m = -m;
3320 n = m*psi;
3321 } else {
3322 n = w;
3323 if (yf < 0.0) n = -n;
3324 m = n*psi;
3325 }
3326 } else if (face == 2) {
3327 m = rho;
3328 if (direct) {
3329 l = w;
3330 if (xf > 0.0) l = -l;
3331 n = -l*psi;
3332 } else {
3333 n = w;
3334 if (yf < 0.0) n = -n;
3335 l = -n*psi;
3336 }
3337 } else if (face == 3) {
3338 l = -rho;
3339 if (direct) {
3340 m = w;
3341 if (xf > 0.0) m = -m;
3342 n = -m*psi;
3343 } else {
3344 n = w;
3345 if (yf < 0.0) n = -n;
3346 m = -n*psi;
3347 }
3348 } else if (face == 4) {
3349 m = -rho;
3350 if (direct) {
3351 l = w;
3352 if (xf < 0.0) l = -l;
3353 n = l*psi;
3354 } else {
3355 n = w;
3356 if (yf < 0.0) n = -n;
3357 l = n*psi;
3358 }
3359 } else {
3360 n = -rho;
3361 if (direct) {
3362 m = w;
3363 if (xf < 0.0) m = -m;
3364 l = m*psi;
3365 } else {
3366 l = w;
3367 if (yf < 0.0) l = -l;
3368 m = l*psi;
3369 }
3370 }
3371
3372 if (l == 0.0 && m == 0.0) {
3373 *phi = 0.0;
3374 } else {
3375 *phi = wcs_atan2d(m, l);
3376 }
3377 *theta = wcs_asind(n);
3378
3379 return 0;
3380 }
3381
3382 /*============================================================================
3383 * TSC: tangential spherical cube projection.
3384 *
3385 * Given and/or returned:
3386 * prj->r0 r0; reset to 180/pi if 0.
3387 * prj->w[0] r0*(pi/4)
3388 * prj->w[1] (4/pi)/r0
3389 *===========================================================================*/
3390
tscset(prj)3391 int tscset(prj)
3392
3393 struct prjprm *prj;
3394
3395 {
3396 if (prj->r0 == 0.0) {
3397 prj->r0 = R2D;
3398 prj->w[0] = 45.0;
3399 prj->w[1] = 1.0/45.0;
3400 } else {
3401 prj->w[0] = prj->r0*PI/4.0;
3402 prj->w[1] = 1.0/prj->w[0];
3403 }
3404
3405 prj->flag = PRJSET;
3406
3407 return 0;
3408 }
3409
3410 /*--------------------------------------------------------------------------*/
3411
tscfwd(phi,theta,prj,x,y)3412 int tscfwd(phi, theta, prj, x, y)
3413
3414 const double phi, theta;
3415 struct prjprm *prj;
3416 double *x, *y;
3417
3418 {
3419 int face;
3420 double costhe, l, m, n, rho, x0, xf, y0, yf;
3421 const double tol = 1.0e-12;
3422
3423 if (prj->flag != PRJSET) {
3424 if (tscset(prj)) return 1;
3425 }
3426
3427 costhe = wcs_cosd(theta);
3428 l = costhe*wcs_cosd(phi);
3429 m = costhe*wcs_sind(phi);
3430 n = wcs_sind(theta);
3431
3432 face = 0;
3433 rho = n;
3434 if (l > rho) {
3435 face = 1;
3436 rho = l;
3437 }
3438 if (m > rho) {
3439 face = 2;
3440 rho = m;
3441 }
3442 if (-l > rho) {
3443 face = 3;
3444 rho = -l;
3445 }
3446 if (-m > rho) {
3447 face = 4;
3448 rho = -m;
3449 }
3450 if (-n > rho) {
3451 face = 5;
3452 rho = -n;
3453 }
3454
3455 if (face == 0) {
3456 xf = m/rho;
3457 yf = -l/rho;
3458 x0 = 0.0;
3459 y0 = 2.0;
3460 } else if (face == 1) {
3461 xf = m/rho;
3462 yf = n/rho;
3463 x0 = 0.0;
3464 y0 = 0.0;
3465 } else if (face == 2) {
3466 xf = -l/rho;
3467 yf = n/rho;
3468 x0 = 2.0;
3469 y0 = 0.0;
3470 } else if (face == 3) {
3471 xf = -m/rho;
3472 yf = n/rho;
3473 x0 = 4.0;
3474 y0 = 0.0;
3475 } else if (face == 4) {
3476 xf = l/rho;
3477 yf = n/rho;
3478 x0 = 6.0;
3479 y0 = 0.0;
3480 } else {
3481 xf = m/rho;
3482 yf = l/rho;
3483 x0 = 0.0;
3484 y0 = -2.0;
3485 }
3486
3487 if (fabs(xf) > 1.0) {
3488 if (fabs(xf) > 1.0+tol) {
3489 return 2;
3490 }
3491 xf = wcs_copysign(1.0,xf);
3492 }
3493 if (fabs(yf) > 1.0) {
3494 if (fabs(yf) > 1.0+tol) {
3495 return 2;
3496 }
3497 yf = wcs_copysign(1.0,yf);
3498 }
3499
3500 *x = prj->w[0]*(xf + x0);
3501 *y = prj->w[0]*(yf + y0);
3502
3503 return 0;
3504 }
3505
3506 /*--------------------------------------------------------------------------*/
3507
tscrev(x,y,prj,phi,theta)3508 int tscrev(x, y, prj, phi, theta)
3509
3510 const double x, y;
3511 struct prjprm *prj;
3512 double *phi, *theta;
3513
3514 {
3515 double l, m, n, xf, yf;
3516
3517 if (prj->flag != PRJSET) {
3518 if (tscset(prj)) return 1;
3519 }
3520
3521 xf = x*prj->w[1];
3522 yf = y*prj->w[1];
3523
3524 /* Check bounds. */
3525 if (fabs(xf) <= 1.0) {
3526 if (fabs(yf) > 3.0) return 2;
3527 } else {
3528 if (fabs(xf) > 7.0) return 2;
3529 if (fabs(yf) > 1.0) return 2;
3530 }
3531
3532 /* Map negative faces to the other side. */
3533 if (xf < -1.0) xf += 8.0;
3534
3535 /* Determine the face. */
3536 if (xf > 5.0) {
3537 /* face = 4 */
3538 xf = xf - 6.0;
3539 m = -1.0/sqrt(1.0 + xf*xf + yf*yf);
3540 l = -m*xf;
3541 n = -m*yf;
3542 } else if (xf > 3.0) {
3543 /* face = 3 */
3544 xf = xf - 4.0;
3545 l = -1.0/sqrt(1.0 + xf*xf + yf*yf);
3546 m = l*xf;
3547 n = -l*yf;
3548 } else if (xf > 1.0) {
3549 /* face = 2 */
3550 xf = xf - 2.0;
3551 m = 1.0/sqrt(1.0 + xf*xf + yf*yf);
3552 l = -m*xf;
3553 n = m*yf;
3554 } else if (yf > 1.0) {
3555 /* face = 0 */
3556 yf = yf - 2.0;
3557 n = 1.0/sqrt(1.0 + xf*xf + yf*yf);
3558 l = -n*yf;
3559 m = n*xf;
3560 } else if (yf < -1.0) {
3561 /* face = 5 */
3562 yf = yf + 2.0;
3563 n = -1.0/sqrt(1.0 + xf*xf + yf*yf);
3564 l = -n*yf;
3565 m = -n*xf;
3566 } else {
3567 /* face = 1 */
3568 l = 1.0/sqrt(1.0 + xf*xf + yf*yf);
3569 m = l*xf;
3570 n = l*yf;
3571 }
3572
3573 if (l == 0.0 && m == 0.0) {
3574 *phi = 0.0;
3575 } else {
3576 *phi = wcs_atan2d(m, l);
3577 }
3578 *theta = wcs_asind(n);
3579
3580 return 0;
3581 }
3582
3583
3584 /*============================================================================
3585 * TNX: IRAF's gnomonic projection.
3586 *
3587 * Given and/or returned:
3588 * prj->r0 r0; reset to 180/pi if 0.
3589 *===========================================================================*/
3590
tnxset(prj)3591 int tnxset(prj)
3592
3593 struct prjprm *prj;
3594
3595 {
3596 if (prj->r0 == 0.0) prj->r0 = R2D;
3597
3598 if (prj->flag == -1) {
3599 prj->flag = -PRJSET;
3600 } else {
3601 prj->flag = PRJSET;
3602 }
3603
3604 return 0;
3605 }
3606
3607 /*--------------------------------------------------------------------------*/
3608
tnxfwd(phi,theta,prj,x,y)3609 int tnxfwd(phi, theta, prj, x, y)
3610
3611 const double phi, theta;
3612 struct prjprm *prj;
3613 double *x, *y;
3614
3615 {
3616 double r, s, xp[2];
3617
3618 if (abs(prj->flag) != PRJSET) {
3619 if(tnxset(prj)) return 1;
3620 }
3621
3622 s = wcs_sind(theta);
3623 if (s == 0.0) return 2;
3624
3625 r = prj->r0*wcs_cosd(theta)/s;
3626 xp[0] = r*wcs_sind(phi);
3627 xp[1] = -r*wcs_cosd(phi);
3628 *x = prj->inv_x? poly_func(prj->inv_x, xp) : xp[0];
3629 *y = prj->inv_y? poly_func(prj->inv_y, xp) : xp[1];
3630
3631 if (prj->flag == PRJSET && s < 0.0) {
3632 return 2;
3633 }
3634
3635 return 0;
3636 }
3637
3638 /*--------------------------------------------------------------------------*/
3639
tnxrev(x,y,prj,phi,theta)3640 int tnxrev(x, y, prj, phi, theta)
3641
3642 const double x, y;
3643 struct prjprm *prj;
3644 double *phi, *theta;
3645
3646 {
3647 double rp,xp,yp;
3648
3649 if (abs(prj->flag) != PRJSET) {
3650 if (tanset(prj)) return 1;
3651 }
3652
3653 xp = x+raw_to_tnxaxis(prj->tnx_lngcor, x, y);
3654 yp = y+raw_to_tnxaxis(prj->tnx_latcor, x, y);
3655 if ((rp = sqrt(xp*xp+yp*yp)) == 0.0) {
3656 *phi = 0.0;
3657 } else {
3658 *phi = wcs_atan2d(xp, -yp);
3659 }
3660 *theta = wcs_atan2d(prj->r0, rp);
3661
3662 return 0;
3663 }
3664
3665 /*--------------------------------------------------------------------------*/
3666
raw_to_pv(struct prjprm * prj,double x,double y,double * xo,double * yo)3667 int raw_to_pv(struct prjprm *prj, double x, double y, double *xo, double *yo)
3668
3669 {
3670 int k;
3671 double *a,*b,
3672 r,r3,r5,r7,xy,x2,x3,x4,x5,x6,x7,y2,y3,y4,y5,y6,y7,xp,yp;
3673
3674 if (abs(prj->flag) != PRJSET) {
3675 if (tanset(prj)) return 1;
3676 }
3677
3678 k=prj->n;
3679 a = prj->p; /* Longitude */
3680 b = prj->p+100; /* Latitude */
3681 xp = *(a++);
3682 xp += *(a++)*x;
3683 yp = *(b++);
3684 yp += *(b++)*y;
3685 if (!--k) goto poly_end;
3686 xp += *(a++)*y;
3687 yp += *(b++)*x;
3688 if (!--k) goto poly_end;
3689 r = sqrt(x*x + y*y);
3690 xp += *(a++)*r;
3691 yp += *(b++)*r;
3692 if (!--k) goto poly_end;
3693 xp += *(a++)*(x2=x*x);
3694 yp += *(b++)*(y2=y*y);
3695 if (!--k) goto poly_end;
3696 xp += *(a++)*(xy=x*y);
3697 yp += *(b++)*xy;
3698 if (!--k) goto poly_end;
3699 xp += *(a++)*y2;
3700 yp += *(b++)*x2;
3701 if (!--k) goto poly_end;
3702 xp += *(a++)*(x3=x*x2);
3703 yp += *(b++)*(y3=y*y2);
3704 if (!--k) goto poly_end;
3705 xp += *(a++)*x2*y;
3706 yp += *(b++)*y2*x;
3707 if (!--k) goto poly_end;
3708 xp += *(a++)*x*y2;
3709 yp += *(b++)*y*x2;
3710 if (!--k) goto poly_end;
3711 xp += *(a++)*y3;
3712 yp += *(b++)*x3;
3713 if (!--k) goto poly_end;
3714 xp += *(a++)*(r3=r*r*r);
3715 yp += *(b++)*r3;
3716 if (!--k) goto poly_end;
3717 xp += *(a++)*(x4=x2*x2);
3718 yp += *(b++)*(y4=y2*y2);
3719 if (!--k) goto poly_end;
3720 xp += *(a++)*x3*y;
3721 yp += *(b++)*y3*x;
3722 if (!--k) goto poly_end;
3723 xp += *(a++)*x2*y2;
3724 yp += *(b++)*x2*y2;
3725 if (!--k) goto poly_end;
3726 xp += *(a++)*x*y3;
3727 yp += *(b++)*y*x3;
3728 if (!--k) goto poly_end;
3729 xp += *(a++)*y4;
3730 yp += *(b++)*x4;
3731 if (!--k) goto poly_end;
3732 xp += *(a++)*(x5=x4*x);
3733 yp += *(b++)*(y5=y4*y);
3734 if (!--k) goto poly_end;
3735 xp += *(a++)*x4*y;
3736 yp += *(b++)*y4*x;
3737 if (!--k) goto poly_end;
3738 xp += *(a++)*x3*y2;
3739 yp += *(b++)*y3*x2;
3740 if (!--k) goto poly_end;
3741 xp += *(a++)*x2*y3;
3742 yp += *(b++)*y2*x3;
3743 if (!--k) goto poly_end;
3744 xp += *(a++)*x*y4;
3745 yp += *(b++)*y*x4;
3746 if (!--k) goto poly_end;
3747 xp += *(a++)*y5;
3748 yp += *(b++)*x5;
3749 if (!--k) goto poly_end;
3750 xp += *(a++)*(r5=r3*r*r);
3751 yp += *(b++)*r5;
3752 if (!--k) goto poly_end;
3753 xp += *(a++)*(x6=x5*x);
3754 yp += *(b++)*(y6=y5*y);
3755 if (!--k) goto poly_end;
3756 xp += *(a++)*x5*y;
3757 yp += *(b++)*y5*x;
3758 if (!--k) goto poly_end;
3759 xp += *(a++)*x4*y2;
3760 yp += *(b++)*y4*x2;
3761 if (!--k) goto poly_end;
3762 xp += *(a++)*x3*y3;
3763 yp += *(b++)*y3*x3;
3764 if (!--k) goto poly_end;
3765 xp += *(a++)*x2*y4;
3766 yp += *(b++)*y4*x2;
3767 if (!--k) goto poly_end;
3768 xp += *(a++)*x*y5;
3769 yp += *(b++)*y*x5;
3770 if (!--k) goto poly_end;
3771 xp += *(a++)*y6;
3772 yp += *(b++)*x6;
3773 if (!--k) goto poly_end;
3774 xp += *(a++)*(x7=x6*x);
3775 yp += *(b++)*(y7=y6*y);
3776 if (!--k) goto poly_end;
3777 xp += *(a++)*x6*y;
3778 yp += *(b++)*y6*x;
3779 if (!--k) goto poly_end;
3780 xp += *(a++)*x5*y2;
3781 yp += *(b++)*y5*x2;
3782 if (!--k) goto poly_end;
3783 xp += *(a++)*x4*y3;
3784 yp += *(b++)*y4*x3;
3785 if (!--k) goto poly_end;
3786 xp += *(a++)*x3*y4;
3787 yp += *(b++)*y3*x4;
3788 if (!--k) goto poly_end;
3789 xp += *(a++)*x2*y5;
3790 yp += *(b++)*y2*x5;
3791 if (!--k) goto poly_end;
3792 xp += *(a++)*x*y6;
3793 yp += *(b++)*y*x6;
3794 if (!--k) goto poly_end;
3795 xp += *(a++)*y7;
3796 yp += *(b++)*x7;
3797 if (!--k) goto poly_end;
3798 xp += *a*(r7=r5*r*r);
3799 yp += *b*r7;
3800
3801 poly_end:
3802
3803 *xo = xp;
3804 *yo = yp;
3805
3806 return 0;
3807 }
3808
3809