1 /*============================================================================
2 *
3 * WCSLIB - an implementation of the FITS WCS proposal.
4 * Copyright (C) 1995-2002, 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., 51 Franklin Street,Fifth Floor, Boston, MA 02110-1301, 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 *
31 * This version of proj.c is based on the version in wcslib-2.9, but has
32 * been modified in the following ways by the Starlink project (e-mail:
33 * ussc@star.rl.ac.uk):
34 * - The copysign macro is now always defined within this file
35 * instead of only being defined if the COPYSIGN macro has previously
36 * been defined.
37 * - Sine values which are slightly larger than 1.0 are now treated
38 * as 1.0 in function astCYPrev.
39 * - The maximum number of projection parameters has been changed from
40 * 10 to 100.
41 * - The maximum number of projection parameters is given by the
42 * WCSLIB_MXPAR macro (defined in proj.h) instead of being hard-wired.
43 * - The names of all functions and structures have been chanegd to avoid
44 * clashes with wcslib. This involves adding "Ast" or "ast" at the
45 * front and changing the capitalisation.
46 * - Include string.h (for strcpy and strcmp prototypes).
47 * - Include stdlib.h (for abs prototype).
48 * - Comment out declarations of npcode and pcodes variables (they
49 * are not needed by AST) in order to avoid clash with similar names
50 * in other modules imported as part of other software systems (e.g.
51 * SkyCat).
52 * - astZPNfwd: Loop from prj->n to zero, not from MAXPAR to zero.
53 * - astZPNfwd: Only return "2" if prj->n is larger than 2.
54 * - Lots of variables are initialised to null values in order to
55 * avoid "use of uninitialised variable" messages from compilers which
56 * are not clever enough to work out that the uninitialised variable is
57 * not in fact ever used.
58 * - Use dynamic rather than static memory for the parameter arrays in
59 * the AstPrjPrm structure.Override astGetObjSize. This is to
60 * reduce the in-memory size of a WcsMap.
61 * - HPX and XPH projections included from a more recent version of WCSLIB,
62 * and modified to use scalar instead of vector positions
63 * - The expressions for xc in astHPXrev and phic in astHPXfwd have
64 * been conditioned differently to the WCSLIB code in order to improve
65 * accuracy of the floor function for arguments very slightly below an
66 * integer value.
67
68 *=============================================================================
69 *
70 * C implementation of the spherical map projections recognized by the FITS
71 * "World Coordinate System" (WCS) convention.
72 *
73 * Summary of routines
74 * -------------------
75 * Each projection is implemented via separate functions for the forward,
76 * *fwd(), and reverse, *rev(), transformation.
77 *
78 * Initialization routines, *set(), compute intermediate values from the
79 * projection parameters but need not be called explicitly - see the
80 * explanation of prj.flag below.
81 *
82 * astPRJset astPRJfwd astPRJrev Driver routines (see below).
83 *
84 * astAZPset astAZPfwd astAZPrev AZP: zenithal/azimuthal perspective
85 * astSZPset astSZPfwd astSZPrev SZP: slant zenithal perspective
86 * astTANset astTANfwd astTANrev TAN: gnomonic
87 * astSTGset astSTGfwd astSTGrev STG: stereographic
88 * astSINset astSINfwd astSINrev SIN: orthographic/synthesis
89 * astARCset astARCfwd astARCrev ARC: zenithal/azimuthal equidistant
90 * astZPNset astZPNfwd astZPNrev ZPN: zenithal/azimuthal polynomial
91 * astZEAset astZEAfwd astZEArev ZEA: zenithal/azimuthal equal area
92 * astAIRset astAIRfwd astAIRrev AIR: Airy
93 * astCYPset astCYPfwd astCYPrev CYP: cylindrical perspective
94 * astCEAset astCEAfwd astCEArev CEA: cylindrical equal area
95 * astCARset astCARfwd astCARrev CAR: Cartesian
96 * astMERset astMERfwd astMERrev MER: Mercator
97 * astSFLset astSFLfwd astSFLrev SFL: Sanson-Flamsteed
98 * astPARset astPARfwd astPARrev PAR: parabolic
99 * astMOLset astMOLfwd astMOLrev MOL: Mollweide
100 * astAITset astAITfwd astAITrev AIT: Hammer-Aitoff
101 * astCOPset astCOPfwd astCOPrev COP: conic perspective
102 * astCOEset astCOEfwd astCOErev COE: conic equal area
103 * astCODset astCODfwd astCODrev COD: conic equidistant
104 * astCOOset astCOOfwd astCOOrev COO: conic orthomorphic
105 * astBONset astBONfwd astBONrev BON: Bonne
106 * astPCOset astPCOfwd astPCOrev PCO: polyconic
107 * astTSCset astTSCfwd astTSCrev TSC: tangential spherical cube
108 * astCSCset astCSCfwd astCSCrev CSC: COBE quadrilateralized spherical cube
109 * astQSCset astQSCfwd astQSCrev QSC: quadrilateralized spherical cube
110 * astHPXset astHPXfwd astHPXrev HPX: HEALPix projection
111 * astXPHset astXPHfwd astXPHrev XPH: HEALPix polar, aka "butterfly"
112 *
113 *
114 * Driver routines; astPRJset(), astPRJfwd() & astPRJrev()
115 * ----------------------------------------------
116 * A set of driver routines are available for use as a generic interface to
117 * the specific projection routines. The interfaces to astPRJfwd() and astPRJrev()
118 * are the same as those of the forward and reverse transformation routines
119 * for the specific projections (see below).
120 *
121 * The interface to astPRJset() differs slightly from that of the initialization
122 * routines for the specific projections and unlike them it must be invoked
123 * explicitly to use astPRJfwd() and astPRJrev().
124 *
125 * Given:
126 * pcode[4] const char
127 * WCS projection code.
128 *
129 * Given and/or returned:
130 * prj AstPrjPrm* Projection parameters (see below).
131 *
132 * Function return value:
133 * int Error status
134 * 0: Success.
135 *
136 *
137 * Initialization routine; *set()
138 * ------------------------------
139 * Initializes members of a AstPrjPrm data structure which hold intermediate
140 * values. Note that this routine need not be called directly; it will be
141 * invoked by astPRJfwd() and astPRJrev() if the "flag" structure member is
142 * anything other than a predefined magic value.
143 *
144 * Given and/or returned:
145 * prj AstPrjPrm* Projection parameters (see below).
146 *
147 * Function return value:
148 * int Error status
149 * 0: Success.
150 * 1: Invalid projection parameters.
151 *
152 * Forward transformation; *fwd()
153 * -----------------------------
154 * Compute (x,y) coordinates in the plane of projection from native spherical
155 * coordinates (phi,theta).
156 *
157 * Given:
158 * phi, const double
159 * theta Longitude and latitude of the projected point in
160 * native spherical coordinates, in degrees.
161 *
162 * Given and returned:
163 * prj AstPrjPrm* Projection parameters (see below).
164 *
165 * Returned:
166 * x,y double* Projected coordinates.
167 *
168 * Function return value:
169 * int Error status
170 * 0: Success.
171 * 1: Invalid projection parameters.
172 * 2: Invalid value of (phi,theta).
173 *
174 * Reverse transformation; *rev()
175 * -----------------------------
176 * Compute native spherical coordinates (phi,theta) from (x,y) coordinates in
177 * the plane of projection.
178 *
179 * Given:
180 * x,y const double
181 * Projected coordinates.
182 *
183 * Given and returned:
184 * prj AstPrjPrm* Projection parameters (see below).
185 *
186 * Returned:
187 * phi, double* Longitude and latitude of the projected point in
188 * theta native spherical coordinates, in degrees.
189 *
190 * Function return value:
191 * int Error status
192 * 0: Success.
193 * 1: Invalid projection parameters.
194 * 2: Invalid value of (x,y).
195 * 1: Invalid projection parameters.
196 *
197 * Projection parameters
198 * ---------------------
199 * The AstPrjPrm struct consists of the following:
200 *
201 * int flag
202 * This flag must be set to zero whenever any of p[] or r0 are set
203 * or changed. This signals the initialization routine to recompute
204 * intermediaries. flag may also be set to -1 to disable strict bounds
205 * checking for the AZP, SZP, TAN, SIN, ZPN, and COP projections.
206 *
207 * double r0
208 * r0; The radius of the generating sphere for the projection, a linear
209 * scaling parameter. If this is zero, it will be reset to the default
210 * value of 180/pi (the value for FITS WCS).
211 *
212 * double p[]
213 * Contains the projection parameters associated with the
214 * longitude axis.
215 *
216 * The remaining members of the AstPrjPrm struct are maintained by the
217 * initialization routines and should not be modified. This is done for the
218 * sake of efficiency and to allow an arbitrary number of contexts to be
219 * maintained simultaneously.
220 *
221 * char code[4]
222 * Three-letter projection code.
223 *
224 * double phi0, theta0
225 * Native longitude and latitude of the reference point, in degrees.
226 *
227 * double w[10]
228 * int n
229 * Intermediate values derived from the projection parameters.
230 *
231 * int (*astPRJfwd)()
232 * int (*astPRJrev)()
233 * Pointers to the forward and reverse projection routines.
234 *
235 * Usage of the p[] array as it applies to each projection is described in
236 * the prologue to each trio of projection routines.
237 *
238 * Argument checking
239 * -----------------
240 * Forward routines:
241 *
242 * The values of phi and theta (the native longitude and latitude)
243 * normally lie in the range [-180,180] for phi, and [-90,90] for theta.
244 * However, all forward projections will accept any value of phi and will
245 * not normalize it.
246 *
247 * The forward projection routines do not explicitly check that theta lies
248 * within the range [-90,90]. They do check for any value of theta which
249 * produces an invalid argument to the projection equations (e.g. leading
250 * to division by zero). The forward routines for AZP, SZP, TAN, SIN,
251 * ZPN, and COP also return error 2 if (phi,theta) corresponds to the
252 * overlapped (far) side of the projection but also return the
253 * corresponding value of (x,y). This strict bounds checking may be
254 * relaxed by setting prj->flag to -1 (rather than 0) when these
255 * projections are initialized.
256 *
257 * Reverse routines:
258 *
259 * Error checking on the projected coordinates (x,y) is limited to that
260 * required to ascertain whether a solution exists. Where a solution does
261 * exist no check is made that the value of phi and theta obtained lie
262 * within the ranges [-180,180] for phi, and [-90,90] for theta.
263 *
264 * Accuracy
265 * --------
266 * Closure to a precision of at least 1E-10 degree of longitude and latitude
267 * has been verified for typical projection parameters on the 1 degree grid
268 * of native longitude and latitude (to within 5 degrees of any latitude
269 * where the projection may diverge).
270 *
271 * Author: Mark Calabretta, Australia Telescope National Facility
272 * $Id$
273 *===========================================================================*/
274
275 /* Set the name of the module we are implementing. This indicates to
276 the header files that define class interfaces that they should make
277 "protected" symbols available. NB, this module is not a proper AST
278 class, but it defines this macro sanyway in order to get the protected
279 symbols defined in memory.h */
280
281 #include <math.h>
282 #include <string.h>
283 #include <stdlib.h>
284 #include "wcsmath.h"
285 #include "wcstrig.h"
286 #include "memory.h"
287 #include "proj.h"
288
289 /* Following variables are not needed in AST and are commented out to
290 avoid name clashes with other software systems (e.g. SkyCat) which
291 defines them.
292
293 int npcode = 28;
294 char pcodes[28][4] =
295 {"AZP", "SZP", "TAN", "STG", "SIN", "ARC", "ZPN", "ZEA", "AIR", "CYP",
296 "CEA", "CAR", "MER", "COP", "COE", "COD", "COO", "SFL", "PAR", "MOL",
297 "AIT", "BON", "PCO", "TSC", "CSC", "QSC", "HPX", "XPH"};
298 */
299
300 const int WCS__AZP = 101;
301 const int WCS__SZP = 102;
302 const int WCS__TAN = 103;
303 const int WCS__STG = 104;
304 const int WCS__SIN = 105;
305 const int WCS__ARC = 106;
306 const int WCS__ZPN = 107;
307 const int WCS__ZEA = 108;
308 const int WCS__AIR = 109;
309 const int WCS__CYP = 201;
310 const int WCS__CEA = 202;
311 const int WCS__CAR = 203;
312 const int WCS__MER = 204;
313 const int WCS__SFL = 301;
314 const int WCS__PAR = 302;
315 const int WCS__MOL = 303;
316 const int WCS__AIT = 401;
317 const int WCS__COP = 501;
318 const int WCS__COE = 502;
319 const int WCS__COD = 503;
320 const int WCS__COO = 504;
321 const int WCS__BON = 601;
322 const int WCS__PCO = 602;
323 const int WCS__TSC = 701;
324 const int WCS__CSC = 702;
325 const int WCS__QSC = 703;
326 const int WCS__HPX = 801;
327 const int WCS__XPH = 802;
328
329 /* Map error number to error message for each function. */
330 const char *astPRJset_errmsg[] = {
331 0,
332 "Invalid projection parameters"};
333
334 const char *astPRJfwd_errmsg[] = {
335 0,
336 "Invalid projection parameters",
337 "Invalid value of (phi,theta)"};
338
339 const char *astPRJrev_errmsg[] = {
340 0,
341 "Invalid projection parameters",
342 "Invalid value of (x,y)"};
343
344
345 #define copysign(X, Y) ((Y) < 0.0 ? -fabs(X) : fabs(X))
346
347
348
349 /*==========================================================================*/
350
astPRJset(pcode,prj)351 int astPRJset(pcode, prj)
352
353 const char pcode[4];
354 struct AstPrjPrm *prj;
355
356 {
357 /* Set pointers to the forward and reverse projection routines. */
358 if (strcmp(pcode, "AZP") == 0) {
359 astAZPset(prj);
360 } else if (strcmp(pcode, "SZP") == 0) {
361 astSZPset(prj);
362 } else if (strcmp(pcode, "TAN") == 0) {
363 astTANset(prj);
364 } else if (strcmp(pcode, "STG") == 0) {
365 astSTGset(prj);
366 } else if (strcmp(pcode, "SIN") == 0) {
367 astSINset(prj);
368 } else if (strcmp(pcode, "ARC") == 0) {
369 astARCset(prj);
370 } else if (strcmp(pcode, "ZPN") == 0) {
371 astZPNset(prj);
372 } else if (strcmp(pcode, "ZEA") == 0) {
373 astZEAset(prj);
374 } else if (strcmp(pcode, "AIR") == 0) {
375 astAIRset(prj);
376 } else if (strcmp(pcode, "CYP") == 0) {
377 astCYPset(prj);
378 } else if (strcmp(pcode, "CEA") == 0) {
379 astCEAset(prj);
380 } else if (strcmp(pcode, "CAR") == 0) {
381 astCARset(prj);
382 } else if (strcmp(pcode, "MER") == 0) {
383 astMERset(prj);
384 } else if (strcmp(pcode, "SFL") == 0) {
385 astSFLset(prj);
386 } else if (strcmp(pcode, "PAR") == 0) {
387 astPARset(prj);
388 } else if (strcmp(pcode, "MOL") == 0) {
389 astMOLset(prj);
390 } else if (strcmp(pcode, "AIT") == 0) {
391 astAITset(prj);
392 } else if (strcmp(pcode, "COP") == 0) {
393 astCOPset(prj);
394 } else if (strcmp(pcode, "COE") == 0) {
395 astCOEset(prj);
396 } else if (strcmp(pcode, "COD") == 0) {
397 astCODset(prj);
398 } else if (strcmp(pcode, "COO") == 0) {
399 astCOOset(prj);
400 } else if (strcmp(pcode, "BON") == 0) {
401 astBONset(prj);
402 } else if (strcmp(pcode, "PCO") == 0) {
403 astPCOset(prj);
404 } else if (strcmp(pcode, "TSC") == 0) {
405 astTSCset(prj);
406 } else if (strcmp(pcode, "CSC") == 0) {
407 astCSCset(prj);
408 } else if (strcmp(pcode, "QSC") == 0) {
409 astQSCset(prj);
410 } else if (strcmp(pcode, "HPX") == 0) {
411 astHPXset(prj);
412 } else if (strcmp(pcode, "XPH") == 0) {
413 astXPHset(prj);
414 } else {
415 /* Unrecognized projection code. */
416 return 1;
417 }
418
419 return 0;
420 }
421
422 /*--------------------------------------------------------------------------*/
423
astPRJfwd(phi,theta,prj,x,y)424 int astPRJfwd(phi, theta, prj, x, y)
425
426 const double phi, theta;
427 struct AstPrjPrm *prj;
428 double *x, *y;
429
430 {
431 return prj->astPRJfwd(phi, theta, prj, x, y);
432 }
433
434 /*--------------------------------------------------------------------------*/
435
astPRJrev(x,y,prj,phi,theta)436 int astPRJrev(x, y, prj, phi, theta)
437
438 const double x, y;
439 struct AstPrjPrm *prj;
440 double *phi, *theta;
441
442 {
443 return prj->astPRJrev(x, y, prj, phi, theta);
444 }
445
446 /*============================================================================
447 * AZP: zenithal/azimuthal perspective projection.
448 *
449 * Given:
450 * prj->p[1] Distance parameter, mu in units of r0.
451 * prj->p[2] Tilt angle, gamma in degrees.
452 *
453 * Given and/or returned:
454 * prj->flag AZP, or -AZP if prj->flag is given < 0.
455 * prj->r0 r0; reset to 180/pi if 0.
456 *
457 * Returned:
458 * prj->code "AZP"
459 * prj->phi0 0.0
460 * prj->theta0 90.0
461 * prj->w[0] r0*(mu+1)
462 * prj->w[1] tan(gamma)
463 * prj->w[2] sec(gamma)
464 * prj->w[3] cos(gamma)
465 * prj->w[4] sin(gamma)
466 * prj->w[5] asin(-1/mu) for |mu| >= 1, -90 otherwise
467 * prj->w[6] mu*cos(gamma)
468 * prj->w[7] 1 if |mu*cos(gamma)| < 1, 0 otherwise
469 * prj->astPRJfwd Pointer to astAZPfwd().
470 * prj->astPRJrev Pointer to astAZPrev().
471 *===========================================================================*/
472
astAZPset(prj)473 int astAZPset(prj)
474
475 struct AstPrjPrm *prj;
476
477 {
478 strcpy(prj->code, "AZP");
479 prj->flag = copysign(WCS__AZP, prj->flag);
480 prj->phi0 = 0.0;
481 prj->theta0 = 90.0;
482
483 if (prj->r0 == 0.0) prj->r0 = R2D;
484
485 prj->w[0] = prj->r0*(prj->p[1] + 1.0);
486 if (prj->w[0] == 0.0) {
487 return 1;
488 }
489
490 prj->w[3] = astCosd(prj->p[2]);
491 if (prj->w[3] == 0.0) {
492 return 1;
493 }
494
495 prj->w[2] = 1.0/prj->w[3];
496 prj->w[4] = astSind(prj->p[2]);
497 prj->w[1] = prj->w[4] / prj->w[3];
498
499 if (fabs(prj->p[1]) > 1.0) {
500 prj->w[5] = astASind(-1.0/prj->p[1]);
501 } else {
502 prj->w[5] = -90.0;
503 }
504
505 prj->w[6] = prj->p[1] * prj->w[3];
506 prj->w[7] = (fabs(prj->w[6]) < 1.0) ? 1.0 : 0.0;
507
508 prj->astPRJfwd = astAZPfwd;
509 prj->astPRJrev = astAZPrev;
510
511 return 0;
512 }
513
514 /*--------------------------------------------------------------------------*/
515
astAZPfwd(phi,theta,prj,x,y)516 int astAZPfwd(phi, theta, prj, x, y)
517
518 const double phi, theta;
519 struct AstPrjPrm *prj;
520 double *x, *y;
521
522 {
523 double a, b, cphi, cthe, r, s, t;
524
525 if (abs(prj->flag) != WCS__AZP) {
526 if (astAZPset(prj)) return 1;
527 }
528
529 cphi = astCosd(phi);
530 cthe = astCosd(theta);
531
532 s = prj->w[1]*cphi;
533 t = (prj->p[1] + astSind(theta)) + cthe*s;
534 if (t == 0.0) {
535 return 2;
536 }
537
538 r = prj->w[0]*cthe/t;
539 *x = r*astSind(phi);
540 *y = -r*cphi*prj->w[2];
541
542 /* Bounds checking. */
543 if (prj->flag > 0) {
544 /* Overlap. */
545 if (theta < prj->w[5]) {
546 return 2;
547 }
548
549 /* Divergence. */
550 if (prj->w[7] > 0.0) {
551 t = prj->p[1] / sqrt(1.0 + s*s);
552
553 if (fabs(t) <= 1.0) {
554 s = astATand(-s);
555 t = astASind(t);
556 a = s - t;
557 b = s + t + 180.0;
558
559 if (a > 90.0) a -= 360.0;
560 if (b > 90.0) b -= 360.0;
561
562 if (theta < ((a > b) ? a : b)) {
563 return 2;
564 }
565 }
566 }
567 }
568
569 return 0;
570 }
571
572 /*--------------------------------------------------------------------------*/
573
astAZPrev(x,y,prj,phi,theta)574 int astAZPrev(x, y, prj, phi, theta)
575
576 const double x, y;
577 struct AstPrjPrm *prj;
578 double *phi, *theta;
579
580 {
581 double a, b, r, s, t, ycosg;
582 const double tol = 1.0e-13;
583
584 if (abs(prj->flag) != WCS__AZP) {
585 if (astAZPset(prj)) return 1;
586 }
587
588 ycosg = y*prj->w[3];
589
590 r = sqrt(x*x + ycosg*ycosg);
591 if (r == 0.0) {
592 *phi = 0.0;
593 *theta = 90.0;
594 } else {
595 *phi = astATan2d(x, -ycosg);
596
597 s = r / (prj->w[0] + y*prj->w[4]);
598 t = s*prj->p[1]/sqrt(s*s + 1.0);
599
600 s = astATan2d(1.0, s);
601
602 if (fabs(t) > 1.0) {
603 t = copysign(90.0,t);
604 if (fabs(t) > 1.0+tol) {
605 return 2;
606 }
607 } else {
608 t = astASind(t);
609 }
610
611 a = s - t;
612 b = s + t + 180.0;
613
614 if (a > 90.0) a -= 360.0;
615 if (b > 90.0) b -= 360.0;
616
617 *theta = (a > b) ? a : b;
618 }
619
620 return 0;
621 }
622
623 /*============================================================================
624 * SZP: slant zenithal perspective projection.
625 *
626 * Given:
627 * prj->p[1] Distance of the point of projection from the centre of the
628 * generating sphere, mu in units of r0.
629 * prj->p[2] Native longitude, phi_c, and ...
630 * prj->p[3] Native latitude, theta_c, on the planewards side of the
631 * intersection of the line through the point of projection
632 * and the centre of the generating sphere, phi_c in degrees.
633 *
634 * Given and/or returned:
635 * prj->flag SZP, or -SZP if prj->flag is given < 0.
636 * prj->r0 r0; reset to 180/pi if 0.
637 *
638 * Returned:
639 * prj->code "SZP"
640 * prj->phi0 0.0
641 * prj->theta0 90.0
642 * prj->w[0] 1/r0
643 * prj->w[1] xp = -mu*cos(theta_c)*sin(phi_c)
644 * prj->w[2] yp = mu*cos(theta_c)*cos(phi_c)
645 * prj->w[3] zp = mu*sin(theta_c) + 1
646 * prj->w[4] r0*xp
647 * prj->w[5] r0*yp
648 * prj->w[6] r0*zp
649 * prj->w[7] (zp - 1)^2
650 * prj->w[8] asin(1-zp) if |1 - zp| < 1, -90 otherwise
651 * prj->astPRJfwd Pointer to astSZPfwd().
652 * prj->astPRJrev Pointer to astSZPrev().
653 *===========================================================================*/
654
astSZPset(prj)655 int astSZPset(prj)
656
657 struct AstPrjPrm *prj;
658
659 {
660 strcpy(prj->code, "SZP");
661 prj->flag = copysign(WCS__SZP, prj->flag);
662 prj->phi0 = 0.0;
663 prj->theta0 = 90.0;
664
665 if (prj->r0 == 0.0) prj->r0 = R2D;
666
667 prj->w[0] = 1.0/prj->r0;
668
669 prj->w[3] = prj->p[1] * astSind(prj->p[3]) + 1.0;
670 if (prj->w[3] == 0.0) {
671 return 1;
672 }
673
674 prj->w[1] = -prj->p[1] * astCosd(prj->p[3]) * astSind(prj->p[2]);
675 prj->w[2] = prj->p[1] * astCosd(prj->p[3]) * astCosd(prj->p[2]);
676 prj->w[4] = prj->r0 * prj->w[1];
677 prj->w[5] = prj->r0 * prj->w[2];
678 prj->w[6] = prj->r0 * prj->w[3];
679 prj->w[7] = (prj->w[3] - 1.0) * prj->w[3] - 1.0;
680
681 if (fabs(prj->w[3] - 1.0) < 1.0) {
682 prj->w[8] = astASind(1.0 - prj->w[3]);
683 } else {
684 prj->w[8] = -90.0;
685 }
686
687 prj->astPRJfwd = astSZPfwd;
688 prj->astPRJrev = astSZPrev;
689
690 return 0;
691 }
692
693 /*--------------------------------------------------------------------------*/
694
astSZPfwd(phi,theta,prj,x,y)695 int astSZPfwd(phi, theta, prj, x, y)
696
697 const double phi, theta;
698 struct AstPrjPrm *prj;
699 double *x, *y;
700
701 {
702 double a, b, cphi, cthe, s, sphi, t;
703
704 if (abs(prj->flag) != WCS__SZP) {
705 if (astSZPset(prj)) return 1;
706 }
707
708 cphi = astCosd(phi);
709 sphi = astSind(phi);
710 cthe = astCosd(theta);
711 s = 1.0 - astSind(theta);
712
713 t = prj->w[3] - s;
714 if (t == 0.0) {
715 return 2;
716 }
717
718 *x = (prj->w[6]*cthe*sphi - prj->w[4]*s)/t;
719 *y = -(prj->w[6]*cthe*cphi + prj->w[5]*s)/t;
720
721 /* Bounds checking. */
722 if (prj->flag > 0) {
723 /* Divergence. */
724 if (theta < prj->w[8]) {
725 return 2;
726 }
727
728 /* Overlap. */
729 if (fabs(prj->p[1]) > 1.0) {
730 s = prj->w[1]*sphi - prj->w[2]*cphi;
731 t = 1.0/sqrt(prj->w[7] + s*s);
732
733 if (fabs(t) <= 1.0) {
734 s = astATan2d(s, prj->w[3] - 1.0);
735 t = astASind(t);
736 a = s - t;
737 b = s + t + 180.0;
738
739 if (a > 90.0) a -= 360.0;
740 if (b > 90.0) b -= 360.0;
741
742 if (theta < ((a > b) ? a : b)) {
743 return 2;
744 }
745 }
746 }
747 }
748
749 return 0;
750 }
751
752 /*--------------------------------------------------------------------------*/
753
astSZPrev(x,y,prj,phi,theta)754 int astSZPrev(x, y, prj, phi, theta)
755
756 const double x, y;
757 struct AstPrjPrm *prj;
758 double *phi, *theta;
759
760 {
761 double a, b, c, d, r2, sth1, sth2, sthe, sxy, t, x1, xp, y1, yp, z;
762 const double tol = 1.0e-13;
763
764 if (abs(prj->flag) != WCS__SZP) {
765 if (astSZPset(prj)) return 1;
766 }
767
768 xp = x*prj->w[0];
769 yp = y*prj->w[0];
770 r2 = xp*xp + yp*yp;
771
772 x1 = (xp - prj->w[1])/prj->w[3];
773 y1 = (yp - prj->w[2])/prj->w[3];
774 sxy = xp*x1 + yp*y1;
775
776 if (r2 < 1.0e-10) {
777 /* Use small angle formula. */
778 z = r2/2.0;
779 *theta = 90.0 - R2D*sqrt(r2/(1.0 + sxy));
780
781 } else {
782 t = x1*x1 + y1*y1;
783 a = t + 1.0;
784 b = sxy - t;
785 c = r2 - sxy - sxy + t - 1.0;
786 d = b*b - a*c;
787
788 /* Check for a solution. */
789 if (d < 0.0) {
790 return 2;
791 }
792 d = sqrt(d);
793
794 /* Choose solution closest to pole. */
795 sth1 = (-b + d)/a;
796 sth2 = (-b - d)/a;
797 sthe = (sth1 > sth2) ? sth1 : sth2;
798 if (sthe > 1.0) {
799 if (sthe-1.0 < tol) {
800 sthe = 1.0;
801 } else {
802 sthe = (sth1 < sth2) ? sth1 : sth2;
803 }
804 }
805
806 if (sthe < -1.0) {
807 if (sthe+1.0 > -tol) {
808 sthe = -1.0;
809 }
810 }
811
812 if (sthe > 1.0 || sthe < -1.0) {
813 return 2;
814 }
815
816 *theta = astASind(sthe);
817
818 z = 1.0 - sthe;
819 }
820
821 *phi = astATan2d(xp - x1*z, -(yp - y1*z));
822
823 return 0;
824 }
825
826 /*============================================================================
827 * TAN: gnomonic projection.
828 *
829 * Given and/or returned:
830 * prj->flag TAN, or -TAN if prj->flag is given < 0.
831 * prj->r0 r0; reset to 180/pi if 0.
832 *
833 * Returned:
834 * prj->code "TAN"
835 * prj->phi0 0.0
836 * prj->theta0 90.0
837 * prj->astPRJfwd Pointer to astTANfwd().
838 * prj->astPRJrev Pointer to astTANrev().
839 *===========================================================================*/
840
astTANset(prj)841 int astTANset(prj)
842
843 struct AstPrjPrm *prj;
844
845 {
846 strcpy(prj->code, "TAN");
847 prj->flag = copysign(WCS__TAN, prj->flag);
848 prj->phi0 = 0.0;
849 prj->theta0 = 90.0;
850
851 if (prj->r0 == 0.0) prj->r0 = R2D;
852
853 prj->astPRJfwd = astTANfwd;
854 prj->astPRJrev = astTANrev;
855
856 return 0;
857 }
858
859 /*--------------------------------------------------------------------------*/
860
astTANfwd(phi,theta,prj,x,y)861 int astTANfwd(phi, theta, prj, x, y)
862
863 const double phi, theta;
864 struct AstPrjPrm *prj;
865 double *x, *y;
866
867 {
868 double r, s;
869
870 if (abs(prj->flag) != WCS__TAN) {
871 if(astTANset(prj)) return 1;
872 }
873
874 s = astSind(theta);
875 if (s == 0.0) {
876 return 2;
877 }
878
879 r = prj->r0*astCosd(theta)/s;
880 *x = r*astSind(phi);
881 *y = -r*astCosd(phi);
882
883 if (prj->flag > 0 && s < 0.0) {
884 return 2;
885 }
886
887 return 0;
888 }
889
890 /*--------------------------------------------------------------------------*/
891
astTANrev(x,y,prj,phi,theta)892 int astTANrev(x, y, prj, phi, theta)
893
894 const double x, y;
895 struct AstPrjPrm *prj;
896 double *phi, *theta;
897
898 {
899 double r;
900
901 if (abs(prj->flag) != WCS__TAN) {
902 if (astTANset(prj)) return 1;
903 }
904
905 r = sqrt(x*x + y*y);
906 if (r == 0.0) {
907 *phi = 0.0;
908 } else {
909 *phi = astATan2d(x, -y);
910 }
911 *theta = astATan2d(prj->r0, r);
912
913 return 0;
914 }
915
916 /*============================================================================
917 * STG: stereographic projection.
918 *
919 * Given and/or returned:
920 * prj->r0 r0; reset to 180/pi if 0.
921 *
922 * Returned:
923 * prj->code "STG"
924 * prj->flag STG
925 * prj->phi0 0.0
926 * prj->theta0 90.0
927 * prj->w[0] 2*r0
928 * prj->w[1] 1/(2*r0)
929 * prj->astPRJfwd Pointer to astSTGfwd().
930 * prj->astPRJrev Pointer to astSTGrev().
931 *===========================================================================*/
932
astSTGset(prj)933 int astSTGset(prj)
934
935 struct AstPrjPrm *prj;
936
937 {
938 strcpy(prj->code, "STG");
939 prj->flag = WCS__STG;
940 prj->phi0 = 0.0;
941 prj->theta0 = 90.0;
942
943 if (prj->r0 == 0.0) {
944 prj->r0 = R2D;
945 prj->w[0] = 360.0/PI;
946 prj->w[1] = PI/360.0;
947 } else {
948 prj->w[0] = 2.0*prj->r0;
949 prj->w[1] = 1.0/prj->w[0];
950 }
951
952 prj->astPRJfwd = astSTGfwd;
953 prj->astPRJrev = astSTGrev;
954
955 return 0;
956 }
957
958 /*--------------------------------------------------------------------------*/
959
astSTGfwd(phi,theta,prj,x,y)960 int astSTGfwd(phi, theta, prj, x, y)
961
962 const double phi, theta;
963 struct AstPrjPrm *prj;
964 double *x, *y;
965
966 {
967 double r, s;
968
969 if (prj->flag != WCS__STG) {
970 if (astSTGset(prj)) return 1;
971 }
972
973 s = 1.0 + astSind(theta);
974 if (s == 0.0) {
975 return 2;
976 }
977
978 r = prj->w[0]*astCosd(theta)/s;
979 *x = r*astSind(phi);
980 *y = -r*astCosd(phi);
981
982 return 0;
983 }
984
985 /*--------------------------------------------------------------------------*/
986
astSTGrev(x,y,prj,phi,theta)987 int astSTGrev(x, y, prj, phi, theta)
988
989 const double x, y;
990 struct AstPrjPrm *prj;
991 double *phi, *theta;
992
993 {
994 double r;
995
996 if (prj->flag != WCS__STG) {
997 if (astSTGset(prj)) return 1;
998 }
999
1000 r = sqrt(x*x + y*y);
1001 if (r == 0.0) {
1002 *phi = 0.0;
1003 } else {
1004 *phi = astATan2d(x, -y);
1005 }
1006 *theta = 90.0 - 2.0*astATand(r*prj->w[1]);
1007
1008 return 0;
1009 }
1010
1011 /*============================================================================
1012 * SIN: orthographic/synthesis projection.
1013 *
1014 * Given:
1015 * prj->p[1:2] Obliqueness parameters, xi and eta.
1016 *
1017 * Given and/or returned:
1018 * prj->flag SIN, or -SIN if prj->flag is given < 0.
1019 * prj->r0 r0; reset to 180/pi if 0.
1020 *
1021 * Returned:
1022 * prj->code "SIN"
1023 * prj->phi0 0.0
1024 * prj->theta0 90.0
1025 * prj->w[0] 1/r0
1026 * prj->w[1] xi**2 + eta**2
1027 * prj->w[2] xi**2 + eta**2 + 1
1028 * prj->w[3] xi**2 + eta**2 - 1
1029 * prj->astPRJfwd Pointer to astSINfwd().
1030 * prj->astPRJrev Pointer to astSINrev().
1031 *===========================================================================*/
1032
astSINset(prj)1033 int astSINset(prj)
1034
1035 struct AstPrjPrm *prj;
1036
1037 {
1038 strcpy(prj->code, "SIN");
1039 prj->flag = copysign(WCS__SIN, prj->flag);
1040 prj->phi0 = 0.0;
1041 prj->theta0 = 90.0;
1042
1043 if (prj->r0 == 0.0) prj->r0 = R2D;
1044
1045 prj->w[0] = 1.0/prj->r0;
1046 prj->w[1] = prj->p[1]*prj->p[1] + prj->p[2]*prj->p[2];
1047 prj->w[2] = prj->w[1] + 1.0;
1048 prj->w[3] = prj->w[1] - 1.0;
1049
1050 prj->astPRJfwd = astSINfwd;
1051 prj->astPRJrev = astSINrev;
1052
1053 return 0;
1054 }
1055
1056 /*--------------------------------------------------------------------------*/
1057
astSINfwd(phi,theta,prj,x,y)1058 int astSINfwd(phi, theta, prj, x, y)
1059
1060 const double phi, theta;
1061 struct AstPrjPrm *prj;
1062 double *x, *y;
1063
1064 {
1065 double cphi, cthe, sphi, t, z;
1066
1067 if (abs(prj->flag) != WCS__SIN) {
1068 if (astSINset(prj)) return 1;
1069 }
1070
1071 t = (90.0 - fabs(theta))*D2R;
1072 if (t < 1.0e-5) {
1073 if (theta > 0.0) {
1074 z = t*t/2.0;
1075 } else {
1076 z = 2.0 - t*t/2.0;
1077 }
1078 cthe = t;
1079 } else {
1080 z = 1.0 - astSind(theta);
1081 cthe = astCosd(theta);
1082 }
1083
1084 cphi = astCosd(phi);
1085 sphi = astSind(phi);
1086 *x = prj->r0*(cthe*sphi + prj->p[1]*z);
1087 *y = -prj->r0*(cthe*cphi - prj->p[2]*z);
1088
1089 /* Validate this solution. */
1090 if (prj->flag > 0) {
1091 if (prj->w[1] == 0.0) {
1092 /* Orthographic projection. */
1093 if (theta < 0.0) {
1094 return 2;
1095 }
1096 } else {
1097 /* "Synthesis" projection. */
1098 t = -astATand(prj->p[1]*sphi - prj->p[2]*cphi);
1099 if (theta < t) {
1100 return 2;
1101 }
1102 }
1103 }
1104
1105 return 0;
1106 }
1107
1108 /*--------------------------------------------------------------------------*/
1109
astSINrev(x,y,prj,phi,theta)1110 int astSINrev (x, y, prj, phi, theta)
1111
1112 const double x, y;
1113 struct AstPrjPrm *prj;
1114 double *phi, *theta;
1115
1116 {
1117 const double tol = 1.0e-13;
1118 double a, b, c, d, r2, sth1, sth2, sthe, sxy, x0, x1, xp, y0, y1, yp, z;
1119
1120 if (abs(prj->flag) != WCS__SIN) {
1121 if (astSINset(prj)) return 1;
1122 }
1123
1124 /* Compute intermediaries. */
1125 x0 = x*prj->w[0];
1126 y0 = y*prj->w[0];
1127 r2 = x0*x0 + y0*y0;
1128
1129 if (prj->w[1] == 0.0) {
1130 /* Orthographic projection. */
1131 if (r2 != 0.0) {
1132 *phi = astATan2d(x0, -y0);
1133 } else {
1134 *phi = 0.0;
1135 }
1136
1137 if (r2 < 0.5) {
1138 *theta = astACosd(sqrt(r2));
1139 } else if (r2 <= 1.0) {
1140 *theta = astASind(sqrt(1.0 - r2));
1141 } else {
1142 return 2;
1143 }
1144
1145 } else {
1146 /* "Synthesis" projection. */
1147 x1 = prj->p[1];
1148 y1 = prj->p[2];
1149 sxy = x0*x1 + y0*y1;
1150
1151 if (r2 < 1.0e-10) {
1152 /* Use small angle formula. */
1153 z = r2/2.0;
1154 *theta = 90.0 - R2D*sqrt(r2/(1.0 + sxy));
1155
1156 } else {
1157 a = prj->w[2];
1158 b = sxy - prj->w[1];
1159 c = r2 - sxy - sxy + prj->w[3];
1160 d = b*b - a*c;
1161
1162 /* Check for a solution. */
1163 if (d < 0.0) {
1164 return 2;
1165 }
1166 d = sqrt(d);
1167
1168 /* Choose solution closest to pole. */
1169 sth1 = (-b + d)/a;
1170 sth2 = (-b - d)/a;
1171 sthe = (sth1 > sth2) ? sth1 : sth2;
1172 if (sthe > 1.0) {
1173 if (sthe-1.0 < tol) {
1174 sthe = 1.0;
1175 } else {
1176 sthe = (sth1 < sth2) ? sth1 : sth2;
1177 }
1178 }
1179
1180 if (sthe < -1.0) {
1181 if (sthe+1.0 > -tol) {
1182 sthe = -1.0;
1183 }
1184 }
1185
1186 if (sthe > 1.0 || sthe < -1.0) {
1187 return 2;
1188 }
1189
1190 *theta = astASind(sthe);
1191 z = 1.0 - sthe;
1192 }
1193
1194 xp = -y0 + prj->p[2]*z;
1195 yp = x0 - prj->p[1]*z;
1196 if (xp == 0.0 && yp == 0.0) {
1197 *phi = 0.0;
1198 } else {
1199 *phi = astATan2d(yp,xp);
1200 }
1201 }
1202
1203 return 0;
1204 }
1205
1206 /*============================================================================
1207 * ARC: zenithal/azimuthal equidistant projection.
1208 *
1209 * Given and/or returned:
1210 * prj->r0 r0; reset to 180/pi if 0.
1211 *
1212 * Returned:
1213 * prj->code "ARC"
1214 * prj->flag ARC
1215 * prj->phi0 0.0
1216 * prj->theta0 90.0
1217 * prj->w[0] r0*(pi/180)
1218 * prj->w[1] (180/pi)/r0
1219 * prj->astPRJfwd Pointer to astARCfwd().
1220 * prj->astPRJrev Pointer to astARCrev().
1221 *===========================================================================*/
1222
astARCset(prj)1223 int astARCset(prj)
1224
1225 struct AstPrjPrm *prj;
1226
1227 {
1228 strcpy(prj->code, "ARC");
1229 prj->flag = WCS__ARC;
1230 prj->phi0 = 0.0;
1231 prj->theta0 = 90.0;
1232
1233 if (prj->r0 == 0.0) {
1234 prj->r0 = R2D;
1235 prj->w[0] = 1.0;
1236 prj->w[1] = 1.0;
1237 } else {
1238 prj->w[0] = prj->r0*D2R;
1239 prj->w[1] = 1.0/prj->w[0];
1240 }
1241
1242 prj->astPRJfwd = astARCfwd;
1243 prj->astPRJrev = astARCrev;
1244
1245 return 0;
1246 }
1247
1248 /*--------------------------------------------------------------------------*/
1249
astARCfwd(phi,theta,prj,x,y)1250 int astARCfwd(phi, theta, prj, x, y)
1251
1252 const double phi, theta;
1253 struct AstPrjPrm *prj;
1254 double *x, *y;
1255
1256 {
1257 double r;
1258
1259 if (prj->flag != WCS__ARC) {
1260 if (astARCset(prj)) return 1;
1261 }
1262
1263 r = prj->w[0]*(90.0 - theta);
1264 *x = r*astSind(phi);
1265 *y = -r*astCosd(phi);
1266
1267 return 0;
1268 }
1269
1270 /*--------------------------------------------------------------------------*/
1271
astARCrev(x,y,prj,phi,theta)1272 int astARCrev(x, y, prj, phi, theta)
1273
1274 const double x, y;
1275 struct AstPrjPrm *prj;
1276 double *phi, *theta;
1277
1278 {
1279 double r;
1280
1281 if (prj->flag != WCS__ARC) {
1282 if (astARCset(prj)) return 1;
1283 }
1284
1285 r = sqrt(x*x + y*y);
1286 if (r == 0.0) {
1287 *phi = 0.0;
1288 } else {
1289 *phi = astATan2d(x, -y);
1290 }
1291 *theta = 90.0 - r*prj->w[1];
1292
1293 return 0;
1294 }
1295
1296 /*============================================================================
1297 * ZPN: zenithal/azimuthal polynomial projection.
1298 *
1299 * Given:
1300 * prj->p[0:WCSLIB_MXPAR-1] Polynomial coefficients.
1301 *
1302 * Given and/or returned:
1303 * prj->flag ZPN, or -ZPN if prj->flag is given < 0.
1304 * prj->r0 r0; reset to 180/pi if 0.
1305 *
1306 * Returned:
1307 * prj->code "ZPN"
1308 * prj->phi0 0.0
1309 * prj->theta0 90.0
1310 * prj->n Degree of the polynomial, N.
1311 * prj->w[0] Co-latitude of the first point of inflection (N > 2).
1312 * prj->w[1] Radius of the first point of inflection (N > 2).
1313 * prj->astPRJfwd Pointer to astZPNfwd().
1314 * prj->astPRJrev Pointer to astZPNrev().
1315 *===========================================================================*/
1316
astZPNset(prj)1317 int astZPNset(prj)
1318
1319 struct AstPrjPrm *prj;
1320
1321 {
1322 int i, j, k, plen;
1323 double d, d1, d2, r, zd, zd1, zd2;
1324 const double tol = 1.0e-13;
1325
1326 strcpy(prj->code, "ZPN");
1327 prj->flag = copysign(WCS__ZPN, prj->flag);
1328 prj->phi0 = 0.0;
1329 prj->theta0 = 90.0;
1330
1331 if (prj->r0 == 0.0) prj->r0 = R2D;
1332
1333 /* Find the highest non-zero coefficient. */
1334 plen = astSizeOf( prj->p )/sizeof( double );
1335 for (k = plen-1; k >= 0 && prj->p[k] == 0.0; k--);
1336 if (k < 0) return 1;
1337
1338 prj->n = k;
1339
1340 if (k >= 3) {
1341 /* Find the point of inflection closest to the pole. */
1342 zd1 = 0.0;
1343 d1 = prj->p[1];
1344 if (d1 <= 0.0) {
1345 return 1;
1346 }
1347
1348 /* Find the point where the derivative first goes negative. */
1349 for (i = 0; i < 180; i++) {
1350 zd2 = i*D2R;
1351 d2 = 0.0;
1352 for (j = k; j > 0; j--) {
1353 d2 = d2*zd2 + j*prj->p[j];
1354 }
1355
1356 if (d2 <= 0.0) break;
1357 zd1 = zd2;
1358 d1 = d2;
1359 }
1360
1361 if (i == 180) {
1362 /* No negative derivative -> no point of inflection. */
1363 zd = PI;
1364 } else {
1365 /* Find where the derivative is zero. */
1366 for (i = 1; i <= 10; i++) {
1367 zd = zd1 - d1*(zd2-zd1)/(d2-d1);
1368
1369 d = 0.0;
1370 for (j = k; j > 0; j--) {
1371 d = d*zd + j*prj->p[j];
1372 }
1373
1374 if (fabs(d) < tol) break;
1375
1376 if (d < 0.0) {
1377 zd2 = zd;
1378 d2 = d;
1379 } else {
1380 zd1 = zd;
1381 d1 = d;
1382 }
1383 }
1384 }
1385
1386 r = 0.0;
1387 for (j = k; j >= 0; j--) {
1388 r = r*zd + prj->p[j];
1389 }
1390 prj->w[0] = zd;
1391 prj->w[1] = r;
1392 }
1393
1394 prj->astPRJfwd = astZPNfwd;
1395 prj->astPRJrev = astZPNrev;
1396
1397 return 0;
1398 }
1399
1400 /*--------------------------------------------------------------------------*/
1401
astZPNfwd(phi,theta,prj,x,y)1402 int astZPNfwd(phi, theta, prj, x, y)
1403
1404 const double phi, theta;
1405 struct AstPrjPrm *prj;
1406 double *x, *y;
1407
1408 {
1409 int j;
1410 double r, s;
1411
1412 if (abs(prj->flag) != WCS__ZPN) {
1413 if (astZPNset(prj)) return 1;
1414 }
1415
1416 s = (90.0 - theta)*D2R;
1417
1418 r = 0.0;
1419 for (j = prj->n; j >= 0; j--) {
1420 r = r*s + prj->p[j];
1421 }
1422 r = prj->r0*r;
1423
1424 *x = r*astSind(phi);
1425 *y = -r*astCosd(phi);
1426
1427 if (prj->flag > 0 && s > prj->w[0] && prj->n > 2 ) {
1428 return 2;
1429 }
1430
1431 return 0;
1432 }
1433
1434 /*--------------------------------------------------------------------------*/
1435
astZPNrev(x,y,prj,phi,theta)1436 int astZPNrev(x, y, prj, phi, theta)
1437
1438 const double x, y;
1439 struct AstPrjPrm *prj;
1440 double *phi, *theta;
1441
1442 {
1443 int i, j, k;
1444 double a, b, c, d, lambda, r, r1, r2, rt, zd, zd1, zd2;
1445 const double tol = 1.0e-13;
1446
1447 if (abs(prj->flag) != WCS__ZPN) {
1448 if (astZPNset(prj)) return 1;
1449 }
1450
1451 k = prj->n;
1452
1453 r = sqrt(x*x + y*y)/prj->r0;
1454
1455 if (k < 1) {
1456 /* Constant - no solution. */
1457 return 1;
1458 } else if (k == 1) {
1459 /* Linear. */
1460 zd = (r - prj->p[0])/prj->p[1];
1461 } else if (k == 2) {
1462 /* Quadratic. */
1463 a = prj->p[2];
1464 b = prj->p[1];
1465 c = prj->p[0] - r;
1466
1467 d = b*b - 4.0*a*c;
1468 if (d < 0.0) {
1469 return 2;
1470 }
1471 d = sqrt(d);
1472
1473 /* Choose solution closest to pole. */
1474 zd1 = (-b + d)/(2.0*a);
1475 zd2 = (-b - d)/(2.0*a);
1476 zd = (zd1<zd2) ? zd1 : zd2;
1477 if (zd < -tol) zd = (zd1>zd2) ? zd1 : zd2;
1478 if (zd < 0.0) {
1479 if (zd < -tol) {
1480 return 2;
1481 }
1482 zd = 0.0;
1483 } else if (zd > PI) {
1484 if (zd > PI+tol) {
1485 return 2;
1486 }
1487 zd = PI;
1488 }
1489 } else {
1490 /* Higher order - solve iteratively. */
1491 zd1 = 0.0;
1492 r1 = prj->p[0];
1493 zd2 = prj->w[0];
1494 r2 = prj->w[1];
1495
1496 if (r < r1) {
1497 if (r < r1-tol) {
1498 return 2;
1499 }
1500 zd = zd1;
1501 } else if (r > r2) {
1502 if (r > r2+tol) {
1503 return 2;
1504 }
1505 zd = zd2;
1506 } else {
1507 /* Disect the interval. */
1508 for (j = 0; j < 100; j++) {
1509 lambda = (r2 - r)/(r2 - r1);
1510 if (lambda < 0.1) {
1511 lambda = 0.1;
1512 } else if (lambda > 0.9) {
1513 lambda = 0.9;
1514 }
1515
1516 zd = zd2 - lambda*(zd2 - zd1);
1517
1518 rt = 0.0;
1519 for (i = k; i >= 0; i--) {
1520 rt = (rt * zd) + prj->p[i];
1521 }
1522
1523 if (rt < r) {
1524 if (r-rt < tol) break;
1525 r1 = rt;
1526 zd1 = zd;
1527 } else {
1528 if (rt-r < tol) break;
1529 r2 = rt;
1530 zd2 = zd;
1531 }
1532
1533 if (fabs(zd2-zd1) < tol) break;
1534 }
1535 }
1536 }
1537
1538 if (r == 0.0) {
1539 *phi = 0.0;
1540 } else {
1541 *phi = astATan2d(x, -y);
1542 }
1543 *theta = 90.0 - zd*R2D;
1544
1545 return 0;
1546 }
1547
1548 /*============================================================================
1549 * ZEA: zenithal/azimuthal equal area projection.
1550 *
1551 * Given and/or returned:
1552 * prj->r0 r0; reset to 180/pi if 0.
1553 *
1554 * Returned:
1555 * prj->code "ZEA"
1556 * prj->flag ZEA
1557 * prj->phi0 0.0
1558 * prj->theta0 90.0
1559 * prj->w[0] 2*r0
1560 * prj->w[1] 1/(2*r0)
1561 * prj->astPRJfwd Pointer to astZEAfwd().
1562 * prj->astPRJrev Pointer to astZEArev().
1563 *===========================================================================*/
1564
astZEAset(prj)1565 int astZEAset(prj)
1566
1567 struct AstPrjPrm *prj;
1568
1569 {
1570 strcpy(prj->code, "ZEA");
1571 prj->flag = WCS__ZEA;
1572 prj->phi0 = 0.0;
1573 prj->theta0 = 90.0;
1574
1575 if (prj->r0 == 0.0) {
1576 prj->r0 = R2D;
1577 prj->w[0] = 360.0/PI;
1578 prj->w[1] = PI/360.0;
1579 } else {
1580 prj->w[0] = 2.0*prj->r0;
1581 prj->w[1] = 1.0/prj->w[0];
1582 }
1583
1584 prj->astPRJfwd = astZEAfwd;
1585 prj->astPRJrev = astZEArev;
1586
1587 return 0;
1588 }
1589
1590 /*--------------------------------------------------------------------------*/
1591
astZEAfwd(phi,theta,prj,x,y)1592 int astZEAfwd(phi, theta, prj, x, y)
1593
1594 const double phi, theta;
1595 struct AstPrjPrm *prj;
1596 double *x, *y;
1597
1598 {
1599 double r;
1600
1601 if (prj->flag != WCS__ZEA) {
1602 if (astZEAset(prj)) return 1;
1603 }
1604
1605 r = prj->w[0]*astSind((90.0 - theta)/2.0);
1606 *x = r*astSind(phi);
1607 *y = -r*astCosd(phi);
1608
1609 return 0;
1610 }
1611
1612 /*--------------------------------------------------------------------------*/
1613
astZEArev(x,y,prj,phi,theta)1614 int astZEArev(x, y, prj, phi, theta)
1615
1616 const double x, y;
1617 struct AstPrjPrm *prj;
1618 double *phi, *theta;
1619
1620 {
1621 double r, s;
1622 const double tol = 1.0e-12;
1623
1624 if (prj->flag != WCS__ZEA) {
1625 if (astZEAset(prj)) return 1;
1626 }
1627
1628 r = sqrt(x*x + y*y);
1629 if (r == 0.0) {
1630 *phi = 0.0;
1631 } else {
1632 *phi = astATan2d(x, -y);
1633 }
1634
1635 s = r*prj->w[1];
1636 if (fabs(s) > 1.0) {
1637 if (fabs(r - prj->w[0]) < tol) {
1638 *theta = -90.0;
1639 } else {
1640 return 2;
1641 }
1642 } else {
1643 *theta = 90.0 - 2.0*astASind(s);
1644 }
1645
1646 return 0;
1647 }
1648
1649 /*============================================================================
1650 * AIR: Airy's projection.
1651 *
1652 * Given:
1653 * prj->p[1] Latitude theta_b within which the error is minimized, in
1654 * degrees.
1655 *
1656 * Given and/or returned:
1657 * prj->r0 r0; reset to 180/pi if 0.
1658 *
1659 * Returned:
1660 * prj->code "AIR"
1661 * prj->flag AIR
1662 * prj->phi0 0.0
1663 * prj->theta0 90.0
1664 * prj->w[0] 2*r0
1665 * prj->w[1] ln(cos(xi_b))/tan(xi_b)**2, where xi_b = (90-theta_b)/2
1666 * prj->w[2] 1/2 - prj->w[1]
1667 * prj->w[3] 2*r0*prj->w[2]
1668 * prj->w[4] tol, cutoff for using small angle approximation, in
1669 * radians.
1670 * prj->w[5] prj->w[2]*tol
1671 * prj->w[6] (180/pi)/prj->w[2]
1672 * prj->astPRJfwd Pointer to astAIRfwd().
1673 * prj->astPRJrev Pointer to astAIRrev().
1674 *===========================================================================*/
1675
astAIRset(prj)1676 int astAIRset(prj)
1677
1678 struct AstPrjPrm *prj;
1679
1680 {
1681 const double tol = 1.0e-4;
1682 double cxi;
1683
1684 strcpy(prj->code, "AIR");
1685 prj->flag = WCS__AIR;
1686 prj->phi0 = 0.0;
1687 prj->theta0 = 90.0;
1688
1689 if (prj->r0 == 0.0) prj->r0 = R2D;
1690
1691 prj->w[0] = 2.0*prj->r0;
1692 if (prj->p[1] == 90.0) {
1693 prj->w[1] = -0.5;
1694 prj->w[2] = 1.0;
1695 } else if (prj->p[1] > -90.0) {
1696 cxi = astCosd((90.0 - prj->p[1])/2.0);
1697 prj->w[1] = log(cxi)*(cxi*cxi)/(1.0-cxi*cxi);
1698 prj->w[2] = 0.5 - prj->w[1];
1699 } else {
1700 return 1;
1701 }
1702
1703 prj->w[3] = prj->w[0] * prj->w[2];
1704 prj->w[4] = tol;
1705 prj->w[5] = prj->w[2]*tol;
1706 prj->w[6] = R2D/prj->w[2];
1707
1708 prj->astPRJfwd = astAIRfwd;
1709 prj->astPRJrev = astAIRrev;
1710
1711 return 0;
1712 }
1713
1714 /*--------------------------------------------------------------------------*/
1715
astAIRfwd(phi,theta,prj,x,y)1716 int astAIRfwd(phi, theta, prj, x, y)
1717
1718 const double phi, theta;
1719 struct AstPrjPrm *prj;
1720 double *x, *y;
1721
1722 {
1723 double cxi, r, txi, xi;
1724
1725 if (prj->flag != WCS__AIR) {
1726 if (astAIRset(prj)) return 1;
1727 }
1728
1729 if (theta == 90.0) {
1730 r = 0.0;
1731 } else if (theta > -90.0) {
1732 xi = D2R*(90.0 - theta)/2.0;
1733 if (xi < prj->w[4]) {
1734 r = xi*prj->w[3];
1735 } else {
1736 cxi = astCosd((90.0 - theta)/2.0);
1737 txi = sqrt(1.0-cxi*cxi)/cxi;
1738 r = -prj->w[0]*(log(cxi)/txi + prj->w[1]*txi);
1739 }
1740 } else {
1741 return 2;
1742 }
1743
1744 *x = r*astSind(phi);
1745 *y = -r*astCosd(phi);
1746
1747 return 0;
1748 }
1749
1750 /*--------------------------------------------------------------------------*/
1751
astAIRrev(x,y,prj,phi,theta)1752 int astAIRrev(x, y, prj, phi, theta)
1753
1754 const double x, y;
1755 struct AstPrjPrm *prj;
1756 double *phi, *theta;
1757
1758 {
1759 int j;
1760 double cxi, lambda, r, r1, r2, rt, txi, x1, x2, xi;
1761 const double tol = 1.0e-12;
1762
1763 if (prj->flag != WCS__AIR) {
1764 if (astAIRset(prj)) return 1;
1765 }
1766
1767 r = sqrt(x*x + y*y)/prj->w[0];
1768
1769 if (r == 0.0) {
1770 xi = 0.0;
1771 } else if (r < prj->w[5]) {
1772 xi = r*prj->w[6];
1773 } else {
1774 /* Find a solution interval. */
1775 x1 = 1.0;
1776 r1 = 0.0;
1777 for (j = 0; j < 30; j++) {
1778 x2 = x1/2.0;
1779 txi = sqrt(1.0-x2*x2)/x2;
1780 r2 = -(log(x2)/txi + prj->w[1]*txi);
1781
1782 if (r2 >= r) break;
1783 x1 = x2;
1784 r1 = r2;
1785 }
1786 if (j == 30) return 2;
1787
1788 for (j = 0; j < 100; j++) {
1789 /* Weighted division of the interval. */
1790 lambda = (r2-r)/(r2-r1);
1791 if (lambda < 0.1) {
1792 lambda = 0.1;
1793 } else if (lambda > 0.9) {
1794 lambda = 0.9;
1795 }
1796 cxi = x2 - lambda*(x2-x1);
1797
1798 txi = sqrt(1.0-cxi*cxi)/cxi;
1799 rt = -(log(cxi)/txi + prj->w[1]*txi);
1800
1801 if (rt < r) {
1802 if (r-rt < tol) break;
1803 r1 = rt;
1804 x1 = cxi;
1805 } else {
1806 if (rt-r < tol) break;
1807 r2 = rt;
1808 x2 = cxi;
1809 }
1810 }
1811 if (j == 100) return 2;
1812
1813 xi = astACosd(cxi);
1814 }
1815
1816 if (r == 0.0) {
1817 *phi = 0.0;
1818 } else {
1819 *phi = astATan2d(x, -y);
1820 }
1821 *theta = 90.0 - 2.0*xi;
1822
1823 return 0;
1824 }
1825
1826 /*============================================================================
1827 * CYP: cylindrical perspective projection.
1828 *
1829 * Given:
1830 * prj->p[1] Distance of point of projection from the centre of the
1831 * generating sphere, mu, in units of r0.
1832 * prj->p[2] Radius of the cylinder of projection, lambda, in units of
1833 * r0.
1834 *
1835 * Given and/or returned:
1836 * prj->r0 r0; reset to 180/pi if 0.
1837 *
1838 * Returned:
1839 * prj->code "CYP"
1840 * prj->flag CYP
1841 * prj->phi0 0.0
1842 * prj->theta0 0.0
1843 * prj->w[0] r0*lambda*(pi/180)
1844 * prj->w[1] (180/pi)/(r0*lambda)
1845 * prj->w[2] r0*(mu + lambda)
1846 * prj->w[3] 1/(r0*(mu + lambda))
1847 * prj->astPRJfwd Pointer to astCYPfwd().
1848 * prj->astPRJrev Pointer to astCYPrev().
1849 *===========================================================================*/
1850
astCYPset(prj)1851 int astCYPset(prj)
1852
1853 struct AstPrjPrm *prj;
1854
1855 {
1856 strcpy(prj->code, "CYP");
1857 prj->flag = WCS__CYP;
1858 prj->phi0 = 0.0;
1859 prj->theta0 = 0.0;
1860
1861 if (prj->r0 == 0.0) {
1862 prj->r0 = R2D;
1863
1864 prj->w[0] = prj->p[2];
1865 if (prj->w[0] == 0.0) {
1866 return 1;
1867 }
1868
1869 prj->w[1] = 1.0/prj->w[0];
1870
1871 prj->w[2] = R2D*(prj->p[1] + prj->p[2]);
1872 if (prj->w[2] == 0.0) {
1873 return 1;
1874 }
1875
1876 prj->w[3] = 1.0/prj->w[2];
1877 } else {
1878 prj->w[0] = prj->r0*prj->p[2]*D2R;
1879 if (prj->w[0] == 0.0) {
1880 return 1;
1881 }
1882
1883 prj->w[1] = 1.0/prj->w[0];
1884
1885 prj->w[2] = prj->r0*(prj->p[1] + prj->p[2]);
1886 if (prj->w[2] == 0.0) {
1887 return 1;
1888 }
1889
1890 prj->w[3] = 1.0/prj->w[2];
1891 }
1892
1893 prj->astPRJfwd = astCYPfwd;
1894 prj->astPRJrev = astCYPrev;
1895
1896 return 0;
1897 }
1898
1899 /*--------------------------------------------------------------------------*/
1900
astCYPfwd(phi,theta,prj,x,y)1901 int astCYPfwd(phi, theta, prj, x, y)
1902
1903 const double phi, theta;
1904 struct AstPrjPrm *prj;
1905 double *x, *y;
1906
1907 {
1908 double s;
1909
1910 if (prj->flag != WCS__CYP) {
1911 if (astCYPset(prj)) return 1;
1912 }
1913
1914 s = prj->p[1] + astCosd(theta);
1915 if (s == 0.0) {
1916 return 2;
1917 }
1918
1919 *x = prj->w[0]*phi;
1920 *y = prj->w[2]*astSind(theta)/s;
1921
1922 return 0;
1923 }
1924
1925 /*--------------------------------------------------------------------------*/
1926
astCYPrev(x,y,prj,phi,theta)1927 int astCYPrev(x, y, prj, phi, theta)
1928
1929 const double x, y;
1930 struct AstPrjPrm *prj;
1931 double *phi, *theta;
1932
1933 {
1934 double eta;
1935 double a;
1936 const double tol = 1.0e-13;
1937
1938 if (prj->flag != WCS__CYP) {
1939 if (astCYPset(prj)) return 1;
1940 }
1941
1942 *phi = x*prj->w[1];
1943 eta = y*prj->w[3];
1944
1945 a = eta*prj->p[1]/sqrt(eta*eta+1.0);
1946 if( fabs( a ) < 1.0 ) {
1947 *theta = astATan2d(eta,1.0) + astASind( a );
1948
1949 } else if( fabs( a ) < 1.0 + tol ) {
1950 if( a > 0.0 ){
1951 *theta = astATan2d(eta,1.0) + 90.0;
1952 } else {
1953 *theta = astATan2d(eta,1.0) - 90.0;
1954 }
1955
1956 } else {
1957 return 2;
1958 }
1959
1960 return 0;
1961 }
1962
1963 /*============================================================================
1964 * CEA: cylindrical equal area projection.
1965 *
1966 * Given:
1967 * prj->p[1] Square of the cosine of the latitude at which the
1968 * projection is conformal, lambda.
1969 *
1970 * Given and/or returned:
1971 * prj->r0 r0; reset to 180/pi if 0.
1972 *
1973 * Returned:
1974 * prj->code "CEA"
1975 * prj->flag CEA
1976 * prj->phi0 0.0
1977 * prj->theta0 0.0
1978 * prj->w[0] r0*(pi/180)
1979 * prj->w[1] (180/pi)/r0
1980 * prj->w[2] r0/lambda
1981 * prj->w[3] lambda/r0
1982 * prj->astPRJfwd Pointer to astCEAfwd().
1983 * prj->astPRJrev Pointer to astCEArev().
1984 *===========================================================================*/
1985
astCEAset(prj)1986 int astCEAset(prj)
1987
1988 struct AstPrjPrm *prj;
1989
1990 {
1991 strcpy(prj->code, "CEA");
1992 prj->flag = WCS__CEA;
1993 prj->phi0 = 0.0;
1994 prj->theta0 = 0.0;
1995
1996 if (prj->r0 == 0.0) {
1997 prj->r0 = R2D;
1998 prj->w[0] = 1.0;
1999 prj->w[1] = 1.0;
2000 if (prj->p[1] <= 0.0 || prj->p[1] > 1.0) {
2001 return 1;
2002 }
2003 prj->w[2] = prj->r0/prj->p[1];
2004 prj->w[3] = prj->p[1]/prj->r0;
2005 } else {
2006 prj->w[0] = prj->r0*D2R;
2007 prj->w[1] = R2D/prj->r0;
2008 if (prj->p[1] <= 0.0 || prj->p[1] > 1.0) {
2009 return 1;
2010 }
2011 prj->w[2] = prj->r0/prj->p[1];
2012 prj->w[3] = prj->p[1]/prj->r0;
2013 }
2014
2015 prj->astPRJfwd = astCEAfwd;
2016 prj->astPRJrev = astCEArev;
2017
2018 return 0;
2019 }
2020
2021 /*--------------------------------------------------------------------------*/
2022
astCEAfwd(phi,theta,prj,x,y)2023 int astCEAfwd(phi, theta, prj, x, y)
2024
2025 const double phi, theta;
2026 struct AstPrjPrm *prj;
2027 double *x, *y;
2028
2029 {
2030 if (prj->flag != WCS__CEA) {
2031 if (astCEAset(prj)) return 1;
2032 }
2033
2034 *x = prj->w[0]*phi;
2035 *y = prj->w[2]*astSind(theta);
2036
2037 return 0;
2038 }
2039
2040 /*--------------------------------------------------------------------------*/
2041
astCEArev(x,y,prj,phi,theta)2042 int astCEArev(x, y, prj, phi, theta)
2043
2044 const double x, y;
2045 struct AstPrjPrm *prj;
2046 double *phi, *theta;
2047
2048 {
2049 double s;
2050 const double tol = 1.0e-13;
2051
2052 if (prj->flag != WCS__CEA) {
2053 if (astCEAset(prj)) return 1;
2054 }
2055
2056 s = y*prj->w[3];
2057 if (fabs(s) > 1.0) {
2058 if (fabs(s) > 1.0+tol) {
2059 return 2;
2060 }
2061 s = copysign(1.0,s);
2062 }
2063
2064 *phi = x*prj->w[1];
2065 *theta = astASind(s);
2066
2067 return 0;
2068 }
2069
2070 /*============================================================================
2071 * CAR: Cartesian projection.
2072 *
2073 * Given and/or returned:
2074 * prj->r0 r0; reset to 180/pi if 0.
2075 *
2076 * Returned:
2077 * prj->code "CAR"
2078 * prj->flag CAR
2079 * prj->phi0 0.0
2080 * prj->theta0 0.0
2081 * prj->w[0] r0*(pi/180)
2082 * prj->w[1] (180/pi)/r0
2083 * prj->astPRJfwd Pointer to astCARfwd().
2084 * prj->astPRJrev Pointer to astCARrev().
2085 *===========================================================================*/
2086
astCARset(prj)2087 int astCARset(prj)
2088
2089 struct AstPrjPrm *prj;
2090
2091 {
2092 strcpy(prj->code, "CAR");
2093 prj->flag = WCS__CAR;
2094 prj->phi0 = 0.0;
2095 prj->theta0 = 0.0;
2096
2097 if (prj->r0 == 0.0) {
2098 prj->r0 = R2D;
2099 prj->w[0] = 1.0;
2100 prj->w[1] = 1.0;
2101 } else {
2102 prj->w[0] = prj->r0*D2R;
2103 prj->w[1] = 1.0/prj->w[0];
2104 }
2105
2106 prj->astPRJfwd = astCARfwd;
2107 prj->astPRJrev = astCARrev;
2108
2109 return 0;
2110 }
2111
2112 /*--------------------------------------------------------------------------*/
2113
astCARfwd(phi,theta,prj,x,y)2114 int astCARfwd(phi, theta, prj, x, y)
2115
2116 const double phi, theta;
2117 struct AstPrjPrm *prj;
2118 double *x, *y;
2119
2120 {
2121 if (prj->flag != WCS__CAR) {
2122 if (astCARset(prj)) return 1;
2123 }
2124
2125 *x = prj->w[0]*phi;
2126 *y = prj->w[0]*theta;
2127
2128 return 0;
2129 }
2130
2131 /*--------------------------------------------------------------------------*/
2132
astCARrev(x,y,prj,phi,theta)2133 int astCARrev(x, y, prj, phi, theta)
2134
2135 const double x, y;
2136 struct AstPrjPrm *prj;
2137 double *phi, *theta;
2138
2139 {
2140 if (prj->flag != WCS__CAR) {
2141 if (astCARset(prj)) return 1;
2142 }
2143
2144 *phi = prj->w[1]*x;
2145 *theta = prj->w[1]*y;
2146
2147 return 0;
2148 }
2149
2150 /*============================================================================
2151 * MER: Mercator's projection.
2152 *
2153 * Given and/or returned:
2154 * prj->r0 r0; reset to 180/pi if 0.
2155 *
2156 * Returned:
2157 * prj->code "MER"
2158 * prj->flag MER
2159 * prj->phi0 0.0
2160 * prj->theta0 0.0
2161 * prj->w[0] r0*(pi/180)
2162 * prj->w[1] (180/pi)/r0
2163 * prj->astPRJfwd Pointer to astMERfwd().
2164 * prj->astPRJrev Pointer to astMERrev().
2165 *===========================================================================*/
2166
astMERset(prj)2167 int astMERset(prj)
2168
2169 struct AstPrjPrm *prj;
2170
2171 {
2172 strcpy(prj->code, "MER");
2173 prj->flag = WCS__MER;
2174 prj->phi0 = 0.0;
2175 prj->theta0 = 0.0;
2176
2177 if (prj->r0 == 0.0) {
2178 prj->r0 = R2D;
2179 prj->w[0] = 1.0;
2180 prj->w[1] = 1.0;
2181 } else {
2182 prj->w[0] = prj->r0*D2R;
2183 prj->w[1] = 1.0/prj->w[0];
2184 }
2185
2186 prj->astPRJfwd = astMERfwd;
2187 prj->astPRJrev = astMERrev;
2188
2189 return 0;
2190 }
2191
2192 /*--------------------------------------------------------------------------*/
2193
astMERfwd(phi,theta,prj,x,y)2194 int astMERfwd(phi, theta, prj, x, y)
2195
2196 const double phi, theta;
2197 struct AstPrjPrm *prj;
2198 double *x, *y;
2199
2200 {
2201 if (prj->flag != WCS__MER) {
2202 if (astMERset(prj)) return 1;
2203 }
2204
2205 if (theta <= -90.0 || theta >= 90.0) {
2206 return 2;
2207 }
2208
2209 *x = prj->w[0]*phi;
2210 *y = prj->r0*log(astTand((90.0+theta)/2.0));
2211
2212 return 0;
2213 }
2214
2215 /*--------------------------------------------------------------------------*/
2216
astMERrev(x,y,prj,phi,theta)2217 int astMERrev(x, y, prj, phi, theta)
2218
2219 const double x, y;
2220 struct AstPrjPrm *prj;
2221 double *phi, *theta;
2222
2223 {
2224 if (prj->flag != WCS__MER) {
2225 if (astMERset(prj)) return 1;
2226 }
2227
2228 *phi = x*prj->w[1];
2229 *theta = 2.0*astATand(exp(y/prj->r0)) - 90.0;
2230
2231 return 0;
2232 }
2233
2234 /*============================================================================
2235 * SFL: Sanson-Flamsteed ("global sinusoid") projection.
2236 *
2237 * Given and/or returned:
2238 * prj->r0 r0; reset to 180/pi if 0.
2239 *
2240 * Returned:
2241 * prj->code "SFL"
2242 * prj->flag SFL
2243 * prj->phi0 0.0
2244 * prj->theta0 0.0
2245 * prj->w[0] r0*(pi/180)
2246 * prj->w[1] (180/pi)/r0
2247 * prj->astPRJfwd Pointer to astSFLfwd().
2248 * prj->astPRJrev Pointer to astSFLrev().
2249 *===========================================================================*/
2250
astSFLset(prj)2251 int astSFLset(prj)
2252
2253 struct AstPrjPrm *prj;
2254
2255 {
2256 strcpy(prj->code, "SFL");
2257 prj->flag = WCS__SFL;
2258 prj->phi0 = 0.0;
2259 prj->theta0 = 0.0;
2260
2261 if (prj->r0 == 0.0) {
2262 prj->r0 = R2D;
2263 prj->w[0] = 1.0;
2264 prj->w[1] = 1.0;
2265 } else {
2266 prj->w[0] = prj->r0*D2R;
2267 prj->w[1] = 1.0/prj->w[0];
2268 }
2269
2270 prj->astPRJfwd = astSFLfwd;
2271 prj->astPRJrev = astSFLrev;
2272
2273 return 0;
2274 }
2275
2276 /*--------------------------------------------------------------------------*/
2277
astSFLfwd(phi,theta,prj,x,y)2278 int astSFLfwd(phi, theta, prj, x, y)
2279
2280 const double phi, theta;
2281 struct AstPrjPrm *prj;
2282 double *x, *y;
2283
2284 {
2285 if (prj->flag != WCS__SFL) {
2286 if (astSFLset(prj)) return 1;
2287 }
2288
2289 *x = prj->w[0]*phi*astCosd(theta);
2290 *y = prj->w[0]*theta;
2291
2292 return 0;
2293 }
2294
2295 /*--------------------------------------------------------------------------*/
2296
astSFLrev(x,y,prj,phi,theta)2297 int astSFLrev(x, y, prj, phi, theta)
2298
2299 const double x, y;
2300 struct AstPrjPrm *prj;
2301 double *phi, *theta;
2302
2303 {
2304 double w;
2305
2306 if (prj->flag != WCS__SFL) {
2307 if (astSFLset(prj)) return 1;
2308 }
2309
2310 w = cos(y/prj->r0);
2311 if (w == 0.0) {
2312 *phi = 0.0;
2313 } else {
2314 *phi = x*prj->w[1]/cos(y/prj->r0);
2315 }
2316 *theta = y*prj->w[1];
2317
2318 return 0;
2319 }
2320
2321 /*============================================================================
2322 * PAR: parabolic projection.
2323 *
2324 * Given and/or returned:
2325 * prj->r0 r0; reset to 180/pi if 0.
2326 *
2327 * Returned:
2328 * prj->code "PAR"
2329 * prj->flag PAR
2330 * prj->phi0 0.0
2331 * prj->theta0 0.0
2332 * prj->w[0] r0*(pi/180)
2333 * prj->w[1] (180/pi)/r0
2334 * prj->w[2] pi*r0
2335 * prj->w[3] 1/(pi*r0)
2336 * prj->astPRJfwd Pointer to astPARfwd().
2337 * prj->astPRJrev Pointer to astPARrev().
2338 *===========================================================================*/
2339
astPARset(prj)2340 int astPARset(prj)
2341
2342 struct AstPrjPrm *prj;
2343
2344 {
2345 strcpy(prj->code, "PAR");
2346 prj->flag = WCS__PAR;
2347 prj->phi0 = 0.0;
2348 prj->theta0 = 0.0;
2349
2350 if (prj->r0 == 0.0) {
2351 prj->r0 = R2D;
2352 prj->w[0] = 1.0;
2353 prj->w[1] = 1.0;
2354 prj->w[2] = 180.0;
2355 prj->w[3] = 1.0/prj->w[2];
2356 } else {
2357 prj->w[0] = prj->r0*D2R;
2358 prj->w[1] = 1.0/prj->w[0];
2359 prj->w[2] = PI*prj->r0;
2360 prj->w[3] = 1.0/prj->w[2];
2361 }
2362
2363 prj->astPRJfwd = astPARfwd;
2364 prj->astPRJrev = astPARrev;
2365
2366 return 0;
2367 }
2368
2369 /*--------------------------------------------------------------------------*/
2370
astPARfwd(phi,theta,prj,x,y)2371 int astPARfwd(phi, theta, prj, x, y)
2372
2373 const double phi, theta;
2374 struct AstPrjPrm *prj;
2375 double *x, *y;
2376
2377 {
2378 double s;
2379
2380 if (prj->flag != WCS__PAR) {
2381 if (astPARset(prj)) return 1;
2382 }
2383
2384 s = astSind(theta/3.0);
2385 *x = prj->w[0]*phi*(1.0 - 4.0*s*s);
2386 *y = prj->w[2]*s;
2387
2388 return 0;
2389 }
2390
2391 /*--------------------------------------------------------------------------*/
2392
astPARrev(x,y,prj,phi,theta)2393 int astPARrev(x, y, prj, phi, theta)
2394
2395 const double x, y;
2396 struct AstPrjPrm *prj;
2397 double *phi, *theta;
2398
2399 {
2400 double s, t;
2401
2402 if (prj->flag != WCS__PAR) {
2403 if (astPARset(prj)) return 1;
2404 }
2405
2406 s = y*prj->w[3];
2407 if (s > 1.0 || s < -1.0) {
2408 return 2;
2409 }
2410
2411 t = 1.0 - 4.0*s*s;
2412 if (t == 0.0) {
2413 if (x == 0.0) {
2414 *phi = 0.0;
2415 } else {
2416 return 2;
2417 }
2418 } else {
2419 *phi = prj->w[1]*x/t;
2420 }
2421
2422 *theta = 3.0*astASind(s);
2423
2424 return 0;
2425 }
2426
2427 /*============================================================================
2428 * MOL: Mollweide's projection.
2429 *
2430 * Given and/or returned:
2431 * prj->r0 r0; reset to 180/pi if 0.
2432 *
2433 * Returned:
2434 * prj->code "MOL"
2435 * prj->flag MOL
2436 * prj->phi0 0.0
2437 * prj->theta0 0.0
2438 * prj->w[0] sqrt(2)*r0
2439 * prj->w[1] sqrt(2)*r0/90
2440 * prj->w[2] 1/(sqrt(2)*r0)
2441 * prj->w[3] 90/r0
2442 * prj->astPRJfwd Pointer to astMOLfwd().
2443 * prj->astPRJrev Pointer to astMOLrev().
2444 *===========================================================================*/
2445
astMOLset(prj)2446 int astMOLset(prj)
2447
2448 struct AstPrjPrm *prj;
2449
2450 {
2451 strcpy(prj->code, "MOL");
2452 prj->flag = WCS__MOL;
2453 prj->phi0 = 0.0;
2454 prj->theta0 = 0.0;
2455
2456 if (prj->r0 == 0.0) prj->r0 = R2D;
2457
2458 prj->w[0] = SQRT2*prj->r0;
2459 prj->w[1] = prj->w[0]/90.0;
2460 prj->w[2] = 1.0/prj->w[0];
2461 prj->w[3] = 90.0/prj->r0;
2462 prj->w[4] = 2.0/PI;
2463
2464 prj->astPRJfwd = astMOLfwd;
2465 prj->astPRJrev = astMOLrev;
2466
2467 return 0;
2468 }
2469
2470 /*--------------------------------------------------------------------------*/
2471
astMOLfwd(phi,theta,prj,x,y)2472 int astMOLfwd(phi, theta, prj, x, y)
2473
2474 const double phi, theta;
2475 struct AstPrjPrm *prj;
2476 double *x, *y;
2477
2478 {
2479 int j;
2480 double gamma, resid, u, v, v0, v1;
2481 const double tol = 1.0e-13;
2482
2483 if (prj->flag != WCS__MOL) {
2484 if (astMOLset(prj)) return 1;
2485 }
2486
2487 if (fabs(theta) == 90.0) {
2488 *x = 0.0;
2489 *y = copysign(prj->w[0],theta);
2490 } else if (theta == 0.0) {
2491 *x = prj->w[1]*phi;
2492 *y = 0.0;
2493 } else {
2494 u = PI*astSind(theta);
2495 v0 = -PI;
2496 v1 = PI;
2497 v = u;
2498 for (j = 0; j < 100; j++) {
2499 resid = (v - u) + sin(v);
2500 if (resid < 0.0) {
2501 if (resid > -tol) break;
2502 v0 = v;
2503 } else {
2504 if (resid < tol) break;
2505 v1 = v;
2506 }
2507 v = (v0 + v1)/2.0;
2508 }
2509
2510 gamma = v/2.0;
2511 *x = prj->w[1]*phi*cos(gamma);
2512 *y = prj->w[0]*sin(gamma);
2513 }
2514
2515 return 0;
2516 }
2517
2518 /*--------------------------------------------------------------------------*/
2519
astMOLrev(x,y,prj,phi,theta)2520 int astMOLrev(x, y, prj, phi, theta)
2521
2522 const double x, y;
2523 struct AstPrjPrm *prj;
2524 double *phi, *theta;
2525
2526 {
2527 double s, y0, z;
2528 const double tol = 1.0e-12;
2529
2530 if (prj->flag != WCS__MOL) {
2531 if (astMOLset(prj)) return 1;
2532 }
2533
2534 y0 = y/prj->r0;
2535 s = 2.0 - y0*y0;
2536 if (s <= tol) {
2537 if (s < -tol) {
2538 return 2;
2539 }
2540 s = 0.0;
2541
2542 if (fabs(x) > tol) {
2543 return 2;
2544 }
2545 *phi = 0.0;
2546 } else {
2547 s = sqrt(s);
2548 *phi = prj->w[3]*x/s;
2549 }
2550
2551 z = y*prj->w[2];
2552 if (fabs(z) > 1.0) {
2553 if (fabs(z) > 1.0+tol) {
2554 return 2;
2555 }
2556 z = copysign(1.0,z) + y0*s/PI;
2557 } else {
2558 z = asin(z)*prj->w[4] + y0*s/PI;
2559 }
2560
2561 if (fabs(z) > 1.0) {
2562 if (fabs(z) > 1.0+tol) {
2563 return 2;
2564 }
2565 z = copysign(1.0,z);
2566 }
2567
2568 *theta = astASind(z);
2569
2570 return 0;
2571 }
2572
2573 /*============================================================================
2574 * AIT: Hammer-Aitoff projection.
2575 *
2576 * Given and/or returned:
2577 * prj->r0 r0; reset to 180/pi if 0.
2578 *
2579 * Returned:
2580 * prj->code "AIT"
2581 * prj->flag AIT
2582 * prj->phi0 0.0
2583 * prj->theta0 0.0
2584 * prj->w[0] 2*r0**2
2585 * prj->w[1] 1/(2*r0)**2
2586 * prj->w[2] 1/(4*r0)**2
2587 * prj->w[3] 1/(2*r0)
2588 * prj->astPRJfwd Pointer to astAITfwd().
2589 * prj->astPRJrev Pointer to astAITrev().
2590 *===========================================================================*/
2591
astAITset(prj)2592 int astAITset(prj)
2593
2594 struct AstPrjPrm *prj;
2595
2596 {
2597 strcpy(prj->code, "AIT");
2598 prj->flag = WCS__AIT;
2599 prj->phi0 = 0.0;
2600 prj->theta0 = 0.0;
2601
2602 if (prj->r0 == 0.0) prj->r0 = R2D;
2603
2604 prj->w[0] = 2.0*prj->r0*prj->r0;
2605 prj->w[1] = 1.0/(2.0*prj->w[0]);
2606 prj->w[2] = prj->w[1]/4.0;
2607 prj->w[3] = 1.0/(2.0*prj->r0);
2608
2609 prj->astPRJfwd = astAITfwd;
2610 prj->astPRJrev = astAITrev;
2611
2612 return 0;
2613 }
2614
2615 /*--------------------------------------------------------------------------*/
2616
astAITfwd(phi,theta,prj,x,y)2617 int astAITfwd(phi, theta, prj, x, y)
2618
2619 const double phi, theta;
2620 struct AstPrjPrm *prj;
2621 double *x, *y;
2622
2623 {
2624 double cthe, w;
2625
2626 if (prj->flag != WCS__AIT) {
2627 if (astAITset(prj)) return 1;
2628 }
2629
2630 cthe = astCosd(theta);
2631 w = sqrt(prj->w[0]/(1.0 + cthe*astCosd(phi/2.0)));
2632 *x = 2.0*w*cthe*astSind(phi/2.0);
2633 *y = w*astSind(theta);
2634
2635 return 0;
2636 }
2637
2638 /*--------------------------------------------------------------------------*/
2639
astAITrev(x,y,prj,phi,theta)2640 int astAITrev(x, y, prj, phi, theta)
2641
2642 const double x, y;
2643 struct AstPrjPrm *prj;
2644 double *phi, *theta;
2645
2646 {
2647 double s, u, xp, yp, z;
2648 const double tol = 1.0e-13;
2649
2650 if (prj->flag != WCS__AIT) {
2651 if (astAITset(prj)) return 1;
2652 }
2653
2654 u = 1.0 - x*x*prj->w[2] - y*y*prj->w[1];
2655 if (u < 0.0) {
2656 if (u < -tol) {
2657 return 2;
2658 }
2659
2660 u = 0.0;
2661 }
2662
2663 z = sqrt(u);
2664 s = z*y/prj->r0;
2665 if (fabs(s) > 1.0) {
2666 if (fabs(s) > 1.0+tol) {
2667 return 2;
2668 }
2669 s = copysign(1.0,s);
2670 }
2671
2672 xp = 2.0*z*z - 1.0;
2673 yp = z*x*prj->w[3];
2674 if (xp == 0.0 && yp == 0.0) {
2675 *phi = 0.0;
2676 } else {
2677 *phi = 2.0*astATan2d(yp, xp);
2678 }
2679 *theta = astASind(s);
2680
2681 return 0;
2682 }
2683
2684 /*============================================================================
2685 * COP: conic perspective projection.
2686 *
2687 * Given:
2688 * prj->p[1] sigma = (theta2+theta1)/2
2689 * prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the
2690 * latitudes of the standard parallels, in degrees.
2691 *
2692 * Given and/or returned:
2693 * prj->flag COP, or -COP if prj->flag is given < 0.
2694 * prj->r0 r0; reset to 180/pi if 0.
2695 *
2696 * Returned:
2697 * prj->code "COP"
2698 * prj->phi0 0.0
2699 * prj->theta0 sigma
2700 * prj->w[0] C = sin(sigma)
2701 * prj->w[1] 1/C
2702 * prj->w[2] Y0 = r0*cos(delta)*cot(sigma)
2703 * prj->w[3] r0*cos(delta)
2704 * prj->w[4] 1/(r0*cos(delta)
2705 * prj->w[5] cot(sigma)
2706 * prj->astPRJfwd Pointer to astCOPfwd().
2707 * prj->astPRJrev Pointer to astCOPrev().
2708 *===========================================================================*/
2709
astCOPset(prj)2710 int astCOPset(prj)
2711
2712 struct AstPrjPrm *prj;
2713
2714 {
2715 strcpy(prj->code, "COP");
2716 prj->flag = copysign(WCS__COP, prj->flag);
2717 prj->phi0 = 0.0;
2718 prj->theta0 = prj->p[1];
2719
2720 if (prj->r0 == 0.0) prj->r0 = R2D;
2721
2722 prj->w[0] = astSind(prj->p[1]);
2723 if (prj->w[0] == 0.0) {
2724 return 1;
2725 }
2726
2727 prj->w[1] = 1.0/prj->w[0];
2728
2729 prj->w[3] = prj->r0*astCosd(prj->p[2]);
2730 if (prj->w[3] == 0.0) {
2731 return 1;
2732 }
2733
2734 prj->w[4] = 1.0/prj->w[3];
2735 prj->w[5] = 1.0/astTand(prj->p[1]);
2736
2737 prj->w[2] = prj->w[3]*prj->w[5];
2738
2739 prj->astPRJfwd = astCOPfwd;
2740 prj->astPRJrev = astCOPrev;
2741
2742 return 0;
2743 }
2744
2745 /*--------------------------------------------------------------------------*/
2746
astCOPfwd(phi,theta,prj,x,y)2747 int astCOPfwd(phi, theta, prj, x, y)
2748
2749 const double phi, theta;
2750 struct AstPrjPrm *prj;
2751 double *x, *y;
2752
2753 {
2754 double a, r, s, t;
2755
2756 if (abs(prj->flag) != WCS__COP) {
2757 if (astCOPset(prj)) return 1;
2758 }
2759
2760 t = theta - prj->p[1];
2761 s = astCosd(t);
2762 if (s == 0.0) {
2763 return 2;
2764 }
2765
2766 a = prj->w[0]*phi;
2767 r = prj->w[2] - prj->w[3]*astSind(t)/s;
2768
2769 *x = r*astSind(a);
2770 *y = prj->w[2] - r*astCosd(a);
2771
2772 if (prj->flag > 0 && r*prj->w[0] < 0.0) {
2773 return 2;
2774 }
2775
2776 return 0;
2777 }
2778
2779 /*--------------------------------------------------------------------------*/
2780
astCOPrev(x,y,prj,phi,theta)2781 int astCOPrev(x, y, prj, phi, theta)
2782
2783 const double x, y;
2784 struct AstPrjPrm *prj;
2785 double *phi, *theta;
2786
2787 {
2788 double a, dy, r;
2789
2790 if (abs(prj->flag) != WCS__COP) {
2791 if (astCOPset(prj)) return 1;
2792 }
2793
2794 dy = prj->w[2] - y;
2795 r = sqrt(x*x + dy*dy);
2796 if (prj->p[1] < 0.0) r = -r;
2797
2798 if (r == 0.0) {
2799 a = 0.0;
2800 } else {
2801 a = astATan2d(x/r, dy/r);
2802 }
2803
2804 *phi = a*prj->w[1];
2805 *theta = prj->p[1] + astATand(prj->w[5] - r*prj->w[4]);
2806
2807 return 0;
2808 }
2809
2810 /*============================================================================
2811 * COE: conic equal area projection.
2812 *
2813 * Given:
2814 * prj->p[1] sigma = (theta2+theta1)/2
2815 * prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the
2816 * latitudes of the standard parallels, in degrees.
2817 *
2818 * Given and/or returned:
2819 * prj->r0 r0; reset to 180/pi if 0.
2820 *
2821 * Returned:
2822 * prj->code "COE"
2823 * prj->flag COE
2824 * prj->phi0 0.0
2825 * prj->theta0 sigma
2826 * prj->w[0] C = (sin(theta1) + sin(theta2))/2
2827 * prj->w[1] 1/C
2828 * prj->w[2] Y0 = chi*sqrt(psi - 2C*astSind(sigma))
2829 * prj->w[3] chi = r0/C
2830 * prj->w[4] psi = 1 + sin(theta1)*sin(theta2)
2831 * prj->w[5] 2C
2832 * prj->w[6] (1 + sin(theta1)*sin(theta2))*(r0/C)**2
2833 * prj->w[7] C/(2*r0**2)
2834 * prj->w[8] chi*sqrt(psi + 2C)
2835 * prj->astPRJfwd Pointer to astCOEfwd().
2836 * prj->astPRJrev Pointer to astCOErev().
2837 *===========================================================================*/
2838
astCOEset(prj)2839 int astCOEset(prj)
2840
2841 struct AstPrjPrm *prj;
2842
2843 {
2844 double theta1, theta2;
2845
2846 strcpy(prj->code, "COE");
2847 prj->flag = WCS__COE;
2848 prj->phi0 = 0.0;
2849 prj->theta0 = prj->p[1];
2850
2851 if (prj->r0 == 0.0) prj->r0 = R2D;
2852
2853 theta1 = prj->p[1] - prj->p[2];
2854 theta2 = prj->p[1] + prj->p[2];
2855
2856 prj->w[0] = (astSind(theta1) + astSind(theta2))/2.0;
2857 if (prj->w[0] == 0.0) {
2858 return 1;
2859 }
2860
2861 prj->w[1] = 1.0/prj->w[0];
2862
2863 prj->w[3] = prj->r0/prj->w[0];
2864 prj->w[4] = 1.0 + astSind(theta1)*astSind(theta2);
2865 prj->w[5] = 2.0*prj->w[0];
2866 prj->w[6] = prj->w[3]*prj->w[3]*prj->w[4];
2867 prj->w[7] = 1.0/(2.0*prj->r0*prj->w[3]);
2868 prj->w[8] = prj->w[3]*sqrt(prj->w[4] + prj->w[5]);
2869
2870 prj->w[2] = prj->w[3]*sqrt(prj->w[4] - prj->w[5]*astSind(prj->p[1]));
2871
2872 prj->astPRJfwd = astCOEfwd;
2873 prj->astPRJrev = astCOErev;
2874
2875 return 0;
2876 }
2877
2878 /*--------------------------------------------------------------------------*/
2879
astCOEfwd(phi,theta,prj,x,y)2880 int astCOEfwd(phi, theta, prj, x, y)
2881
2882 const double phi, theta;
2883 struct AstPrjPrm *prj;
2884 double *x, *y;
2885
2886 {
2887 double a, r;
2888
2889 if (prj->flag != WCS__COE) {
2890 if (astCOEset(prj)) return 1;
2891 }
2892
2893 a = phi*prj->w[0];
2894 if (theta == -90.0) {
2895 r = prj->w[8];
2896 } else {
2897 r = prj->w[3]*sqrt(prj->w[4] - prj->w[5]*astSind(theta));
2898 }
2899
2900 *x = r*astSind(a);
2901 *y = prj->w[2] - r*astCosd(a);
2902
2903 return 0;
2904 }
2905
2906 /*--------------------------------------------------------------------------*/
2907
astCOErev(x,y,prj,phi,theta)2908 int astCOErev(x, y, prj, phi, theta)
2909
2910 const double x, y;
2911 struct AstPrjPrm *prj;
2912 double *phi, *theta;
2913
2914 {
2915 double a, dy, r, w;
2916 const double tol = 1.0e-12;
2917
2918 if (prj->flag != WCS__COE) {
2919 if (astCOEset(prj)) return 1;
2920 }
2921
2922 dy = prj->w[2] - y;
2923 r = sqrt(x*x + dy*dy);
2924 if (prj->p[1] < 0.0) r = -r;
2925
2926 if (r == 0.0) {
2927 a = 0.0;
2928 } else {
2929 a = astATan2d(x/r, dy/r);
2930 }
2931
2932 *phi = a*prj->w[1];
2933 if (fabs(r - prj->w[8]) < tol) {
2934 *theta = -90.0;
2935 } else {
2936 w = (prj->w[6] - r*r)*prj->w[7];
2937 if (fabs(w) > 1.0) {
2938 if (fabs(w-1.0) < tol) {
2939 *theta = 90.0;
2940 } else if (fabs(w+1.0) < tol) {
2941 *theta = -90.0;
2942 } else {
2943 return 2;
2944 }
2945 } else {
2946 *theta = astASind(w);
2947 }
2948 }
2949
2950 return 0;
2951 }
2952
2953 /*============================================================================
2954 * COD: conic equidistant projection.
2955 *
2956 * Given:
2957 * prj->p[1] sigma = (theta2+theta1)/2
2958 * prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the
2959 * latitudes of the standard parallels, in degrees.
2960 *
2961 * Given and/or returned:
2962 * prj->r0 r0; reset to 180/pi if 0.
2963 *
2964 * Returned:
2965 * prj->code "COD"
2966 * prj->flag COD
2967 * prj->phi0 0.0
2968 * prj->theta0 sigma
2969 * prj->w[0] C = r0*sin(sigma)*sin(delta)/delta
2970 * prj->w[1] 1/C
2971 * prj->w[2] Y0 = delta*cot(delta)*cot(sigma)
2972 * prj->w[3] Y0 + sigma
2973 * prj->astPRJfwd Pointer to astCODfwd().
2974 * prj->astPRJrev Pointer to astCODrev().
2975 *===========================================================================*/
2976
astCODset(prj)2977 int astCODset(prj)
2978
2979 struct AstPrjPrm *prj;
2980
2981 {
2982 strcpy(prj->code, "COD");
2983 prj->flag = WCS__COD;
2984 prj->phi0 = 0.0;
2985 prj->theta0 = prj->p[1];
2986
2987 if (prj->r0 == 0.0) prj->r0 = R2D;
2988
2989 if (prj->p[2] == 0.0) {
2990 prj->w[0] = prj->r0*astSind(prj->p[1])*D2R;
2991 } else {
2992 prj->w[0] = prj->r0*astSind(prj->p[1])*astSind(prj->p[2])/prj->p[2];
2993 }
2994
2995 if (prj->w[0] == 0.0) {
2996 return 1;
2997 }
2998
2999 prj->w[1] = 1.0/prj->w[0];
3000 prj->w[2] = prj->r0*astCosd(prj->p[2])*astCosd(prj->p[1])/prj->w[0];
3001 prj->w[3] = prj->w[2] + prj->p[1];
3002
3003 prj->astPRJfwd = astCODfwd;
3004 prj->astPRJrev = astCODrev;
3005
3006 return 0;
3007 }
3008
3009 /*--------------------------------------------------------------------------*/
3010
astCODfwd(phi,theta,prj,x,y)3011 int astCODfwd(phi, theta, prj, x, y)
3012
3013 const double phi, theta;
3014 struct AstPrjPrm *prj;
3015 double *x, *y;
3016
3017 {
3018 double a, r;
3019
3020 if (prj->flag != WCS__COD) {
3021 if (astCODset(prj)) return 1;
3022 }
3023
3024 a = prj->w[0]*phi;
3025 r = prj->w[3] - theta;
3026
3027 *x = r*astSind(a);
3028 *y = prj->w[2] - r*astCosd(a);
3029
3030 return 0;
3031 }
3032
3033 /*--------------------------------------------------------------------------*/
3034
astCODrev(x,y,prj,phi,theta)3035 int astCODrev(x, y, prj, phi, theta)
3036
3037 const double x, y;
3038 struct AstPrjPrm *prj;
3039 double *phi, *theta;
3040
3041 {
3042 double a, dy, r;
3043
3044 if (prj->flag != WCS__COD) {
3045 if (astCODset(prj)) return 1;
3046 }
3047
3048 dy = prj->w[2] - y;
3049 r = sqrt(x*x + dy*dy);
3050 if (prj->p[1] < 0.0) r = -r;
3051
3052 if (r == 0.0) {
3053 a = 0.0;
3054 } else {
3055 a = astATan2d(x/r, dy/r);
3056 }
3057
3058 *phi = a*prj->w[1];
3059 *theta = prj->w[3] - r;
3060
3061 return 0;
3062 }
3063
3064 /*============================================================================
3065 * COO: conic orthomorphic projection.
3066 *
3067 * Given:
3068 * prj->p[1] sigma = (theta2+theta1)/2
3069 * prj->p[2] delta = (theta2-theta1)/2, where theta1 and theta2 are the
3070 * latitudes of the standard parallels, in degrees.
3071 *
3072 * Given and/or returned:
3073 * prj->r0 r0; reset to 180/pi if 0.
3074 *
3075 * Returned:
3076 * prj->code "COO"
3077 * prj->flag COO
3078 * prj->phi0 0.0
3079 * prj->theta0 sigma
3080 * prj->w[0] C = ln(cos(theta2)/cos(theta1))/ln(tan(tau2)/tan(tau1))
3081 * where tau1 = (90 - theta1)/2
3082 * tau2 = (90 - theta2)/2
3083 * prj->w[1] 1/C
3084 * prj->w[2] Y0 = psi*tan((90-sigma)/2)**C
3085 * prj->w[3] psi = (r0*cos(theta1)/C)/tan(tau1)**C
3086 * prj->w[4] 1/psi
3087 * prj->astPRJfwd Pointer to astCOOfwd().
3088 * prj->astPRJrev Pointer to astCOOrev().
3089 *===========================================================================*/
3090
astCOOset(prj)3091 int astCOOset(prj)
3092
3093 struct AstPrjPrm *prj;
3094
3095 {
3096 double cos1, cos2, tan1, tan2, theta1, theta2;
3097
3098 strcpy(prj->code, "COO");
3099 prj->flag = WCS__COO;
3100 prj->phi0 = 0.0;
3101 prj->theta0 = prj->p[1];
3102
3103 if (prj->r0 == 0.0) prj->r0 = R2D;
3104
3105 theta1 = prj->p[1] - prj->p[2];
3106 theta2 = prj->p[1] + prj->p[2];
3107
3108 tan1 = astTand((90.0 - theta1)/2.0);
3109 cos1 = astCosd(theta1);
3110
3111 if (theta1 == theta2) {
3112 prj->w[0] = astSind(theta1);
3113 } else {
3114 tan2 = astTand((90.0 - theta2)/2.0);
3115 cos2 = astCosd(theta2);
3116 prj->w[0] = log(cos2/cos1)/log(tan2/tan1);
3117 }
3118 if (prj->w[0] == 0.0) {
3119 return 1;
3120 }
3121
3122 prj->w[1] = 1.0/prj->w[0];
3123
3124 prj->w[3] = prj->r0*(cos1/prj->w[0])/pow(tan1,prj->w[0]);
3125 if (prj->w[3] == 0.0) {
3126 return 1;
3127 }
3128 prj->w[2] = prj->w[3]*pow(astTand((90.0 - prj->p[1])/2.0),prj->w[0]);
3129 prj->w[4] = 1.0/prj->w[3];
3130
3131 prj->astPRJfwd = astCOOfwd;
3132 prj->astPRJrev = astCOOrev;
3133
3134 return 0;
3135 }
3136
3137 /*--------------------------------------------------------------------------*/
3138
astCOOfwd(phi,theta,prj,x,y)3139 int astCOOfwd(phi, theta, prj, x, y)
3140
3141 const double phi, theta;
3142 struct AstPrjPrm *prj;
3143 double *x, *y;
3144
3145 {
3146 double a, r;
3147
3148 if (prj->flag != WCS__COO) {
3149 if (astCOOset(prj)) return 1;
3150 }
3151
3152 a = prj->w[0]*phi;
3153 if (theta == -90.0) {
3154 if (prj->w[0] < 0.0) {
3155 r = 0.0;
3156 } else {
3157 return 2;
3158 }
3159 } else {
3160 r = prj->w[3]*pow(astTand((90.0 - theta)/2.0),prj->w[0]);
3161 }
3162
3163 *x = r*astSind(a);
3164 *y = prj->w[2] - r*astCosd(a);
3165
3166 return 0;
3167 }
3168
3169 /*--------------------------------------------------------------------------*/
3170
astCOOrev(x,y,prj,phi,theta)3171 int astCOOrev(x, y, prj, phi, theta)
3172
3173 const double x, y;
3174 struct AstPrjPrm *prj;
3175 double *phi, *theta;
3176
3177 {
3178 double a, dy, r;
3179
3180 if (prj->flag != WCS__COO) {
3181 if (astCOOset(prj)) return 1;
3182 }
3183
3184 dy = prj->w[2] - y;
3185 r = sqrt(x*x + dy*dy);
3186 if (prj->p[1] < 0.0) r = -r;
3187
3188 if (r == 0.0) {
3189 a = 0.0;
3190 } else {
3191 a = astATan2d(x/r, dy/r);
3192 }
3193
3194 *phi = a*prj->w[1];
3195 if (r == 0.0) {
3196 if (prj->w[0] < 0.0) {
3197 *theta = -90.0;
3198 } else {
3199 return 2;
3200 }
3201 } else {
3202 *theta = 90.0 - 2.0*astATand(pow(r*prj->w[4],prj->w[1]));
3203 }
3204
3205 return 0;
3206 }
3207
3208 /*============================================================================
3209 * BON: Bonne's projection.
3210 *
3211 * Given:
3212 * prj->p[1] Bonne conformal latitude, theta1, in degrees.
3213 *
3214 * Given and/or returned:
3215 * prj->r0 r0; reset to 180/pi if 0.
3216 *
3217 * Returned:
3218 * prj->code "BON"
3219 * prj->flag BON
3220 * prj->phi0 0.0
3221 * prj->theta0 0.0
3222 * prj->w[1] r0*pi/180
3223 * prj->w[2] Y0 = r0*(cot(theta1) + theta1*pi/180)
3224 * prj->astPRJfwd Pointer to astBONfwd().
3225 * prj->astPRJrev Pointer to astBONrev().
3226 *===========================================================================*/
3227
astBONset(prj)3228 int astBONset(prj)
3229
3230 struct AstPrjPrm *prj;
3231
3232 {
3233 strcpy(prj->code, "BON");
3234 prj->flag = WCS__BON;
3235 prj->phi0 = 0.0;
3236 prj->theta0 = 0.0;
3237
3238 if (prj->r0 == 0.0) {
3239 prj->r0 = R2D;
3240 prj->w[1] = 1.0;
3241 prj->w[2] = prj->r0*astCosd(prj->p[1])/astSind(prj->p[1]) + prj->p[1];
3242 } else {
3243 prj->w[1] = prj->r0*D2R;
3244 prj->w[2] = prj->r0*(astCosd(prj->p[1])/astSind(prj->p[1]) + prj->p[1]*D2R);
3245 }
3246
3247 prj->astPRJfwd = astBONfwd;
3248 prj->astPRJrev = astBONrev;
3249
3250 return 0;
3251 }
3252
3253 /*--------------------------------------------------------------------------*/
3254
astBONfwd(phi,theta,prj,x,y)3255 int astBONfwd(phi, theta, prj, x, y)
3256
3257 const double phi, theta;
3258 struct AstPrjPrm *prj;
3259 double *x, *y;
3260
3261 {
3262 double a, r;
3263
3264 if (prj->p[1] == 0.0) {
3265 /* Sanson-Flamsteed. */
3266 return astSFLfwd(phi, theta, prj, x, y);
3267 }
3268
3269 if (prj->flag != WCS__BON) {
3270 if (astBONset(prj)) return 1;
3271 }
3272
3273 r = prj->w[2] - theta*prj->w[1];
3274 a = prj->r0*phi*astCosd(theta)/r;
3275
3276 *x = r*astSind(a);
3277 *y = prj->w[2] - r*astCosd(a);
3278
3279 return 0;
3280 }
3281
3282 /*--------------------------------------------------------------------------*/
3283
astBONrev(x,y,prj,phi,theta)3284 int astBONrev(x, y, prj, phi, theta)
3285
3286 const double x, y;
3287 struct AstPrjPrm *prj;
3288 double *phi, *theta;
3289
3290 {
3291 double a, cthe, dy, r;
3292
3293 if (prj->p[1] == 0.0) {
3294 /* Sanson-Flamsteed. */
3295 return astSFLrev(x, y, prj, phi, theta);
3296 }
3297
3298 if (prj->flag != WCS__BON) {
3299 if (astBONset(prj)) return 1;
3300 }
3301
3302 dy = prj->w[2] - y;
3303 r = sqrt(x*x + dy*dy);
3304 if (prj->p[1] < 0.0) r = -r;
3305
3306 if (r == 0.0) {
3307 a = 0.0;
3308 } else {
3309 a = astATan2d(x/r, dy/r);
3310 }
3311
3312 *theta = (prj->w[2] - r)/prj->w[1];
3313 cthe = astCosd(*theta);
3314 if (cthe == 0.0) {
3315 *phi = 0.0;
3316 } else {
3317 *phi = a*(r/prj->r0)/cthe;
3318 }
3319
3320 return 0;
3321 }
3322
3323 /*============================================================================
3324 * PCO: polyconic projection.
3325 *
3326 * Given and/or returned:
3327 * prj->r0 r0; reset to 180/pi if 0.
3328 *
3329 * Returned:
3330 * prj->code "PCO"
3331 * prj->flag PCO
3332 * prj->phi0 0.0
3333 * prj->theta0 0.0
3334 * prj->w[0] r0*(pi/180)
3335 * prj->w[1] 1/r0
3336 * prj->w[2] 2*r0
3337 * prj->astPRJfwd Pointer to astPCOfwd().
3338 * prj->astPRJrev Pointer to astPCOrev().
3339 *===========================================================================*/
3340
astPCOset(prj)3341 int astPCOset(prj)
3342
3343 struct AstPrjPrm *prj;
3344
3345 {
3346 strcpy(prj->code, "PCO");
3347 prj->flag = WCS__PCO;
3348 prj->phi0 = 0.0;
3349 prj->theta0 = 0.0;
3350
3351 if (prj->r0 == 0.0) {
3352 prj->r0 = R2D;
3353 prj->w[0] = 1.0;
3354 prj->w[1] = 1.0;
3355 prj->w[2] = 360.0/PI;
3356 } else {
3357 prj->w[0] = prj->r0*D2R;
3358 prj->w[1] = 1.0/prj->w[0];
3359 prj->w[2] = 2.0*prj->r0;
3360 }
3361
3362 prj->astPRJfwd = astPCOfwd;
3363 prj->astPRJrev = astPCOrev;
3364
3365 return 0;
3366 }
3367
3368 /*--------------------------------------------------------------------------*/
3369
astPCOfwd(phi,theta,prj,x,y)3370 int astPCOfwd(phi, theta, prj, x, y)
3371
3372 const double phi, theta;
3373 struct AstPrjPrm *prj;
3374 double *x, *y;
3375
3376 {
3377 double a, cthe, cotthe, sthe;
3378
3379 if (prj->flag != WCS__PCO) {
3380 if (astPCOset(prj)) return 1;
3381 }
3382
3383 cthe = astCosd(theta);
3384 sthe = astSind(theta);
3385 a = phi*sthe;
3386
3387 if (sthe == 0.0) {
3388 *x = prj->w[0]*phi;
3389 *y = 0.0;
3390 } else {
3391 cotthe = cthe/sthe;
3392 *x = prj->r0*cotthe*astSind(a);
3393 *y = prj->r0*(cotthe*(1.0 - astCosd(a)) + theta*D2R);
3394 }
3395
3396 return 0;
3397 }
3398
3399 /*--------------------------------------------------------------------------*/
3400
astPCOrev(x,y,prj,phi,theta)3401 int astPCOrev(x, y, prj, phi, theta)
3402
3403 const double x, y;
3404 struct AstPrjPrm *prj;
3405 double *phi, *theta;
3406
3407 {
3408 int j;
3409 double f, fneg, fpos, lambda, tanthe, theneg, thepos, w, xp, xx, ymthe, yp;
3410 const double tol = 1.0e-12;
3411
3412 if (prj->flag != WCS__PCO) {
3413 if (astPCOset(prj)) return 1;
3414 }
3415
3416 w = fabs(y*prj->w[1]);
3417 if (w < tol) {
3418 *phi = x*prj->w[1];
3419 *theta = 0.0;
3420 } else if (fabs(w-90.0) < tol) {
3421 *phi = 0.0;
3422 *theta = copysign(90.0,y);
3423 } else {
3424 /* Iterative solution using weighted division of the interval. */
3425 if (y > 0.0) {
3426 thepos = 90.0;
3427 } else {
3428 thepos = -90.0;
3429 }
3430 theneg = 0.0;
3431
3432 xx = x*x;
3433 ymthe = y - prj->w[0]*thepos;
3434 fpos = xx + ymthe*ymthe;
3435 fneg = -999.0;
3436
3437 for (j = 0; j < 64; j++) {
3438 if (fneg < -100.0) {
3439 /* Equal division of the interval. */
3440 *theta = (thepos+theneg)/2.0;
3441 } else {
3442 /* Weighted division of the interval. */
3443 lambda = fpos/(fpos-fneg);
3444 if (lambda < 0.1) {
3445 lambda = 0.1;
3446 } else if (lambda > 0.9) {
3447 lambda = 0.9;
3448 }
3449 *theta = thepos - lambda*(thepos-theneg);
3450 }
3451
3452 /* Compute the residue. */
3453 ymthe = y - prj->w[0]*(*theta);
3454 tanthe = astTand(*theta);
3455 f = xx + ymthe*(ymthe - prj->w[2]/tanthe);
3456
3457 /* Check for convergence. */
3458 if (fabs(f) < tol) break;
3459 if (fabs(thepos-theneg) < tol) break;
3460
3461 /* Redefine the interval. */
3462 if (f > 0.0) {
3463 thepos = *theta;
3464 fpos = f;
3465 } else {
3466 theneg = *theta;
3467 fneg = f;
3468 }
3469 }
3470
3471 xp = prj->r0 - ymthe*tanthe;
3472 yp = x*tanthe;
3473 if (xp == 0.0 && yp == 0.0) {
3474 *phi = 0.0;
3475 } else {
3476 *phi = astATan2d(yp, xp)/astSind(*theta);
3477 }
3478 }
3479
3480 return 0;
3481 }
3482
3483 /*============================================================================
3484 * TSC: tangential spherical cube projection.
3485 *
3486 * Given and/or returned:
3487 * prj->r0 r0; reset to 180/pi if 0.
3488 *
3489 * Returned:
3490 * prj->code "TSC"
3491 * prj->flag TSC
3492 * prj->phi0 0.0
3493 * prj->theta0 0.0
3494 * prj->w[0] r0*(pi/4)
3495 * prj->w[1] (4/pi)/r0
3496 * prj->astPRJfwd Pointer to astTSCfwd().
3497 * prj->astPRJrev Pointer to astTSCrev().
3498 *===========================================================================*/
3499
astTSCset(prj)3500 int astTSCset(prj)
3501
3502 struct AstPrjPrm *prj;
3503
3504 {
3505 strcpy(prj->code, "TSC");
3506 prj->flag = WCS__TSC;
3507 prj->phi0 = 0.0;
3508 prj->theta0 = 0.0;
3509
3510 if (prj->r0 == 0.0) {
3511 prj->r0 = R2D;
3512 prj->w[0] = 45.0;
3513 prj->w[1] = 1.0/45.0;
3514 } else {
3515 prj->w[0] = prj->r0*PI/4.0;
3516 prj->w[1] = 1.0/prj->w[0];
3517 }
3518
3519 prj->astPRJfwd = astTSCfwd;
3520 prj->astPRJrev = astTSCrev;
3521
3522 return 0;
3523 }
3524
3525 /*--------------------------------------------------------------------------*/
3526
astTSCfwd(phi,theta,prj,x,y)3527 int astTSCfwd(phi, theta, prj, x, y)
3528
3529 const double phi, theta;
3530 struct AstPrjPrm *prj;
3531 double *x, *y;
3532
3533 {
3534 int face;
3535 double cthe, l, m, n, rho, x0, xf, y0, yf;
3536 const double tol = 1.0e-12;
3537
3538 x0 = 0.0;
3539 xf = 0.0;
3540 y0 = 0.0;
3541 yf = 0.0;
3542
3543 if (prj->flag != WCS__TSC) {
3544 if (astTSCset(prj)) return 1;
3545 }
3546
3547 cthe = astCosd(theta);
3548 l = cthe*astCosd(phi);
3549 m = cthe*astSind(phi);
3550 n = astSind(theta);
3551
3552 face = 0;
3553 rho = n;
3554 if (l > rho) {
3555 face = 1;
3556 rho = l;
3557 }
3558 if (m > rho) {
3559 face = 2;
3560 rho = m;
3561 }
3562 if (-l > rho) {
3563 face = 3;
3564 rho = -l;
3565 }
3566 if (-m > rho) {
3567 face = 4;
3568 rho = -m;
3569 }
3570 if (-n > rho) {
3571 face = 5;
3572 rho = -n;
3573 }
3574
3575 if (face == 0) {
3576 xf = m/rho;
3577 yf = -l/rho;
3578 x0 = 0.0;
3579 y0 = 2.0;
3580 } else if (face == 1) {
3581 xf = m/rho;
3582 yf = n/rho;
3583 x0 = 0.0;
3584 y0 = 0.0;
3585 } else if (face == 2) {
3586 xf = -l/rho;
3587 yf = n/rho;
3588 x0 = 2.0;
3589 y0 = 0.0;
3590 } else if (face == 3) {
3591 xf = -m/rho;
3592 yf = n/rho;
3593 x0 = 4.0;
3594 y0 = 0.0;
3595 } else if (face == 4) {
3596 xf = l/rho;
3597 yf = n/rho;
3598 x0 = 6.0;
3599 y0 = 0.0;
3600 } else if (face == 5) {
3601 xf = m/rho;
3602 yf = l/rho;
3603 x0 = 0.0;
3604 y0 = -2.0;
3605 }
3606
3607 if (fabs(xf) > 1.0) {
3608 if (fabs(xf) > 1.0+tol) {
3609 return 2;
3610 }
3611 xf = copysign(1.0,xf);
3612 }
3613 if (fabs(yf) > 1.0) {
3614 if (fabs(yf) > 1.0+tol) {
3615 return 2;
3616 }
3617 yf = copysign(1.0,yf);
3618 }
3619
3620 *x = prj->w[0]*(xf + x0);
3621 *y = prj->w[0]*(yf + y0);
3622
3623 return 0;
3624 }
3625
3626 /*--------------------------------------------------------------------------*/
3627
astTSCrev(x,y,prj,phi,theta)3628 int astTSCrev(x, y, prj, phi, theta)
3629
3630 const double x, y;
3631 struct AstPrjPrm *prj;
3632 double *phi, *theta;
3633
3634 {
3635 double l, m, n, xf, yf;
3636
3637 if (prj->flag != WCS__TSC) {
3638 if (astTSCset(prj)) return 1;
3639 }
3640
3641 xf = x*prj->w[1];
3642 yf = y*prj->w[1];
3643
3644 /* Check bounds. */
3645 if (fabs(xf) <= 1.0) {
3646 if (fabs(yf) > 3.0) return 2;
3647 } else {
3648 if (fabs(xf) > 7.0) return 2;
3649 if (fabs(yf) > 1.0) return 2;
3650 }
3651
3652 /* Map negative faces to the other side. */
3653 if (xf < -1.0) xf += 8.0;
3654
3655 /* Determine the face. */
3656 if (xf > 5.0) {
3657 /* face = 4 */
3658 xf = xf - 6.0;
3659 m = -1.0/sqrt(1.0 + xf*xf + yf*yf);
3660 l = -m*xf;
3661 n = -m*yf;
3662 } else if (xf > 3.0) {
3663 /* face = 3 */
3664 xf = xf - 4.0;
3665 l = -1.0/sqrt(1.0 + xf*xf + yf*yf);
3666 m = l*xf;
3667 n = -l*yf;
3668 } else if (xf > 1.0) {
3669 /* face = 2 */
3670 xf = xf - 2.0;
3671 m = 1.0/sqrt(1.0 + xf*xf + yf*yf);
3672 l = -m*xf;
3673 n = m*yf;
3674 } else if (yf > 1.0) {
3675 /* face = 0 */
3676 yf = yf - 2.0;
3677 n = 1.0/sqrt(1.0 + xf*xf + yf*yf);
3678 l = -n*yf;
3679 m = n*xf;
3680 } else if (yf < -1.0) {
3681 /* face = 5 */
3682 yf = yf + 2.0;
3683 n = -1.0/sqrt(1.0 + xf*xf + yf*yf);
3684 l = -n*yf;
3685 m = -n*xf;
3686 } else {
3687 /* face = 1 */
3688 l = 1.0/sqrt(1.0 + xf*xf + yf*yf);
3689 m = l*xf;
3690 n = l*yf;
3691 }
3692
3693 if (l == 0.0 && m == 0.0) {
3694 *phi = 0.0;
3695 } else {
3696 *phi = astATan2d(m, l);
3697 }
3698 *theta = astASind(n);
3699
3700 return 0;
3701 }
3702
3703 /*============================================================================
3704 * CSC: COBE quadrilateralized spherical cube projection.
3705 *
3706 * Given and/or returned:
3707 * prj->r0 r0; reset to 180/pi if 0.
3708 *
3709 * Returned:
3710 * prj->code "CSC"
3711 * prj->flag CSC
3712 * prj->phi0 0.0
3713 * prj->theta0 0.0
3714 * prj->w[0] r0*(pi/4)
3715 * prj->w[1] (4/pi)/r0
3716 * prj->astPRJfwd Pointer to astCSCfwd().
3717 * prj->astPRJrev Pointer to astCSCrev().
3718 *===========================================================================*/
3719
astCSCset(prj)3720 int astCSCset(prj)
3721
3722 struct AstPrjPrm *prj;
3723
3724 {
3725 strcpy(prj->code, "CSC");
3726 prj->flag = WCS__CSC;
3727 prj->phi0 = 0.0;
3728 prj->theta0 = 0.0;
3729
3730 if (prj->r0 == 0.0) {
3731 prj->r0 = R2D;
3732 prj->w[0] = 45.0;
3733 prj->w[1] = 1.0/45.0;
3734 } else {
3735 prj->w[0] = prj->r0*PI/4.0;
3736 prj->w[1] = 1.0/prj->w[0];
3737 }
3738
3739 prj->astPRJfwd = astCSCfwd;
3740 prj->astPRJrev = astCSCrev;
3741
3742 return 0;
3743 }
3744
3745 /*--------------------------------------------------------------------------*/
3746
astCSCfwd(phi,theta,prj,x,y)3747 int astCSCfwd(phi, theta, prj, x, y)
3748
3749 const double phi, theta;
3750 struct AstPrjPrm *prj;
3751 double *x, *y;
3752
3753 {
3754 int face;
3755 float cthe, eta, l, m, n, rho, xi;
3756 const float tol = 1.0e-7;
3757
3758 float a, a2, a2b2, a4, ab, b, b2, b4, ca2, cb2, x0, xf, y0, yf;
3759 const float gstar = 1.37484847732;
3760 const float mm = 0.004869491981;
3761 const float gamma = -0.13161671474;
3762 const float omega1 = -0.159596235474;
3763 const float d0 = 0.0759196200467;
3764 const float d1 = -0.0217762490699;
3765 const float c00 = 0.141189631152;
3766 const float c10 = 0.0809701286525;
3767 const float c01 = -0.281528535557;
3768 const float c11 = 0.15384112876;
3769 const float c20 = -0.178251207466;
3770 const float c02 = 0.106959469314;
3771
3772 eta = 0.0;
3773 xi = 0.0;
3774 x0 = 0.0;
3775 y0 = 0.0;
3776
3777 if (prj->flag != WCS__CSC) {
3778 if (astCSCset(prj)) return 1;
3779 }
3780
3781 cthe = astCosd(theta);
3782 l = cthe*astCosd(phi);
3783 m = cthe*astSind(phi);
3784 n = astSind(theta);
3785
3786 face = 0;
3787 rho = n;
3788 if (l > rho) {
3789 face = 1;
3790 rho = l;
3791 }
3792 if (m > rho) {
3793 face = 2;
3794 rho = m;
3795 }
3796 if (-l > rho) {
3797 face = 3;
3798 rho = -l;
3799 }
3800 if (-m > rho) {
3801 face = 4;
3802 rho = -m;
3803 }
3804 if (-n > rho) {
3805 face = 5;
3806 rho = -n;
3807 }
3808
3809 if (face == 0) {
3810 xi = m;
3811 eta = -l;
3812 x0 = 0.0;
3813 y0 = 2.0;
3814 } else if (face == 1) {
3815 xi = m;
3816 eta = n;
3817 x0 = 0.0;
3818 y0 = 0.0;
3819 } else if (face == 2) {
3820 xi = -l;
3821 eta = n;
3822 x0 = 2.0;
3823 y0 = 0.0;
3824 } else if (face == 3) {
3825 xi = -m;
3826 eta = n;
3827 x0 = 4.0;
3828 y0 = 0.0;
3829 } else if (face == 4) {
3830 xi = l;
3831 eta = n;
3832 x0 = 6.0;
3833 y0 = 0.0;
3834 } else if (face == 5) {
3835 xi = m;
3836 eta = l;
3837 x0 = 0.0;
3838 y0 = -2.0;
3839 }
3840
3841 a = xi/rho;
3842 b = eta/rho;
3843
3844 a2 = a*a;
3845 b2 = b*b;
3846 ca2 = 1.0 - a2;
3847 cb2 = 1.0 - b2;
3848
3849 /* Avoid floating underflows. */
3850 ab = fabs(a*b);
3851 a4 = (a2 > 1.0e-16) ? a2*a2 : 0.0;
3852 b4 = (b2 > 1.0e-16) ? b2*b2 : 0.0;
3853 a2b2 = (ab > 1.0e-16) ? a2*b2 : 0.0;
3854
3855 xf = a*(a2 + ca2*(gstar + b2*(gamma*ca2 + mm*a2 +
3856 cb2*(c00 + c10*a2 + c01*b2 + c11*a2b2 + c20*a4 + c02*b4)) +
3857 a2*(omega1 - ca2*(d0 + d1*a2))));
3858 yf = b*(b2 + cb2*(gstar + a2*(gamma*cb2 + mm*b2 +
3859 ca2*(c00 + c10*b2 + c01*a2 + c11*a2b2 + c20*b4 + c02*a4)) +
3860 b2*(omega1 - cb2*(d0 + d1*b2))));
3861
3862 if (fabs(xf) > 1.0) {
3863 if (fabs(xf) > 1.0+tol) {
3864 return 2;
3865 }
3866 xf = copysign(1.0,xf);
3867 }
3868 if (fabs(yf) > 1.0) {
3869 if (fabs(yf) > 1.0+tol) {
3870 return 2;
3871 }
3872 yf = copysign(1.0,yf);
3873 }
3874
3875 *x = prj->w[0]*(x0 + xf);
3876 *y = prj->w[0]*(y0 + yf);
3877
3878 return 0;
3879 }
3880
3881 /*--------------------------------------------------------------------------*/
3882
astCSCrev(x,y,prj,phi,theta)3883 int astCSCrev(x, y, prj, phi, theta)
3884
3885 const double x, y;
3886 struct AstPrjPrm *prj;
3887 double *phi, *theta;
3888
3889 {
3890 int face;
3891 float l, m, n;
3892
3893 float a, b, xf, xx, yf, yy, z0, z1, z2, z3, z4, z5, z6;
3894 const float p00 = -0.27292696;
3895 const float p10 = -0.07629969;
3896 const float p20 = -0.22797056;
3897 const float p30 = 0.54852384;
3898 const float p40 = -0.62930065;
3899 const float p50 = 0.25795794;
3900 const float p60 = 0.02584375;
3901 const float p01 = -0.02819452;
3902 const float p11 = -0.01471565;
3903 const float p21 = 0.48051509;
3904 const float p31 = -1.74114454;
3905 const float p41 = 1.71547508;
3906 const float p51 = -0.53022337;
3907 const float p02 = 0.27058160;
3908 const float p12 = -0.56800938;
3909 const float p22 = 0.30803317;
3910 const float p32 = 0.98938102;
3911 const float p42 = -0.83180469;
3912 const float p03 = -0.60441560;
3913 const float p13 = 1.50880086;
3914 const float p23 = -0.93678576;
3915 const float p33 = 0.08693841;
3916 const float p04 = 0.93412077;
3917 const float p14 = -1.41601920;
3918 const float p24 = 0.33887446;
3919 const float p05 = -0.63915306;
3920 const float p15 = 0.52032238;
3921 const float p06 = 0.14381585;
3922
3923 l = 0.0;
3924 m = 0.0;
3925 n = 0.0;
3926
3927 if (prj->flag != WCS__CSC) {
3928 if (astCSCset(prj)) return 1;
3929 }
3930
3931 xf = x*prj->w[1];
3932 yf = y*prj->w[1];
3933
3934 /* Check bounds. */
3935 if (fabs(xf) <= 1.0) {
3936 if (fabs(yf) > 3.0) return 2;
3937 } else {
3938 if (fabs(xf) > 7.0) return 2;
3939 if (fabs(yf) > 1.0) return 2;
3940 }
3941
3942 /* Map negative faces to the other side. */
3943 if (xf < -1.0) xf += 8.0;
3944
3945 /* Determine the face. */
3946 if (xf > 5.0) {
3947 face = 4;
3948 xf = xf - 6.0;
3949 } else if (xf > 3.0) {
3950 face = 3;
3951 xf = xf - 4.0;
3952 } else if (xf > 1.0) {
3953 face = 2;
3954 xf = xf - 2.0;
3955 } else if (yf > 1.0) {
3956 face = 0;
3957 yf = yf - 2.0;
3958 } else if (yf < -1.0) {
3959 face = 5;
3960 yf = yf + 2.0;
3961 } else {
3962 face = 1;
3963 }
3964
3965 xx = xf*xf;
3966 yy = yf*yf;
3967
3968 z0 = p00 + xx*(p10 + xx*(p20 + xx*(p30 + xx*(p40 + xx*(p50 + xx*(p60))))));
3969 z1 = p01 + xx*(p11 + xx*(p21 + xx*(p31 + xx*(p41 + xx*(p51)))));
3970 z2 = p02 + xx*(p12 + xx*(p22 + xx*(p32 + xx*(p42))));
3971 z3 = p03 + xx*(p13 + xx*(p23 + xx*(p33)));
3972 z4 = p04 + xx*(p14 + xx*(p24));
3973 z5 = p05 + xx*(p15);
3974 z6 = p06;
3975
3976 a = z0 + yy*(z1 + yy*(z2 + yy*(z3 + yy*(z4 + yy*(z5 + yy*z6)))));
3977 a = xf + xf*(1.0 - xx)*a;
3978
3979 z0 = p00 + yy*(p10 + yy*(p20 + yy*(p30 + yy*(p40 + yy*(p50 + yy*(p60))))));
3980 z1 = p01 + yy*(p11 + yy*(p21 + yy*(p31 + yy*(p41 + yy*(p51)))));
3981 z2 = p02 + yy*(p12 + yy*(p22 + yy*(p32 + yy*(p42))));
3982 z3 = p03 + yy*(p13 + yy*(p23 + yy*(p33)));
3983 z4 = p04 + yy*(p14 + yy*(p24));
3984 z5 = p05 + yy*(p15);
3985 z6 = p06;
3986
3987 b = z0 + xx*(z1 + xx*(z2 + xx*(z3 + xx*(z4 + xx*(z5 + xx*z6)))));
3988 b = yf + yf*(1.0 - yy)*b;
3989
3990 if (face == 0) {
3991 n = 1.0/sqrt(a*a + b*b + 1.0);
3992 l = -b*n;
3993 m = a*n;
3994 } else if (face == 1) {
3995 l = 1.0/sqrt(a*a + b*b + 1.0);
3996 m = a*l;
3997 n = b*l;
3998 } else if (face == 2) {
3999 m = 1.0/sqrt(a*a + b*b + 1.0);
4000 l = -a*m;
4001 n = b*m;
4002 } else if (face == 3) {
4003 l = -1.0/sqrt(a*a + b*b + 1.0);
4004 m = a*l;
4005 n = -b*l;
4006 } else if (face == 4) {
4007 m = -1.0/sqrt(a*a + b*b + 1.0);
4008 l = -a*m;
4009 n = -b*m;
4010 } else if (face == 5) {
4011 n = -1.0/sqrt(a*a + b*b + 1.0);
4012 l = -b*n;
4013 m = -a*n;
4014 }
4015
4016 if (l == 0.0 && m == 0.0) {
4017 *phi = 0.0;
4018 } else {
4019 *phi = astATan2d(m, l);
4020 }
4021 *theta = astASind(n);
4022
4023 return 0;
4024 }
4025
4026 /*============================================================================
4027 * QSC: quadrilaterilized spherical cube projection.
4028 *
4029 * Given and/or returned:
4030 * prj->r0 r0; reset to 180/pi if 0.
4031 *
4032 * Returned:
4033 * prj->code "QSC"
4034 * prj->flag QSC
4035 * prj->phi0 0.0
4036 * prj->theta0 0.0
4037 * prj->w[0] r0*(pi/4)
4038 * prj->w[1] (4/pi)/r0
4039 * prj->astPRJfwd Pointer to astQSCfwd().
4040 * prj->astPRJrev Pointer to astQSCrev().
4041 *===========================================================================*/
4042
astQSCset(prj)4043 int astQSCset(prj)
4044
4045 struct AstPrjPrm *prj;
4046
4047 {
4048 strcpy(prj->code, "QSC");
4049 prj->flag = WCS__QSC;
4050 prj->phi0 = 0.0;
4051 prj->theta0 = 0.0;
4052
4053 if (prj->r0 == 0.0) {
4054 prj->r0 = R2D;
4055 prj->w[0] = 45.0;
4056 prj->w[1] = 1.0/45.0;
4057 } else {
4058 prj->w[0] = prj->r0*PI/4.0;
4059 prj->w[1] = 1.0/prj->w[0];
4060 }
4061
4062 prj->astPRJfwd = astQSCfwd;
4063 prj->astPRJrev = astQSCrev;
4064
4065 return 0;
4066 }
4067
4068 /*--------------------------------------------------------------------------*/
4069
astQSCfwd(phi,theta,prj,x,y)4070 int astQSCfwd(phi, theta, prj, x, y)
4071
4072 const double phi, theta;
4073 struct AstPrjPrm *prj;
4074 double *x, *y;
4075
4076 {
4077 int face;
4078 double cthe, eta, l, m, n, omega, p, rho, rhu, t, tau, x0, xf, xi, y0, yf;
4079 const double tol = 1.0e-12;
4080
4081 eta = 0.0;
4082 x0 = 0.0;
4083 xf = 0.0;
4084 xi = 0.0;
4085 y0 = 0.0;
4086 yf = 0.0;
4087
4088 if (prj->flag != WCS__QSC) {
4089 if (astQSCset(prj)) return 1;
4090 }
4091
4092 if (fabs(theta) == 90.0) {
4093 *x = 0.0;
4094 *y = copysign(2.0*prj->w[0],theta);
4095 return 0;
4096 }
4097
4098 cthe = astCosd(theta);
4099 l = cthe*astCosd(phi);
4100 m = cthe*astSind(phi);
4101 n = astSind(theta);
4102
4103 face = 0;
4104 rho = n;
4105 if (l > rho) {
4106 face = 1;
4107 rho = l;
4108 }
4109 if (m > rho) {
4110 face = 2;
4111 rho = m;
4112 }
4113 if (-l > rho) {
4114 face = 3;
4115 rho = -l;
4116 }
4117 if (-m > rho) {
4118 face = 4;
4119 rho = -m;
4120 }
4121 if (-n > rho) {
4122 face = 5;
4123 rho = -n;
4124 }
4125
4126 rhu = 1.0 - rho;
4127
4128 if (face == 0) {
4129 xi = m;
4130 eta = -l;
4131 if (rhu < 1.0e-8) {
4132 /* Small angle formula. */
4133 t = (90.0 - theta)*D2R;
4134 rhu = t*t/2.0;
4135 }
4136 x0 = 0.0;
4137 y0 = 2.0;
4138 } else if (face == 1) {
4139 xi = m;
4140 eta = n;
4141 if (rhu < 1.0e-8) {
4142 /* Small angle formula. */
4143 t = theta*D2R;
4144 p = fmod(phi,360.0);
4145 if (p < -180.0) p += 360.0;
4146 if (p > 180.0) p -= 360.0;
4147 p *= D2R;
4148 rhu = (p*p + t*t)/2.0;
4149 }
4150 x0 = 0.0;
4151 y0 = 0.0;
4152 } else if (face == 2) {
4153 xi = -l;
4154 eta = n;
4155 if (rhu < 1.0e-8) {
4156 /* Small angle formula. */
4157 t = theta*D2R;
4158 p = fmod(phi,360.0);
4159 if (p < -180.0) p += 360.0;
4160 p = (90.0 - p)*D2R;
4161 rhu = (p*p + t*t)/2.0;
4162 }
4163 x0 = 2.0;
4164 y0 = 0.0;
4165 } else if (face == 3) {
4166 xi = -m;
4167 eta = n;
4168 if (rhu < 1.0e-8) {
4169 /* Small angle formula. */
4170 t = theta*D2R;
4171 p = fmod(phi,360.0);
4172 if (p < 0.0) p += 360.0;
4173 p = (180.0 - p)*D2R;
4174 rhu = (p*p + t*t)/2.0;
4175 }
4176 x0 = 4.0;
4177 y0 = 0.0;
4178 } else if (face == 4) {
4179 xi = l;
4180 eta = n;
4181 if (rhu < 1.0e-8) {
4182 /* Small angle formula. */
4183 t = theta*D2R;
4184 p = fmod(phi,360.0);
4185 if (p > 180.0) p -= 360.0;
4186 p *= (90.0 + p)*D2R;
4187 rhu = (p*p + t*t)/2.0;
4188 }
4189 x0 = 6;
4190 y0 = 0.0;
4191 } else if (face == 5) {
4192 xi = m;
4193 eta = l;
4194 if (rhu < 1.0e-8) {
4195 /* Small angle formula. */
4196 t = (90.0 + theta)*D2R;
4197 rhu = t*t/2.0;
4198 }
4199 x0 = 0.0;
4200 y0 = -2;
4201 }
4202
4203 if (xi == 0.0 && eta == 0.0) {
4204 xf = 0.0;
4205 yf = 0.0;
4206 } else if (-xi >= fabs(eta)) {
4207 omega = eta/xi;
4208 tau = 1.0 + omega*omega;
4209 xf = -sqrt(rhu/(1.0-1.0/sqrt(1.0+tau)));
4210 yf = (xf/15.0)*(astATand(omega) - astASind(omega/sqrt(tau+tau)));
4211 } else if (xi >= fabs(eta)) {
4212 omega = eta/xi;
4213 tau = 1.0 + omega*omega;
4214 xf = sqrt(rhu/(1.0-1.0/sqrt(1.0+tau)));
4215 yf = (xf/15.0)*(astATand(omega) - astASind(omega/sqrt(tau+tau)));
4216 } else if (-eta > fabs(xi)) {
4217 omega = xi/eta;
4218 tau = 1.0 + omega*omega;
4219 yf = -sqrt(rhu/(1.0-1.0/sqrt(1.0+tau)));
4220 xf = (yf/15.0)*(astATand(omega) - astASind(omega/sqrt(tau+tau)));
4221 } else if (eta > fabs(xi)) {
4222 omega = xi/eta;
4223 tau = 1.0 + omega*omega;
4224 yf = sqrt(rhu/(1.0-1.0/sqrt(1.0+tau)));
4225 xf = (yf/15.0)*(astATand(omega) - astASind(omega/sqrt(tau+tau)));
4226 }
4227
4228 if (fabs(xf) > 1.0) {
4229 if (fabs(xf) > 1.0+tol) {
4230 return 2;
4231 }
4232 xf = copysign(1.0,xf);
4233 }
4234 if (fabs(yf) > 1.0) {
4235 if (fabs(yf) > 1.0+tol) {
4236 return 2;
4237 }
4238 yf = copysign(1.0,yf);
4239 }
4240
4241 *x = prj->w[0]*(xf + x0);
4242 *y = prj->w[0]*(yf + y0);
4243
4244
4245 return 0;
4246 }
4247
4248 /*--------------------------------------------------------------------------*/
4249
astQSCrev(x,y,prj,phi,theta)4250 int astQSCrev(x, y, prj, phi, theta)
4251
4252 const double x, y;
4253 struct AstPrjPrm *prj;
4254 double *phi, *theta;
4255
4256 {
4257 int direct, face;
4258 double l, m, n, omega, rho, rhu, tau, xf, yf, w;
4259 const double tol = 1.0e-12;
4260
4261 l = 0.0;
4262 m = 0.0;
4263 n = 0.0;
4264
4265 if (prj->flag != WCS__QSC) {
4266 if (astQSCset(prj)) return 1;
4267 }
4268
4269 xf = x*prj->w[1];
4270 yf = y*prj->w[1];
4271
4272 /* Check bounds. */
4273 if (fabs(xf) <= 1.0) {
4274 if (fabs(yf) > 3.0) return 2;
4275 } else {
4276 if (fabs(xf) > 7.0) return 2;
4277 if (fabs(yf) > 1.0) return 2;
4278 }
4279
4280 /* Map negative faces to the other side. */
4281 if (xf < -1.0) xf += 8.0;
4282
4283 /* Determine the face. */
4284 if (xf > 5.0) {
4285 face = 4;
4286 xf = xf - 6.0;
4287 } else if (xf > 3.0) {
4288 face = 3;
4289 xf = xf - 4.0;
4290 } else if (xf > 1.0) {
4291 face = 2;
4292 xf = xf - 2.0;
4293 } else if (yf > 1.0) {
4294 face = 0;
4295 yf = yf - 2.0;
4296 } else if (yf < -1.0) {
4297 face = 5;
4298 yf = yf + 2.0;
4299 } else {
4300 face = 1;
4301 }
4302
4303 direct = (fabs(xf) > fabs(yf));
4304 if (direct) {
4305 if (xf == 0.0) {
4306 omega = 0.0;
4307 tau = 1.0;
4308 rho = 1.0;
4309 rhu = 0.0;
4310 } else {
4311 w = 15.0*yf/xf;
4312 omega = astSind(w)/(astCosd(w) - SQRT2INV);
4313 tau = 1.0 + omega*omega;
4314 rhu = xf*xf*(1.0 - 1.0/sqrt(1.0 + tau));
4315 rho = 1.0 - rhu;
4316 }
4317 } else {
4318 if (yf == 0.0) {
4319 omega = 0.0;
4320 tau = 1.0;
4321 rho = 1.0;
4322 rhu = 0.0;
4323 } else {
4324 w = 15.0*xf/yf;
4325 omega = astSind(w)/(astCosd(w) - SQRT2INV);
4326 tau = 1.0 + omega*omega;
4327 rhu = yf*yf*(1.0 - 1.0/sqrt(1.0 + tau));
4328 rho = 1.0 - rhu;
4329 }
4330 }
4331
4332 if (rho < -1.0) {
4333 if (rho < -1.0-tol) {
4334 return 2;
4335 }
4336
4337 rho = -1.0;
4338 rhu = 2.0;
4339 w = 0.0;
4340 } else {
4341 w = sqrt(rhu*(2.0-rhu)/tau);
4342 }
4343
4344 if (face == 0) {
4345 n = rho;
4346 if (direct) {
4347 m = w;
4348 if (xf < 0.0) m = -m;
4349 l = -m*omega;
4350 } else {
4351 l = w;
4352 if (yf > 0.0) l = -l;
4353 m = -l*omega;
4354 }
4355 } else if (face == 1) {
4356 l = rho;
4357 if (direct) {
4358 m = w;
4359 if (xf < 0.0) m = -m;
4360 n = m*omega;
4361 } else {
4362 n = w;
4363 if (yf < 0.0) n = -n;
4364 m = n*omega;
4365 }
4366 } else if (face == 2) {
4367 m = rho;
4368 if (direct) {
4369 l = w;
4370 if (xf > 0.0) l = -l;
4371 n = -l*omega;
4372 } else {
4373 n = w;
4374 if (yf < 0.0) n = -n;
4375 l = -n*omega;
4376 }
4377 } else if (face == 3) {
4378 l = -rho;
4379 if (direct) {
4380 m = w;
4381 if (xf > 0.0) m = -m;
4382 n = -m*omega;
4383 } else {
4384 n = w;
4385 if (yf < 0.0) n = -n;
4386 m = -n*omega;
4387 }
4388 } else if (face == 4) {
4389 m = -rho;
4390 if (direct) {
4391 l = w;
4392 if (xf < 0.0) l = -l;
4393 n = l*omega;
4394 } else {
4395 n = w;
4396 if (yf < 0.0) n = -n;
4397 l = n*omega;
4398 }
4399 } else if (face == 5) {
4400 n = -rho;
4401 if (direct) {
4402 m = w;
4403 if (xf < 0.0) m = -m;
4404 l = m*omega;
4405 } else {
4406 l = w;
4407 if (yf < 0.0) l = -l;
4408 m = l*omega;
4409 }
4410 }
4411
4412 if (l == 0.0 && m == 0.0) {
4413 *phi = 0.0;
4414 } else {
4415 *phi = astATan2d(m, l);
4416 }
4417 *theta = astASind(n);
4418
4419 return 0;
4420 }
4421
4422 /*============================================================================
4423 * HPX: HEALPix projection.
4424 *
4425 * Given:
4426 * prj->p[1] H - the number of facets in longitude.
4427 * prj->p[2] K - the number of facets in latitude
4428 *
4429 * Given and/or returned:
4430 * prj->r0 Reset to 180/pi if 0.
4431 * prj->phi0 Reset to 0.0
4432 * prj->theta0 Reset to 0.0
4433 *
4434 * Returned:
4435 * prj->flag HPX
4436 * prj->code "HPX"
4437 * prj->n True if K is odd.
4438 * prj->w[0] r0*(pi/180)
4439 * prj->w[1] (180/pi)/r0
4440 * prj->w[2] (K-1)/K
4441 * prj->w[3] 90*K/H
4442 * prj->w[4] (K+1)/2
4443 * prj->w[5] 90*(K-1)/H
4444 * prj->w[6] 180/H
4445 * prj->w[7] H/360
4446 * prj->w[8] (90*K/H)*r0*(pi/180)
4447 * prj->w[9] (180/H)*r0*(pi/180)
4448 * prj->astPRJfwd Pointer to astHPXfwd().
4449 * prj->astPRJrev Pointer to astHPXrev().
4450
4451
4452 *===========================================================================*/
4453
astHPXset(prj)4454 int astHPXset(prj)
4455
4456 struct AstPrjPrm *prj;
4457
4458 {
4459 strcpy(prj->code, "HPX");
4460 prj->flag = WCS__HPX;
4461 prj->phi0 = 0.0;
4462 prj->theta0 = 0.0;
4463
4464 prj->n = ((int)prj->p[2])%2;
4465
4466 if (prj->r0 == 0.0) {
4467 prj->r0 = R2D;
4468 prj->w[0] = 1.0;
4469 prj->w[1] = 1.0;
4470 } else {
4471 prj->w[0] = prj->r0*D2R;
4472 prj->w[1] = R2D/prj->r0;
4473 }
4474
4475 prj->w[2] = (prj->p[2] - 1.0) / prj->p[2];
4476 prj->w[3] = 90.0 * prj->p[2] / prj->p[1];
4477 prj->w[4] = (prj->p[2] + 1.0) / 2.0;
4478 prj->w[5] = 90.0 * (prj->p[2] - 1.0) / prj->p[1];
4479 prj->w[6] = 180.0 / prj->p[1];
4480 prj->w[7] = prj->p[1] / 360.0;
4481 prj->w[8] = prj->w[3] * prj->w[0];
4482 prj->w[9] = prj->w[6] * prj->w[0];
4483
4484 prj->astPRJfwd = astHPXfwd;
4485 prj->astPRJrev = astHPXrev;
4486
4487 return 0;
4488 }
4489
4490 /*--------------------------------------------------------------------------*/
4491
astHPXfwd(phi,theta,prj,x,y)4492 int astHPXfwd(phi, theta, prj, x, y)
4493
4494 const double phi, theta;
4495 struct AstPrjPrm *prj;
4496 double *x, *y;
4497
4498 {
4499 double abssin, sigma, sinthe, phic;
4500 int hodd;
4501
4502 if( prj->flag != WCS__HPX ) {
4503 if( astHPXset( prj ) ) return 1;
4504 }
4505
4506 sinthe = astSind( theta );
4507 abssin = fabs( sinthe );
4508
4509 /* Equatorial zone */
4510 if( abssin <= prj->w[2] ) {
4511 *x = prj->w[0] * phi;
4512 *y = prj->w[8] * sinthe;
4513
4514 /* Polar zone */
4515 } else {
4516
4517 /* DSB - The expression for phic is conditioned differently to the
4518 WCSLIB code in order to improve accuracy of the floor function for
4519 arguments very slightly below an integer value. */
4520 hodd = ((int)prj->p[1]) % 2;
4521 if( !prj->n && theta <= 0.0 ) hodd = 1 - hodd;
4522 if( hodd ) {
4523 phic = -180.0 + (2.0*floor( prj->w[7] * phi + 1/2 ) + prj->p[1] ) * prj->w[6];
4524 } else {
4525 phic = -180.0 + (2.0*floor( prj->w[7] * phi ) + prj->p[1] + 1 ) * prj->w[6];
4526 }
4527
4528 sigma = sqrt( prj->p[2]*( 1.0 - abssin ));
4529
4530 *x = prj->w[0] *( phic + ( phi - phic )*sigma );
4531
4532 *y = prj->w[9] * ( prj->w[4] - sigma );
4533 if( theta < 0 ) *y = -*y;
4534
4535 }
4536
4537 return 0;
4538 }
4539
4540 /*--------------------------------------------------------------------------*/
4541
astHPXrev(x,y,prj,phi,theta)4542 int astHPXrev(x, y, prj, phi, theta)
4543
4544 const double x, y;
4545 struct AstPrjPrm *prj;
4546 double *phi, *theta;
4547
4548 {
4549 double absy, sigma, t, yr, xc;
4550 int hodd;
4551
4552 if (prj->flag != WCS__HPX) {
4553 if (astHPXset(prj)) return 1;
4554 }
4555
4556 yr = prj->w[1]*y;
4557 absy = fabs( yr );
4558
4559 /* Equatorial zone */
4560 if( absy <= prj->w[5] ) {
4561 *phi = prj->w[1] * x;
4562 t = yr/prj->w[3];
4563 if( t < -1.0 || t > 1.0 ) {
4564 return 2;
4565 } else {
4566 *theta = astASind( t );
4567 }
4568
4569 /* Polar zone */
4570 } else if( absy <= 90 ){
4571
4572
4573 /* DSB - The expression for xc is conditioned differently to the
4574 WCSLIB code in order to improve accuracy of the floor function for
4575 arguments very slightly below an integer value. */
4576 hodd = ((int)prj->p[1]) % 2;
4577 if( !prj->n && yr <= 0.0 ) hodd = 1 - hodd;
4578 if( hodd ) {
4579 xc = -180.0 + (2.0*floor( prj->w[7] * x + 1/2 ) + prj->p[1] ) * prj->w[6];
4580 } else {
4581 xc = -180.0 + (2.0*floor( prj->w[7] * x ) + prj->p[1] + 1 ) * prj->w[6];
4582 }
4583
4584 sigma = prj->w[4] - absy / prj->w[6];
4585
4586 if( sigma == 0.0 ) {
4587 return 2;
4588 } else {
4589
4590 t = ( x - xc )/sigma;
4591 if( fabs( t ) <= prj->w[6] ) {
4592 *phi = prj->w[1] *( xc + t );
4593 } else {
4594 return 2;
4595 }
4596 }
4597
4598 t = 1.0 - sigma*sigma/prj->p[2];
4599 if( t < -1.0 || t > 1.0 ) {
4600 return 2;
4601 } else {
4602 *theta = astASind( t );
4603 if( y < 0 ) *theta = -*theta;
4604 }
4605
4606 } else {
4607 return 2;
4608 }
4609
4610 return 0;
4611 }
4612
4613 /*============================================================================
4614 * XPH: HEALPix polar, aka "butterfly" projection.
4615 *
4616 * Given and/or returned:
4617 * prj->r0 Reset to 180/pi if 0.
4618 * prj->phi0 Reset to 0.0 if undefined.
4619 * prj->theta0 Reset to 0.0 if undefined.
4620 *
4621 * Returned:
4622 * prj->flag XPH
4623 * prj->code "XPH"
4624 * prj->w[0] r0*(pi/180)/sqrt(2)
4625 * prj->w[1] (180/pi)/r0/sqrt(2)
4626 * prj->w[2] 2/3
4627 * prj->w[3] tol (= 1e-4)
4628 * prj->w[4] sqrt(2/3)*(180/pi)
4629 * prj->w[5] 90 - tol*sqrt(2/3)*(180/pi)
4630 * prj->w[6] sqrt(3/2)*(pi/180)
4631 * prj->astPRJfwd Pointer to astXPHfwd().
4632 * prj->astPRJrev Pointer to astXPHrev().
4633 *===========================================================================*/
4634
astXPHset(prj)4635 int astXPHset(prj)
4636
4637 struct AstPrjPrm *prj;
4638
4639 {
4640 strcpy(prj->code, "XPH");
4641 prj->flag = WCS__XPH;
4642
4643 if (prj->r0 == 0.0) {
4644 prj->r0 = R2D;
4645 prj->w[0] = 1.0;
4646 prj->w[1] = 1.0;
4647 } else {
4648 prj->w[0] = prj->r0*D2R;
4649 prj->w[1] = R2D/prj->r0;
4650 }
4651
4652 prj->w[0] /= sqrt(2.0);
4653 prj->w[1] /= sqrt(2.0);
4654 prj->w[2] = 2.0/3.0;
4655 prj->w[3] = 1e-4;
4656 prj->w[4] = sqrt(prj->w[2])*R2D;
4657 prj->w[5] = 90.0 - prj->w[3]*prj->w[4];
4658 prj->w[6] = sqrt(1.5)*D2R;
4659
4660 prj->astPRJfwd = astXPHfwd;
4661 prj->astPRJrev = astXPHrev;
4662
4663 return 0;
4664 }
4665
4666 /*--------------------------------------------------------------------------*/
4667
astXPHfwd(phi,theta,prj,x,y)4668 int astXPHfwd(phi, theta, prj, x, y)
4669
4670 const double phi, theta;
4671 struct AstPrjPrm *prj;
4672 double *x, *y;
4673
4674 {
4675 double abssin, chi, eta, psi, sigma, sinthe, xi;
4676
4677 if (prj->flag != WCS__XPH) {
4678 if (astXPHset(prj)) return 1;
4679 }
4680
4681 /* Do phi dependence. */
4682 chi = phi;
4683 if (180.0 <= fabs(chi)) {
4684 chi = fmod(chi, 360.0);
4685 if (chi < -180.0) {
4686 chi += 360.0;
4687 } else if (180.0 <= chi) {
4688 chi -= 360.0;
4689 }
4690 }
4691
4692 /* phi is also recomputed from chi to avoid rounding problems. */
4693 chi += 180.0;
4694 psi = fmod(chi, 90.0);
4695
4696 /* y is used to hold phi (rounded). */
4697 *x = psi;
4698 *y = chi - 180.0;
4699
4700 /* Do theta dependence. */
4701 sinthe = astSind(theta);
4702 abssin = fabs(sinthe);
4703
4704 if (abssin <= prj->w[2]) {
4705 /* Equatorial regime. */
4706 xi = *x;
4707 eta = 67.5 * sinthe;
4708
4709 } else {
4710 /* Polar regime. */
4711 if (theta < prj->w[5]) {
4712 sigma = sqrt(3.0*(1.0 - abssin));
4713 } else {
4714 sigma = (90.0 - theta)*prj->w[6];
4715 }
4716
4717 xi = 45.0 + (*x - 45.0)*sigma;
4718 eta = 45.0 * (2.0 - sigma);
4719 if (theta < 0.0) eta = -eta;
4720 }
4721
4722 xi -= 45.0;
4723 eta -= 90.0;
4724
4725 /* Recall that y holds phi. */
4726 if (*y < -90.0) {
4727 *x = prj->w[0]*(-xi + eta);
4728 *y = prj->w[0]*(-xi - eta);
4729
4730 } else if (*y < 0.0) {
4731 *x = prj->w[0]*(+xi + eta);
4732 *y = prj->w[0]*(-xi + eta);
4733
4734 } else if (*y < 90.0) {
4735 *x = prj->w[0]*( xi - eta);
4736 *y = prj->w[0]*( xi + eta);
4737
4738 } else {
4739 *x = prj->w[0]*(-xi - eta);
4740 *y = prj->w[0]*( xi - eta);
4741 }
4742
4743 return 0;
4744
4745 }
4746
4747 /*--------------------------------------------------------------------------*/
4748
astXPHrev(x,y,prj,phi,theta)4749 int astXPHrev(x, y, prj, phi, theta)
4750
4751 const double x, y;
4752 struct AstPrjPrm *prj;
4753 double *phi, *theta;
4754
4755 {
4756 double abseta, eta, eta1, sigma, xi, xi1, xr, yr;
4757 const double tol = 1.0e-12;
4758
4759 if (prj->flag != WCS__XPH) {
4760 if (astXPHset(prj)) return 1;
4761 }
4762
4763
4764 xr = x*prj->w[1];
4765 yr = y*prj->w[1];
4766 if (xr <= 0.0 && 0.0 < yr) {
4767 xi1 = -xr - yr;
4768 eta1 = xr - yr;
4769 *phi = -180.0;
4770 } else if (xr < 0.0 && yr <= 0.0) {
4771 xi1 = xr - yr;
4772 eta1 = xr + yr;
4773 *phi = -90.0;
4774 } else if (0.0 <= xr && yr < 0.0) {
4775 xi1 = xr + yr;
4776 eta1 = -xr + yr;
4777 *phi = 0.0;
4778 } else {
4779 xi1 = -xr + yr;
4780 eta1 = -xr - yr;
4781 *phi = 90.0;
4782 }
4783
4784 xi = xi1 + 45.0;
4785 eta = eta1 + 90.0;
4786 abseta = fabs(eta);
4787
4788 if (abseta <= 90.0) {
4789 if (abseta <= 45.0) {
4790 /* Equatorial regime. */
4791 *phi += xi;
4792 *theta = astASind(eta/67.5);
4793
4794 /* Bounds checking. */
4795 if (45.0+tol < fabs(xi1)) return 2;
4796
4797 } else {
4798 /* Polar regime. */
4799 sigma = (90.0 - abseta) / 45.0;
4800
4801 /* Ensure an exact result for points on the boundary. */
4802 if (xr == 0.0) {
4803 if (yr <= 0.0) {
4804 *phi = 0.0;
4805 } else {
4806 *phi = 180.0;
4807 }
4808 } else if (yr == 0.0) {
4809 if (xr < 0.0) {
4810 *phi = -90.0;
4811 } else {
4812 *phi = 90.0;
4813 }
4814 } else {
4815 *phi += 45.0 + xi1/sigma;
4816 }
4817
4818 if (sigma < prj->w[3]) {
4819 *theta = 90.0 - sigma*prj->w[4];
4820 } else {
4821 *theta = astASind(1.0 - sigma*sigma/3.0);
4822 }
4823 if (eta < 0.0) *theta = -(*theta);
4824
4825 /* Bounds checking. */
4826 if (eta < -45.0 && eta+90.0+tol < fabs(xi1)) return 2;
4827 }
4828
4829 } else {
4830 /* Beyond latitude range. */
4831 *phi = 0.0;
4832 *theta = 0.0;
4833 return 2;
4834 }
4835
4836 return 0;
4837 }
4838
4839
4840