1 /*============================================================================
2 WCSLIB 7.7 - an implementation of the FITS WCS standard.
3 Copyright (C) 1995-2021, Mark Calabretta
4
5 This file is part of WCSLIB.
6
7 WCSLIB is free software: you can redistribute it and/or modify it under the
8 terms of the GNU Lesser General Public License as published by the Free
9 Software Foundation, either version 3 of the License, or (at your option)
10 any later version.
11
12 WCSLIB is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 more details.
16
17 You should have received a copy of the GNU Lesser General Public License
18 along with WCSLIB. If not, see http://www.gnu.org/licenses.
19
20 Author: Mark Calabretta, Australia Telescope National Facility, CSIRO.
21 http://www.atnf.csiro.au/people/Mark.Calabretta
22 $Id: cel.c,v 7.7 2021/07/12 06:36:49 mcalabre Exp $
23 *===========================================================================*/
24
25 #include <math.h>
26 #include <stdio.h>
27 #include <stdlib.h>
28
29 #include "wcserr.h"
30 #include "wcsmath.h"
31 #include "wcsprintf.h"
32 #include "wcstrig.h"
33 #include "sph.h"
34 #include "cel.h"
35
36 const int CELSET = 137;
37
38 // Map status return value to message.
39 const char *cel_errmsg[] = {
40 "Success",
41 "Null celprm pointer passed",
42 "Invalid projection parameters",
43 "Invalid coordinate transformation parameters",
44 "Ill-conditioned coordinate transformation parameters",
45 "One or more of the (x,y) coordinates were invalid",
46 "One or more of the (lng,lat) coordinates were invalid"};
47
48 // Map error returns for lower-level routines.
49 const int cel_prjerr[] = {
50 CELERR_SUCCESS, // 0: PRJERR_SUCCESS
51 CELERR_NULL_POINTER, // 1: PRJERR_NULL_POINTER
52 CELERR_BAD_PARAM, // 2: PRJERR_BAD_PARAM
53 CELERR_BAD_PIX, // 3: PRJERR_BAD_PIX
54 CELERR_BAD_WORLD // 4: PRJERR_BAD_WORLD
55 };
56
57 // Convenience macro for invoking wcserr_set().
58 #define CEL_ERRMSG(status) WCSERR_SET(status), cel_errmsg[status]
59
60 //----------------------------------------------------------------------------
61
celini(struct celprm * cel)62 int celini(struct celprm *cel)
63
64 {
65 register int k;
66
67 if (cel == 0x0) return CELERR_NULL_POINTER;
68
69 cel->flag = 0;
70
71 cel->offset = 0;
72 cel->phi0 = UNDEFINED;
73 cel->theta0 = UNDEFINED;
74 cel->ref[0] = 0.0;
75 cel->ref[1] = 0.0;
76 cel->ref[2] = UNDEFINED;
77 cel->ref[3] = +90.0;
78
79 for (k = 0; k < 5; cel->euler[k++] = 0.0);
80 cel->latpreq = -1;
81
82 cel->err = 0x0;
83
84 return cel_prjerr[prjini(&(cel->prj))];
85 }
86
87 //----------------------------------------------------------------------------
88
celfree(struct celprm * cel)89 int celfree(struct celprm *cel)
90
91 {
92 if (cel == 0x0) return CELERR_NULL_POINTER;
93
94 wcserr_clear(&(cel->err));
95
96 return cel_prjerr[prjfree(&(cel->prj))];
97 }
98
99 //----------------------------------------------------------------------------
100
celsize(const struct celprm * cel,int sizes[2])101 int celsize(const struct celprm *cel, int sizes[2])
102
103 {
104 if (cel == 0x0) {
105 sizes[0] = sizes[1] = 0;
106 return CELERR_SUCCESS;
107 }
108
109 // Base size, in bytes.
110 sizes[0] = sizeof(struct celprm);
111
112 // Total size of allocated memory, in bytes.
113 sizes[1] = 0;
114
115 int exsizes[2];
116
117 // celprm::prj.
118 prjsize(&(cel->prj), exsizes);
119 sizes[1] += exsizes[1];
120
121 // celprm::err.
122 wcserr_size(cel->err, exsizes);
123 sizes[1] += exsizes[0] + exsizes[1];
124
125 return CELERR_SUCCESS;
126 }
127
128 //----------------------------------------------------------------------------
129
celprt(const struct celprm * cel)130 int celprt(const struct celprm *cel)
131
132 {
133 int i;
134
135 if (cel == 0x0) return CELERR_NULL_POINTER;
136
137 wcsprintf(" flag: %d\n", cel->flag);
138 wcsprintf(" offset: %d\n", cel->offset);
139 if (undefined(cel->phi0)) {
140 wcsprintf(" phi0: UNDEFINED\n");
141 } else {
142 wcsprintf(" phi0: %9f\n", cel->phi0);
143 }
144 if (undefined(cel->theta0)) {
145 wcsprintf(" theta0: UNDEFINED\n");
146 } else {
147 wcsprintf(" theta0: %9f\n", cel->theta0);
148 }
149 wcsprintf(" ref:");
150 for (i = 0; i < 4; i++) {
151 wcsprintf(" %#- 11.5g", cel->ref[i]);
152 }
153 wcsprintf("\n");
154 wcsprintf(" prj: (see below)\n");
155
156 wcsprintf(" euler:");
157 for (i = 0; i < 5; i++) {
158 wcsprintf(" %#- 11.5g", cel->euler[i]);
159 }
160 wcsprintf("\n");
161 wcsprintf(" latpreq: %d", cel->latpreq);
162 if (cel->latpreq == 0) {
163 wcsprintf(" (not required)\n");
164 } else if (cel->latpreq == 1) {
165 wcsprintf(" (disambiguation)\n");
166 } else if (cel->latpreq == 2) {
167 wcsprintf(" (specification)\n");
168 } else {
169 wcsprintf(" (UNDEFINED)\n");
170 }
171 wcsprintf(" isolat: %d\n", cel->isolat);
172
173 WCSPRINTF_PTR(" err: ", cel->err, "\n");
174 if (cel->err) {
175 wcserr_prt(cel->err, " ");
176 }
177
178 wcsprintf("\n");
179 wcsprintf(" prj.*\n");
180 prjprt(&(cel->prj));
181
182 return 0;
183 }
184
185 //----------------------------------------------------------------------------
186
celperr(const struct celprm * cel,const char * prefix)187 int celperr(const struct celprm *cel, const char *prefix)
188
189 {
190 if (cel == 0x0) return CELERR_NULL_POINTER;
191
192 if (cel->err && wcserr_prt(cel->err, prefix) == 0) {
193 wcserr_prt(cel->prj.err, prefix);
194 }
195
196 return 0;
197 }
198
199
200 //----------------------------------------------------------------------------
201
celset(struct celprm * cel)202 int celset(struct celprm *cel)
203
204 {
205 static const char *function = "celset";
206
207 int status;
208 const double tol = 1.0e-10;
209 double clat0, cphip, cthe0, lat0, lng0, phip, slat0, slz, sphip, sthe0;
210 double latp, latp1, latp2, lngp;
211 double u, v, x, y, z;
212 struct prjprm *celprj;
213 struct wcserr **err;
214
215 if (cel == 0x0) return CELERR_NULL_POINTER;
216 err = &(cel->err);
217
218 // Initialize the projection driver routines.
219 celprj = &(cel->prj);
220 if (cel->offset) {
221 celprj->phi0 = cel->phi0;
222 celprj->theta0 = cel->theta0;
223 } else {
224 // Ensure that these are undefined - no fiducial offset.
225 celprj->phi0 = UNDEFINED;
226 celprj->theta0 = UNDEFINED;
227 }
228
229 if ((status = prjset(celprj))) {
230 return wcserr_set(CEL_ERRMSG(cel_prjerr[status]));
231 }
232
233 // Defaults set by the projection routines.
234 if (undefined(cel->phi0)) {
235 cel->phi0 = celprj->phi0;
236 }
237
238 if (undefined(cel->theta0)) {
239 cel->theta0 = celprj->theta0;
240
241 } else if (fabs(cel->theta0) > 90.0) {
242 if (fabs(cel->theta0) > 90.0 + tol) {
243 return wcserr_set(WCSERR_SET(CELERR_BAD_COORD_TRANS),
244 "Invalid coordinate transformation parameters: theta0 > 90");
245 }
246
247 if (cel->theta0 > 90.0) {
248 cel->theta0 = 90.0;
249 } else {
250 cel->theta0 = -90.0;
251 }
252 }
253
254
255 lng0 = cel->ref[0];
256 lat0 = cel->ref[1];
257 phip = cel->ref[2];
258 latp = cel->ref[3];
259
260 // Set default for native longitude of the celestial pole?
261 if (undefined(phip) || phip == 999.0) {
262 phip = (lat0 < cel->theta0) ? 180.0 : 0.0;
263 phip += cel->phi0;
264
265 if (phip < -180.0) {
266 phip += 360.0;
267 } else if (phip > 180.0) {
268 phip -= 360.0;
269 }
270
271 cel->ref[2] = phip;
272 }
273
274
275 // Compute celestial coordinates of the native pole.
276 cel->latpreq = 0;
277 if (cel->theta0 == 90.0) {
278 // Fiducial point at the native pole.
279 lngp = lng0;
280 latp = lat0;
281
282 } else {
283 // Fiducial point away from the native pole.
284 sincosd(lat0, &slat0, &clat0);
285 sincosd(cel->theta0, &sthe0, &cthe0);
286
287 if (phip == cel->phi0) {
288 sphip = 0.0;
289 cphip = 1.0;
290
291 u = cel->theta0;
292 v = 90.0 - lat0;
293
294 } else {
295 sincosd(phip - cel->phi0, &sphip, &cphip);
296
297 x = cthe0*cphip;
298 y = sthe0;
299 z = sqrt(x*x + y*y);
300 if (z == 0.0) {
301 if (slat0 != 0.0) {
302 return wcserr_set(WCSERR_SET(CELERR_BAD_COORD_TRANS),
303 "Invalid coordinate description:\n"
304 "lat0 == 0 is required for |phip - phi0| = 90 and theta0 == 0");
305 }
306
307 // latp determined solely by LATPOLEa in this case.
308 cel->latpreq = 2;
309 if (latp > 90.0) {
310 latp = 90.0;
311 } else if (latp < -90.0) {
312 latp = -90.0;
313 }
314
315 // Avert a spurious compiler warning.
316 u = v = 0.0;
317
318 } else {
319 slz = slat0/z;
320 if (fabs(slz) > 1.0) {
321 if ((fabs(slz) - 1.0) < tol) {
322 if (slz > 0.0) {
323 slz = 1.0;
324 } else {
325 slz = -1.0;
326 }
327 } else {
328 return wcserr_set(WCSERR_SET(CELERR_BAD_COORD_TRANS),
329 "Invalid coordinate description:\n|lat0| <= %.3f is required "
330 "for these values of phip, phi0, and theta0", asind(z));
331 }
332 }
333
334 u = atan2d(y,x);
335 v = acosd(slz);
336 }
337 }
338
339 if (cel->latpreq == 0) {
340 latp1 = u + v;
341 if (latp1 > 180.0) {
342 latp1 -= 360.0;
343 } else if (latp1 < -180.0) {
344 latp1 += 360.0;
345 }
346
347 latp2 = u - v;
348 if (latp2 > 180.0) {
349 latp2 -= 360.0;
350 } else if (latp2 < -180.0) {
351 latp2 += 360.0;
352 }
353
354 if (fabs(latp1) < 90.0+tol &&
355 fabs(latp2) < 90.0+tol) {
356 // There are two valid solutions for latp.
357 cel->latpreq = 1;
358 }
359
360 if (fabs(latp-latp1) < fabs(latp-latp2)) {
361 if (fabs(latp1) < 90.0+tol) {
362 latp = latp1;
363 } else {
364 latp = latp2;
365 }
366 } else {
367 if (fabs(latp2) < 90.0+tol) {
368 latp = latp2;
369 } else {
370 latp = latp1;
371 }
372 }
373
374 // Account for rounding error.
375 if (fabs(latp) < 90.0+tol) {
376 if (latp > 90.0) {
377 latp = 90.0;
378 } else if (latp < -90.0) {
379 latp = -90.0;
380 }
381 }
382 }
383
384 z = cosd(latp)*clat0;
385 if (fabs(z) < tol) {
386 if (fabs(clat0) < tol) {
387 // Celestial pole at the fiducial point.
388 lngp = lng0;
389
390 } else if (latp > 0.0) {
391 // Celestial north pole at the native pole.
392 lngp = lng0 + phip - cel->phi0 - 180.0;
393
394 } else {
395 // Celestial south pole at the native pole.
396 lngp = lng0 - phip + cel->phi0;
397 }
398
399 } else {
400 x = (sthe0 - sind(latp)*slat0)/z;
401 y = sphip*cthe0/clat0;
402 if (x == 0.0 && y == 0.0) {
403 // Sanity check (shouldn't be possible).
404 return wcserr_set(WCSERR_SET(CELERR_BAD_COORD_TRANS),
405 "Invalid coordinate transformation parameters, internal error");
406 }
407 lngp = lng0 - atan2d(y,x);
408 }
409
410 // Make celestial longitude of the native pole the same sign as at the
411 // fiducial point.
412 if (lng0 >= 0.0) {
413 if (lngp < 0.0) {
414 lngp += 360.0;
415 } else if (lngp > 360.0) {
416 lngp -= 360.0;
417 }
418 } else {
419 if (lngp > 0.0) {
420 lngp -= 360.0;
421 } else if (lngp < -360.0) {
422 lngp += 360.0;
423 }
424 }
425 }
426
427 // Reset LATPOLEa.
428 cel->ref[3] = latp;
429
430 // Set the Euler angles.
431 cel->euler[0] = lngp;
432 cel->euler[1] = 90.0 - latp;
433 cel->euler[2] = phip;
434 sincosd(cel->euler[1], &cel->euler[4], &cel->euler[3]);
435 cel->isolat = (cel->euler[4] == 0.0);
436 cel->flag = CELSET;
437
438 // Check for ill-conditioned parameters.
439 if (fabs(latp) > 90.0+tol) {
440 return wcserr_set(WCSERR_SET(CELERR_ILL_COORD_TRANS),
441 "Ill-conditioned coordinate transformation parameters\nNo valid "
442 "solution for latp for these values of phip, phi0, and theta0");
443 }
444
445 return 0;
446 }
447
448 //----------------------------------------------------------------------------
449
celx2s(struct celprm * cel,int nx,int ny,int sxy,int sll,const double x[],const double y[],double phi[],double theta[],double lng[],double lat[],int stat[])450 int celx2s(
451 struct celprm *cel,
452 int nx,
453 int ny,
454 int sxy,
455 int sll,
456 const double x[],
457 const double y[],
458 double phi[],
459 double theta[],
460 double lng[],
461 double lat[],
462 int stat[])
463
464 {
465 static const char *function = "celx2s";
466
467 int istat, nphi, status = 0;
468 struct prjprm *celprj;
469 struct wcserr **err;
470
471 // Initialize.
472 if (cel == 0x0) return CELERR_NULL_POINTER;
473 err = &(cel->err);
474
475 if (cel->flag != CELSET) {
476 if ((status = celset(cel))) return status;
477 }
478
479 // Apply spherical deprojection.
480 celprj = &(cel->prj);
481 if ((istat = celprj->prjx2s(celprj, nx, ny, sxy, 1, x, y, phi, theta,
482 stat))) {
483 if (istat) {
484 status = wcserr_set(CEL_ERRMSG(cel_prjerr[istat]));
485 if (status != CELERR_BAD_PIX) {
486 return status;
487 }
488 }
489 }
490
491 nphi = (ny > 0) ? (nx*ny) : nx;
492
493 // Compute celestial coordinates.
494 sphx2s(cel->euler, nphi, 0, 1, sll, phi, theta, lng, lat);
495
496 return status;
497 }
498
499 //----------------------------------------------------------------------------
500
cels2x(struct celprm * cel,int nlng,int nlat,int sll,int sxy,const double lng[],const double lat[],double phi[],double theta[],double x[],double y[],int stat[])501 int cels2x(
502 struct celprm *cel,
503 int nlng,
504 int nlat,
505 int sll,
506 int sxy,
507 const double lng[],
508 const double lat[],
509 double phi[],
510 double theta[],
511 double x[],
512 double y[],
513 int stat[])
514
515 {
516 static const char *function = "cels2x";
517
518 int istat, nphi, ntheta, status = 0;
519 struct prjprm *celprj;
520 struct wcserr **err;
521
522 // Initialize.
523 if (cel == 0x0) return CELERR_NULL_POINTER;
524 err = &(cel->err);
525
526 if (cel->flag != CELSET) {
527 if ((status = celset(cel))) return status;
528 }
529
530 // Compute native coordinates.
531 sphs2x(cel->euler, nlng, nlat, sll, 1, lng, lat, phi, theta);
532
533 if (cel->isolat) {
534 // Constant celestial latitude -> constant native latitude.
535 nphi = nlng;
536 ntheta = nlat;
537 } else {
538 nphi = (nlat > 0) ? (nlng*nlat) : nlng;
539 ntheta = 0;
540 }
541
542 // Apply the spherical projection.
543 celprj = &(cel->prj);
544 if ((istat = celprj->prjs2x(celprj, nphi, ntheta, 1, sxy, phi, theta, x, y,
545 stat))) {
546 if (istat) {
547 status = wcserr_set(CEL_ERRMSG(cel_prjerr[istat]));
548 if (status != CELERR_BAD_WORLD) {
549 return status;
550 }
551 }
552 }
553
554 return status;
555 }
556