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
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2 of the License, or (at your option) any later version.
10 *
11 * This library is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with this library; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 * Correspondence concerning WCSLIB may be directed to:
21 * Internet email: mcalabre@atnf.csiro.au
22 * Postal address: Dr. Mark Calabretta,
23 * Australia Telescope National Facility,
24 * P.O. Box 76,
25 * Epping, NSW, 2121,
26 * AUSTRALIA
27 *
28 *=============================================================================
29 *
30 * C routines which implement the FITS World Coordinate System (WCS)
31 * convention.
32 *
33 * Summary of routines
34 * -------------------
35 * These routines are provided as drivers for the lower level spherical
36 * coordinate transformation and projection routines. There are separate
37 * driver routines for the forward, celfwd(), and reverse, celrev(),
38 * transformations.
39 *
40 * An initialization routine, celset(), computes intermediate values from
41 * the transformation parameters but need not be called explicitly - see the
42 * explanation of cel.flag below.
43 *
44 *
45 * Initialization routine; celset()
46 * --------------------------------
47 * Initializes members of a celprm data structure which hold intermediate
48 * values. Note that this routine need not be called directly; it will be
49 * invoked by celfwd() and celrev() if the "flag" structure member is
50 * anything other than a predefined magic value.
51 *
52 * Given:
53 * pcode[4] const char
54 * WCS projection code (see below).
55 *
56 * Given and returned:
57 * cel celprm* Spherical coordinate transformation parameters
58 * (see below).
59 * prj prjprm* Projection parameters (usage is described in the
60 * prologue to "proj.c").
61 *
62 * Function return value:
63 * int Error status
64 * 0: Success.
65 * 1: Invalid coordinate transformation parameters.
66 * 2: Ill-conditioned coordinate transformation
67 * parameters.
68 *
69 * Forward transformation; celfwd()
70 * --------------------------------
71 * Compute (x,y) coordinates in the plane of projection from celestial
72 * coordinates (lng,lat).
73 *
74 * Given:
75 * pcode[4] const char
76 * WCS projection code (see below).
77 * lng,lat const double
78 * Celestial longitude and latitude of the projected
79 * point, in degrees.
80 *
81 * Given and returned:
82 * cel celprm* Spherical coordinate transformation parameters
83 * (see below).
84 *
85 * Returned:
86 * phi, double* Longitude and latitude in the native coordinate
87 * theta system of the projection, in degrees.
88 *
89 * Given and returned:
90 * prj prjprm* Projection parameters (usage is described in the
91 * prologue to "proj.c").
92 *
93 * Returned:
94 * x,y double* Projected coordinates, "degrees".
95 *
96 * Function return value:
97 * int Error status
98 * 0: Success.
99 * 1: Invalid coordinate transformation parameters.
100 * 2: Invalid projection parameters.
101 * 3: Invalid value of (lng,lat).
102 *
103 * Reverse transformation; celrev()
104 * --------------------------------
105 * Compute the celestial coordinates (lng,lat) of the point with projected
106 * coordinates (x,y).
107 *
108 * Given:
109 * pcode[4] const char
110 * WCS projection code (see below).
111 * x,y const double
112 * Projected coordinates, "degrees".
113 *
114 * Given and returned:
115 * prj prjprm* Projection parameters (usage is described in the
116 * prologue to "proj.c").
117 *
118 * Returned:
119 * phi, double* Longitude and latitude in the native coordinate
120 * theta system of the projection, in degrees.
121 *
122 * Given and returned:
123 * cel celprm* Spherical coordinate transformation parameters
124 * (see below).
125 *
126 * Returned:
127 * lng,lat double* Celestial longitude and latitude of the projected
128 * point, in degrees.
129 *
130 * Function return value:
131 * int Error status
132 * 0: Success.
133 * 1: Invalid coordinate transformation parameters.
134 * 2: Invalid projection parameters.
135 * 3: Invalid value of (x,y).
136 *
137 * Coordinate transformation parameters
138 * ------------------------------------
139 * The celprm struct consists of the following:
140 *
141 * int flag
142 * The celprm struct contains pointers to the forward and reverse
143 * projection routines as well as intermediaries computed from the
144 * reference coordinates (see below). Whenever the projection code
145 * (pcode) or any of ref[4] are set or changed then this flag must be
146 * set to zero to signal the initialization routine, celset(), to
147 * redetermine the function pointers and recompute intermediaries.
148 * Once this has been done pcode itself is ignored.
149 *
150 * double ref[4]
151 * The first pair of values should be set to the celestial longitude
152 * and latitude (usually right ascension and declination) of the
153 * reference point of the projection. These are given by the CRVALn
154 * keywords in FITS.
155 *
156 * The second pair of values are the native longitude of the celestial
157 * pole and the celestial latitude of the native pole and correspond to
158 * FITS keywords LONPOLE and LATPOLE.
159 *
160 * LONPOLE defaults to 0 degrees if the celestial latitude of the
161 * reference point of the projection is greater than the native
162 * latitude, otherwise 180 degrees. (This is the condition for the
163 * celestial latitude to increase in the same direction as the native
164 * latitude at the reference point.) ref[2] may be set to 999.0 to
165 * indicate that the correct default should be substituted.
166 *
167 * In some circumstances the celestial latitude of the native pole may
168 * be determined by the first three values only to within a sign and
169 * LATPOLE is used to choose between the two solutions. LATPOLE is
170 * set in ref[3] and the solution closest to this value is used to
171 * reset ref[3]. It is therefore legitimate, for example, to set
172 * ref[3] to 999.0 to choose the more northerly solution - the default
173 * if the LATPOLE card is omitted from the FITS header. For the
174 * special case where the reference point of the projection is at
175 * native latitude zero, its celestial latitude is zero, and
176 * LONPOLE = +/- 90 then the celestial latitude of the pole is not
177 * determined by the first three reference values and LATPOLE
178 * specifies it completely.
179 *
180 * The remaining members of the celprm struct are maintained by the
181 * initialization routines and should not be modified. This is done for the
182 * sake of efficiency and to allow an arbitrary number of contexts to be
183 * maintained simultaneously.
184 *
185 * double euler[5]
186 * Euler angles and associated intermediaries derived from the
187 * coordinate reference values.
188 *
189 *
190 * WCS projection codes
191 * --------------------
192 * Zenithals/azimuthals:
193 * AZP: zenithal/azimuthal perspective
194 * TAN: gnomonic
195 * STG: stereographic
196 * SIN: synthesis (generalized orthographic)
197 * ARC: zenithal/azimuthal equidistant
198 * ZPN: zenithal/azimuthal polynomial
199 * ZEA: zenithal/azimuthal equal area
200 * AIR: Airy
201 *
202 * Cylindricals:
203 * CYP: cylindrical perspective
204 * CEA: cylindrical equal area
205 * CAR: Cartesian
206 * MER: Mercator
207 *
208 * Pseudo-cylindricals:
209 * SFL: Sanson-Flamsteed
210 * PAR: parabolic
211 * MOL: Mollweide
212 *
213 * Conventional:
214 * AIT: Hammer-Aitoff
215 *
216 * Conics:
217 * COP: conic perspective
218 * COD: conic equidistant
219 * COE: conic equal area
220 * COO: conic orthomorphic
221 *
222 * Polyconics:
223 * BON: Bonne
224 * PCO: polyconic
225 *
226 * Quad-cubes:
227 * TSC: tangential spherical cube
228 * CSC: COBE quadrilateralized spherical cube
229 * QSC: quadrilateralized spherical cube
230 *
231 * Author: Mark Calabretta, Australia Telescope National Facility
232 * $Id: cel.c,v 2.14 2002/04/03 01:25:29 mcalabre Exp $
233 *===========================================================================*/
234
235 #include <math.h>
236 #include <string.h>
237 #include "wcslib.h"
238
239 /* Map error number to error message for each function. */
240 const char *celset_errmsg[] = {
241 0,
242 "Invalid coordinate transformation parameters",
243 "Ill-conditioned coordinate transformation parameters"};
244
245 const char *celfwd_errmsg[] = {
246 0,
247 "Invalid coordinate transformation parameters",
248 "Invalid projection parameters",
249 "Invalid value of (lng,lat)"};
250
251 const char *celrev_errmsg[] = {
252 0,
253 "Invalid coordinate transformation parameters",
254 "Invalid projection parameters",
255 "Invalid value of (x,y)"};
256
257
258 int
celset(pcode,cel,prj)259 celset(pcode, cel, prj)
260
261 const char pcode[4];
262 struct celprm *cel;
263 struct prjprm *prj;
264
265 {
266 int dophip;
267 const double tol = 1.0e-10;
268 double clat0, cphip, cthe0, slat0, sphip, sthe0;
269 double latp, latp1, latp2;
270 double u, v, x, y, z;
271
272 /* Initialize the projection driver routines. */
273 if (prjset(pcode, prj)) {
274 return 1;
275 }
276
277 /* Set default for native longitude of the celestial pole? */
278 dophip = (cel->ref[2] == 999.0);
279
280 /* Compute celestial coordinates of the native pole. */
281 if (prj->theta0 == 90.0) {
282 /* Reference point is at the native pole. */
283
284 if (dophip) {
285 /* Set default for longitude of the celestial pole. */
286 cel->ref[2] = 180.0;
287 }
288
289 latp = cel->ref[1];
290 cel->ref[3] = latp;
291
292 cel->euler[0] = cel->ref[0];
293 cel->euler[1] = 90.0 - latp;
294 } else {
295 /* Reference point away from the native pole. */
296
297 /* Set default for longitude of the celestial pole. */
298 if (dophip) {
299 cel->ref[2] = (cel->ref[1] < prj->theta0) ? 180.0 : 0.0;
300 }
301
302 clat0 = cosdeg (cel->ref[1]);
303 slat0 = sindeg (cel->ref[1]);
304 cphip = cosdeg (cel->ref[2]);
305 sphip = sindeg (cel->ref[2]);
306 cthe0 = cosdeg (prj->theta0);
307 sthe0 = sindeg (prj->theta0);
308
309 x = cthe0*cphip;
310 y = sthe0;
311 z = sqrt(x*x + y*y);
312 if (z == 0.0) {
313 if (slat0 != 0.0) {
314 return 1;
315 }
316
317 /* latp determined by LATPOLE in this case. */
318 latp = cel->ref[3];
319 } else {
320 if (fabs(slat0/z) > 1.0) {
321 return 1;
322 }
323
324 u = atan2deg (y,x);
325 v = acosdeg (slat0/z);
326
327 latp1 = u + v;
328 if (latp1 > 180.0) {
329 latp1 -= 360.0;
330 } else if (latp1 < -180.0) {
331 latp1 += 360.0;
332 }
333
334 latp2 = u - v;
335 if (latp2 > 180.0) {
336 latp2 -= 360.0;
337 } else if (latp2 < -180.0) {
338 latp2 += 360.0;
339 }
340
341 if (fabs(cel->ref[3]-latp1) < fabs(cel->ref[3]-latp2)) {
342 if (fabs(latp1) < 90.0+tol) {
343 latp = latp1;
344 } else {
345 latp = latp2;
346 }
347 } else {
348 if (fabs(latp2) < 90.0+tol) {
349 latp = latp2;
350 } else {
351 latp = latp1;
352 }
353 }
354
355 cel->ref[3] = latp;
356 }
357
358 cel->euler[1] = 90.0 - latp;
359
360 z = cosdeg (latp)*clat0;
361 if (fabs(z) < tol) {
362 if (fabs(clat0) < tol) {
363 /* Celestial pole at the reference point. */
364 cel->euler[0] = cel->ref[0];
365 cel->euler[1] = 90.0 - prj->theta0;
366 } else if (latp > 0.0) {
367 /* Celestial pole at the native north pole.*/
368 cel->euler[0] = cel->ref[0] + cel->ref[2] - 180.0;
369 cel->euler[1] = 0.0;
370 } else if (latp < 0.0) {
371 /* Celestial pole at the native south pole. */
372 cel->euler[0] = cel->ref[0] - cel->ref[2];
373 cel->euler[1] = 180.0;
374 }
375 } else {
376 x = (sthe0 - sindeg (latp)*slat0)/z;
377 y = sphip*cthe0/clat0;
378 if (x == 0.0 && y == 0.0) {
379 return 1;
380 }
381 cel->euler[0] = cel->ref[0] - atan2deg (y,x);
382 }
383
384 /* Make euler[0] the same sign as ref[0]. */
385 if (cel->ref[0] >= 0.0) {
386 if (cel->euler[0] < 0.0) cel->euler[0] += 360.0;
387 } else {
388 if (cel->euler[0] > 0.0) cel->euler[0] -= 360.0;
389 }
390 }
391
392 cel->euler[2] = cel->ref[2];
393 cel->euler[3] = cosdeg (cel->euler[1]);
394 cel->euler[4] = sindeg (cel->euler[1]);
395 cel->flag = CELSET;
396
397 /* Check for ill-conditioned parameters. */
398 if (fabs(latp) > 90.0+tol) {
399 return 2;
400 }
401
402 return 0;
403 }
404
405 /*--------------------------------------------------------------------------*/
406
407 int
celfwd(pcode,lng,lat,cel,phi,theta,prj,x,y)408 celfwd(pcode, lng, lat, cel, phi, theta, prj, x, y)
409
410 const char pcode[4];
411 const double lng, lat;
412 struct celprm *cel;
413 double *phi, *theta;
414 struct prjprm *prj;
415 double *x, *y;
416
417 {
418 int err;
419
420 if (cel->flag != CELSET) {
421 if (celset(pcode, cel, prj)) return 1;
422 }
423
424 /* Compute native coordinates. */
425 sphfwd(lng, lat, cel->euler, phi, theta);
426
427 /* Apply forward projection. */
428 if ((err = prj->prjfwd(*phi, *theta, prj, x, y))) {
429 return err == 1 ? 2 : 3;
430 }
431
432 return 0;
433 }
434
435 /*--------------------------------------------------------------------------*/
436
437 int
celrev(pcode,x,y,prj,phi,theta,cel,lng,lat)438 celrev(pcode, x, y, prj, phi, theta, cel, lng, lat)
439
440 const char pcode[4];
441 const double x, y;
442 struct prjprm *prj;
443 double *phi, *theta;
444 struct celprm *cel;
445 double *lng, *lat;
446
447 {
448 int err;
449
450 if (cel->flag != CELSET) {
451 if(celset(pcode, cel, prj)) return 1;
452 }
453
454 /* Apply reverse projection. */
455 if ((err = prj->prjrev(x, y, prj, phi, theta))) {
456 return err == 1 ? 2 : 3;
457 }
458
459 /* Compute native coordinates. */
460 sphrev(*phi, *theta, cel->euler, lng, lat);
461
462 return 0;
463 }
464
465 /* Dec 20 1999 Doug Mink - Change cosd() and sind() to cosdeg() and sindeg()
466 * Dec 20 1999 Doug Mink - Include wcslib.h, which includes wcsmath.h and cel.h
467 *
468 * Dec 18 2000 Doug Mink - Include string.h for strcmp()
469 *
470 * Mar 20 2001 Doug Mink - Add () around err assignments in if statements
471 * Sep 19 2001 Doug Mink - Add above changes to WCSLIB-2.7 cel.c
472 *
473 * Mar 12 2002 Doug Mink - Add changes to WCSLIB-2.8.2 cel.c
474 */
475