1 /*============================================================================
2 WCSLIB 7.3 - an implementation of the FITS WCS standard.
3 Copyright (C) 1995-2020, 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 Direct correspondence concerning WCSLIB to mark@calabretta.id.au
21
22 Author: Mark Calabretta, Australia Telescope National Facility, CSIRO.
23 http://www.atnf.csiro.au/people/Mark.Calabretta
24 $Id: wcsfix.c,v 7.3.1.2 2020/08/17 12:10:44 mcalabre Exp mcalabre $
25 *===========================================================================*/
26
27 #include <math.h>
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <string.h>
31
32 #include "wcserr.h"
33 #include "wcsmath.h"
34 #include "wcstrig.h"
35 #include "wcsutil.h"
36 #include "lin.h"
37 #include "sph.h"
38 #include "wcs.h"
39 #include "wcsunits.h"
40 #include "wcsfix.h"
41
42 extern const int WCSSET;
43
44 // Maximum number of coordinate axes that can be handled.
45 #define NMAX 16
46
47 // Map status return value to message.
48 const char *wcsfix_errmsg[] = {
49 "Success",
50 "Null wcsprm pointer passed",
51 "Memory allocation failed",
52 "Linear transformation matrix is singular",
53 "Inconsistent or unrecognized coordinate axis types",
54 "Invalid parameter value",
55 "Invalid coordinate transformation parameters",
56 "Ill-conditioned coordinate transformation parameters",
57 "All of the corner pixel coordinates are invalid",
58 "Could not determine reference pixel coordinate",
59 "Could not determine reference pixel value"};
60
61 // Map error returns for lower-level routines.
62 const int fix_linerr[] = {
63 FIXERR_SUCCESS, // 0: LINERR_SUCCESS
64 FIXERR_NULL_POINTER, // 1: LINERR_NULL_POINTER
65 FIXERR_MEMORY, // 2: LINERR_MEMORY
66 FIXERR_SINGULAR_MTX, // 3: LINERR_SINGULAR_MTX
67 FIXERR_BAD_PARAM, // 4: LINERR_DISTORT_INIT
68 FIXERR_NO_REF_PIX_COORD, // 5: LINERR_DISTORT
69 FIXERR_NO_REF_PIX_VAL // 6: LINERR_DEDISTORT
70 };
71
72 const int fix_wcserr[] = {
73 FIXERR_SUCCESS, // 0: WCSERR_SUCCESS
74 FIXERR_NULL_POINTER, // 1: WCSERR_NULL_POINTER
75 FIXERR_MEMORY, // 2: WCSERR_MEMORY
76 FIXERR_SINGULAR_MTX, // 3: WCSERR_SINGULAR_MTX
77 FIXERR_BAD_CTYPE, // 4: WCSERR_BAD_CTYPE
78 FIXERR_BAD_PARAM, // 5: WCSERR_BAD_PARAM
79 FIXERR_BAD_COORD_TRANS, // 6: WCSERR_BAD_COORD_TRANS
80 FIXERR_ILL_COORD_TRANS, // 7: WCSERR_ILL_COORD_TRANS
81 FIXERR_BAD_CORNER_PIX, // 8: WCSERR_BAD_PIX
82 FIXERR_NO_REF_PIX_VAL, // 9: WCSERR_BAD_WORLD
83 FIXERR_NO_REF_PIX_VAL // 10: WCSERR_BAD_WORLD_COORD
84 // ...others not used
85 };
86
87 // Convenience macro for invoking wcserr_set().
88 #define WCSFIX_ERRMSG(status) WCSERR_SET(status), wcsfix_errmsg[status]
89
90 //----------------------------------------------------------------------------
91
wcsfix(int ctrl,const int naxis[],struct wcsprm * wcs,int stat[])92 int wcsfix(int ctrl, const int naxis[], struct wcsprm *wcs, int stat[])
93
94 {
95 int status = 0;
96
97 if ((stat[CDFIX] = cdfix(wcs)) > 0) {
98 status = 1;
99 }
100
101 if ((stat[DATFIX] = datfix(wcs)) > 0) {
102 status = 1;
103 }
104
105 if ((stat[OBSFIX] = obsfix(0, wcs)) > 0) {
106 status = 1;
107 }
108
109 if ((stat[UNITFIX] = unitfix(ctrl, wcs)) > 0) {
110 status = 1;
111 }
112
113 if ((stat[SPCFIX] = spcfix(wcs)) > 0) {
114 status = 1;
115 }
116
117 if ((stat[CELFIX] = celfix(wcs)) > 0) {
118 status = 1;
119 }
120
121 if ((stat[CYLFIX] = cylfix(naxis, wcs)) > 0) {
122 status = 1;
123 }
124
125 return status;
126 }
127
128 //----------------------------------------------------------------------------
129
wcsfixi(int ctrl,const int naxis[],struct wcsprm * wcs,int stat[],struct wcserr info[])130 int wcsfixi(
131 int ctrl,
132 const int naxis[],
133 struct wcsprm *wcs,
134 int stat[],
135 struct wcserr info[])
136
137 {
138 int ifix, status = 0;
139 struct wcserr err;
140
141 // Handling the status values returned from the sub-fixers is trickier than
142 // it might seem, especially considering that wcs->err may contain an error
143 // status on input which should be preserved if no translation errors occur.
144 // The simplest way seems to be to save a copy of wcs->err and clear it
145 // before each sub-fixer. The last real error to occur, excluding
146 // informative messages, is the one returned.
147
148 // To get informative messages from spcfix() it must precede celfix() and
149 // cylfix(). The latter call wcsset() which also translates AIPS-convention
150 // spectral axes.
151 wcserr_copy(wcs->err, &err);
152
153 for (ifix = CDFIX; ifix < NWCSFIX; ifix++) {
154 // Clear (delete) wcs->err.
155 wcserr_clear(&(wcs->err));
156
157 switch (ifix) {
158 case CDFIX:
159 stat[ifix] = cdfix(wcs);
160 break;
161 case DATFIX:
162 stat[ifix] = datfix(wcs);
163 break;
164 case OBSFIX:
165 stat[ifix] = obsfix(0, wcs);
166 break;
167 case UNITFIX:
168 stat[ifix] = unitfix(ctrl, wcs);
169 break;
170 case SPCFIX:
171 stat[ifix] = spcfix(wcs);
172 break;
173 case CELFIX:
174 stat[ifix] = celfix(wcs);
175 break;
176 case CYLFIX:
177 stat[ifix] = cylfix(naxis, wcs);
178 break;
179 default:
180 continue;
181 }
182
183 if (stat[ifix] == FIXERR_NO_CHANGE) {
184 // No change => no message.
185 wcserr_copy(0x0, info+ifix);
186
187 } else if (stat[ifix] == 0) {
188 // Successful translation, but there may be an informative message.
189 if (wcs->err && wcs->err->status < 0) {
190 wcserr_copy(wcs->err, info+ifix);
191 } else {
192 wcserr_copy(0x0, info+ifix);
193 }
194
195 } else {
196 // An informative message or error message.
197 wcserr_copy(wcs->err, info+ifix);
198
199 if ((status = (stat[ifix] > 0))) {
200 // It was an error, replace the previous one.
201 wcserr_copy(wcs->err, &err);
202 }
203 }
204 }
205
206 // Restore the last error to occur.
207 if (err.status) {
208 wcserr_copy(&err, wcs->err);
209 } else {
210 wcserr_clear(&(wcs->err));
211 }
212
213 return status;
214 }
215
216 //----------------------------------------------------------------------------
217
cdfix(struct wcsprm * wcs)218 int cdfix(struct wcsprm *wcs)
219
220 {
221 int i, k, naxis, status = FIXERR_NO_CHANGE;
222 double *cd;
223
224 if (wcs == 0x0) return FIXERR_NULL_POINTER;
225
226 if ((wcs->altlin & 1) || !(wcs->altlin & 2)) {
227 // Either we have PCi_ja or there are no CDi_ja.
228 return FIXERR_NO_CHANGE;
229 }
230
231 naxis = wcs->naxis;
232 status = FIXERR_NO_CHANGE;
233 for (i = 0; i < naxis; i++) {
234 // Row of zeros?
235 cd = wcs->cd + i * naxis;
236 for (k = 0; k < naxis; k++, cd++) {
237 if (*cd != 0.0) goto next;
238 }
239
240 // Column of zeros?
241 cd = wcs->cd + i;
242 for (k = 0; k < naxis; k++, cd += naxis) {
243 if (*cd != 0.0) goto next;
244 }
245
246 cd = wcs->cd + i * (naxis + 1);
247 *cd = 1.0;
248 status = 0;
249
250 next: ;
251 }
252
253 return status;
254 }
255
256 //----------------------------------------------------------------------------
257
parse_date(const char * buf,int * hour,int * minute,double * sec)258 static int parse_date(const char *buf, int *hour, int *minute, double *sec)
259
260 {
261 char ctmp[72];
262
263 if (sscanf(buf, "%2d:%2d:%s", hour, minute, ctmp) < 3 ||
264 wcsutil_str2double(ctmp, sec)) {
265 return 1;
266 }
267
268 return 0;
269 }
270
271
write_date(char * buf,int hour,int minute,double sec)272 static void write_date(char *buf, int hour, int minute, double sec)
273
274 {
275 char ctmp[32];
276
277 wcsutil_double2str(ctmp, "%04.1f", sec);
278 sprintf(buf, "T%.2d:%.2d:%s", hour, minute, ctmp);
279 }
280
281
newline(char ** cp)282 static char *newline(char **cp)
283
284 {
285 size_t k;
286
287 if ((k = strlen(*cp))) {
288 *cp += k;
289 strcat(*cp, ".\n");
290 *cp += 2;
291 }
292
293 return *cp;
294 }
295
296
datfix(struct wcsprm * wcs)297 int datfix(struct wcsprm *wcs)
298
299 {
300 static const char *function = "datfix";
301
302 // MJD of J2000.0 and B1900.0.
303 const double mjd2000 = 51544.5;
304 const double mjd1900 = 15019.81352;
305
306 // Days per Julian year and per tropical year.
307 const double djy = 365.25;
308 const double dty = 365.242198781;
309
310 const char *dateid;
311 char *cp, *date, infomsg[512], orig_date[72];
312 int day, dd, hour = 0, i, jd, minute = 0, month, msec, n4, status, year;
313 double bepoch, jepoch, mjd[2], mjdsum, mjdtmp, sec = 0.0, t, *wcsmjd;
314 struct wcserr **err;
315
316 if (wcs == 0x0) return FIXERR_NULL_POINTER;
317 err = &(wcs->err);
318
319 cp = infomsg;
320 *cp = '\0';
321 status = FIXERR_NO_CHANGE;
322
323 for (i = 0; i < 5; i++) {
324 // MJDREF is split into integer and fractional parts, wheres MJDOBS and
325 // the rest are a single value.
326 if (i == 0) {
327 // Note, DATEREF and MJDREF, not DATE-REF and MJD-REF (sigh).
328 dateid = "REF";
329 date = wcs->dateref;
330 wcsmjd = wcs->mjdref;
331 } else if (i == 1) {
332 dateid = "-OBS";
333 date = wcs->dateobs;
334 wcsmjd = &(wcs->mjdobs);
335 } else if (i == 2) {
336 dateid = "-BEG";
337 date = wcs->datebeg;
338 wcsmjd = &(wcs->mjdbeg);
339 } else if (i == 3) {
340 dateid = "-AVG";
341 date = wcs->dateavg;
342 wcsmjd = &(wcs->mjdavg);
343 } else if (i == 4) {
344 dateid = "-END";
345 date = wcs->dateend;
346 wcsmjd = &(wcs->mjdend);
347 }
348
349 strncpy(orig_date, date, 72);
350
351 if (date[0] == '\0') {
352 // Fill in DATE from MJD if possible.
353
354 if (i == 1 && undefined(*wcsmjd)) {
355 // See if we have jepoch or bepoch.
356 if (!undefined(wcs->jepoch)) {
357 *wcsmjd = mjd2000 + (wcs->jepoch - 2000.0)*djy;
358 sprintf(newline(&cp), "Set MJD-OBS to %.6f from JEPOCH", *wcsmjd);
359
360 } else if (!undefined(wcs->bepoch)) {
361 *wcsmjd = mjd1900 + (wcs->bepoch - 1900.0)*dty;
362 sprintf(newline(&cp), "Set MJD-OBS to %.6f from BEPOCH", *wcsmjd);
363 }
364 }
365
366 if (undefined(*wcsmjd)) {
367 // No date information was provided.
368
369 } else {
370 // Calendar date from MJD, with allowance for MJD < 0.
371 if (i == 0) {
372 // MJDREF is already split into integer and fractional parts.
373 mjd[0] = wcsmjd[0];
374 mjd[1] = wcsmjd[1];
375 if (1.0 < mjd[1]) {
376 // Ensure the fractional part lies between 0 and +1.
377 t = floor(mjd[1]);
378 mjd[0] += t;
379 mjd[1] -= t;
380 }
381 } else {
382 // Split it into integer and fractional parts.
383 mjd[0] = floor(*wcsmjd);
384 mjd[1] = *wcsmjd - mjd[0];
385 }
386
387 jd = 2400001 + (int)mjd[0];
388
389 n4 = 4*(jd + ((2*((4*jd - 17918)/146097)*3)/4 + 1)/2 - 37);
390 dd = 10*(((n4-237)%1461)/4) + 5;
391
392 year = n4/1461 - 4712;
393 month = (2 + dd/306)%12 + 1;
394 day = (dd%306)/10 + 1;
395 sprintf(date, "%.4d-%.2d-%.2d", year, month, day);
396
397 // Write time part only if non-zero.
398 if (0.0 < (t = mjd[1])) {
399 t *= 24.0;
400 hour = (int)t;
401 t = 60.0 * (t - hour);
402 minute = (int)t;
403 sec = 60.0 * (t - minute);
404
405 // Round to 1ms.
406 dd = 60000*(60*hour + minute) + (int)(1000*(sec+0.0005));
407 hour = dd / 3600000;
408 dd -= 3600000 * hour;
409 minute = dd / 60000;
410 msec = dd - 60000 * minute;
411 sprintf(date+10, "T%.2d:%.2d:%.2d", hour, minute, msec/1000);
412
413 // Write fractions of a second only if non-zero.
414 if (msec%1000) {
415 sprintf(date+19, ".%.3d", msec%1000);
416 }
417 }
418 }
419
420 } else {
421 if (strlen(date) < 8) {
422 // Can't be a valid date.
423 status = FIXERR_BAD_PARAM;
424 sprintf(newline(&cp), "Invalid DATE%s format '%s' is too short",
425 dateid, date);
426 continue;
427 }
428
429 // Identify the date format.
430 if (date[4] == '-' && date[7] == '-') {
431 // Standard year-2000 form: CCYY-MM-DD[Thh:mm:ss[.sss...]]
432 if (sscanf(date, "%4d-%2d-%2d", &year, &month, &day) < 3) {
433 status = FIXERR_BAD_PARAM;
434 sprintf(newline(&cp), "Invalid DATE%s format '%s'", dateid, date);
435 continue;
436 }
437
438 if (date[10] == 'T') {
439 if (parse_date(date+11, &hour, &minute, &sec)) {
440 status = FIXERR_BAD_PARAM;
441 sprintf(newline(&cp), "Invalid time in DATE%s '%s'", dateid,
442 date+11);
443 continue;
444 }
445 } else if (date[10] == ' ') {
446 hour = 0;
447 minute = 0;
448 sec = 0.0;
449 if (parse_date(date+11, &hour, &minute, &sec)) {
450 write_date(date+10, hour, minute, sec);
451 } else {
452 date[10] = 'T';
453 }
454 }
455
456 } else if (date[4] == '/' && date[7] == '/') {
457 // Also allow CCYY/MM/DD[Thh:mm:ss[.sss...]]
458 if (sscanf(date, "%4d/%2d/%2d", &year, &month, &day) < 3) {
459 status = FIXERR_BAD_PARAM;
460 sprintf(newline(&cp), "Invalid DATE%s format '%s'", dateid, date);
461 continue;
462 }
463
464 if (date[10] == 'T') {
465 if (parse_date(date+11, &hour, &minute, &sec)) {
466 status = FIXERR_BAD_PARAM;
467 sprintf(newline(&cp), "Invalid time in DATE%s '%s'", dateid,
468 date+11);
469 continue;
470 }
471 } else if (date[10] == ' ') {
472 hour = 0;
473 minute = 0;
474 sec = 0.0;
475 if (parse_date(date+11, &hour, &minute, &sec)) {
476 write_date(date+10, hour, minute, sec);
477 } else {
478 date[10] = 'T';
479 }
480 }
481
482 // Looks ok, fix it up.
483 date[4] = '-';
484 date[7] = '-';
485
486 } else {
487 if (i == 1 && date[2] == '/' && date[5] == '/') {
488 // Old format DATE-OBS date: DD/MM/YY, also allowing DD/MM/CCYY.
489 if (sscanf(date, "%2d/%2d/%4d", &day, &month, &year) < 3) {
490 status = FIXERR_BAD_PARAM;
491 sprintf(newline(&cp), "Invalid DATE%s format '%s'", dateid,
492 date);
493 continue;
494 }
495
496 } else if (i == 1 && date[2] == '-' && date[5] == '-') {
497 // Also recognize DD-MM-YY and DD-MM-CCYY
498 if (sscanf(date, "%2d-%2d-%4d", &day, &month, &year) < 3) {
499 status = FIXERR_BAD_PARAM;
500 sprintf(newline(&cp), "Invalid DATE%s format '%s'", dateid,
501 date);
502 continue;
503 }
504
505 } else {
506 // Not a valid date format.
507 status = FIXERR_BAD_PARAM;
508 sprintf(newline(&cp), "Invalid DATE%s format '%s'", dateid, date);
509 continue;
510 }
511
512 if (year < 100) year += 1900;
513
514 // Doesn't have a time.
515 sprintf(date, "%.4d-%.2d-%.2d", year, month, day);
516 }
517
518 // Compute MJD.
519 mjd[0] = (double)((1461*(year - (12-month)/10 + 4712))/4
520 + (306*((month+9)%12) + 5)/10
521 - (3*((year - (12-month)/10 + 4900)/100))/4
522 + day - 2399904);
523 mjd[1] = (hour + (minute + sec/60.0)/60.0)/24.0;
524 mjdsum = mjd[0] + mjd[1];
525
526 if (undefined(*wcsmjd)) {
527 if (i == 0) {
528 wcsmjd[0] = mjd[0];
529 wcsmjd[1] = mjd[1];
530 } else {
531 *wcsmjd = mjdsum;
532 }
533 sprintf(newline(&cp), "Set MJD%s to %.6f from DATE%s", dateid,
534 mjdsum, dateid);
535
536 } else {
537 // Check for consistency.
538 if (i == 0) {
539 mjdtmp = wcsmjd[0] + wcsmjd[1];
540 } else {
541 mjdtmp = *wcsmjd;
542 }
543
544 if (0.001 < fabs(mjdsum - mjdtmp)) {
545 status = FIXERR_BAD_PARAM;
546 sprintf(newline(&cp),
547 "Invalid parameter values: MJD%s and DATE%s are inconsistent",
548 dateid, dateid);
549 }
550 }
551
552 if (i == 1) {
553 if (!undefined(wcs->jepoch)) {
554 // Check consistency of JEPOCH.
555 jepoch = 2000.0 + (*wcsmjd - mjd2000) / djy;
556
557 if (0.000002 < fabs(jepoch - wcs->jepoch)) {
558 // Informational only, no error.
559 sprintf(newline(&cp), "JEPOCH is inconsistent with DATE-OBS");
560 }
561 }
562
563 if (!undefined(wcs->bepoch)) {
564 // Check consistency of BEPOCH.
565 bepoch = 1900.0 + (*wcsmjd - mjd1900) / dty;
566
567 if (0.000002 < fabs(bepoch - wcs->bepoch)) {
568 // Informational only, no error.
569 sprintf(newline(&cp), "BEPOCH is inconsistent with DATE-OBS");
570 }
571 }
572 }
573 }
574
575 if (strncmp(orig_date, date, 72)) {
576 if (orig_date[0] == '\0') {
577 sprintf(newline(&cp), "Set DATE%s to '%s' from MJD%s", dateid, date,
578 dateid);
579 } else {
580 sprintf(newline(&cp), "Changed DATE%s from '%s' to '%s'", dateid,
581 orig_date, date);
582 }
583
584 if (status == FIXERR_NO_CHANGE) status = 0;
585 }
586 }
587
588 if (*infomsg) {
589 wcserr_set(WCSERR_SET(FIXERR_DATE_FIX), infomsg);
590 }
591
592 return status;
593 }
594
595 //----------------------------------------------------------------------------
596
obsfix(int ctrl,struct wcsprm * wcs)597 int obsfix(int ctrl, struct wcsprm *wcs)
598
599 {
600 static const char *function = "obsfix";
601
602 // IAU(1976) ellipsoid (as prescribed by WCS Paper VII).
603 const double a = 6378140.0, f = 1.0 / 298.2577;
604 const double e2 = (2.0 - f)*f;
605
606 char *cp, infomsg[256];
607 int havelbh = 7, havexyz = 7, i, status;
608 size_t k;
609 double coslat, coslng, d, hgt, lat, lng, n, r2, rho, sinlat, sinlng, x, y,
610 z, zeta;
611 struct wcserr **err;
612
613 if (wcs == 0x0) return FIXERR_NULL_POINTER;
614 err = &(wcs->err);
615
616 // Set masks for checking partially-defined coordinate triplets.
617 havexyz -= 1*undefined(wcs->obsgeo[0]);
618 havexyz -= 2*undefined(wcs->obsgeo[1]);
619 havexyz -= 4*undefined(wcs->obsgeo[2]);
620 havelbh -= 1*undefined(wcs->obsgeo[3]);
621 havelbh -= 2*undefined(wcs->obsgeo[4]);
622 havelbh -= 4*undefined(wcs->obsgeo[5]);
623
624 if (ctrl == 2) {
625 // Make no changes.
626 if (0 < havexyz && havexyz < 7) {
627 return wcserr_set(WCSERR_SET(FIXERR_BAD_PARAM),
628 "Partially undefined Cartesian coordinate triplet");
629 }
630
631 if (0 < havelbh && havelbh < 7) {
632 return wcserr_set(WCSERR_SET(FIXERR_BAD_PARAM),
633 "Partially undefined Geodetic coordinate triplet");
634 }
635
636 if (havexyz == 0 || havelbh == 0) {
637 return FIXERR_NO_CHANGE;
638 }
639 }
640
641 if (havexyz == 0 && havelbh == 0) {
642 return FIXERR_NO_CHANGE;
643 }
644
645
646 infomsg[0] = '\0';
647 status = FIXERR_NO_CHANGE;
648
649 if (havelbh == 7) {
650 // Compute (x,y,z) from (lng,lat,hgt).
651 sincosd(wcs->obsgeo[3], &sinlng, &coslng);
652 sincosd(wcs->obsgeo[4], &sinlat, &coslat);
653 n = a / sqrt(1.0 - e2*sinlat*sinlat);
654 rho = n + wcs->obsgeo[5];
655
656 x = rho*coslng*coslat;
657 y = rho*sinlng*coslat;
658 z = (rho - n*e2)*sinlat;
659
660 if (havexyz < 7) {
661 // One or more of the Cartesian elements was undefined.
662 status = 0;
663 cp = infomsg;
664
665 if (ctrl == 1 || !(havexyz & 1)) {
666 wcs->obsgeo[0] = x;
667 sprintf(cp, "%s OBSGEO-X to %12.3f from OBSGEO-[LBH]",
668 (havexyz & 1) ? "Reset" : "Set", x);
669 }
670
671 if (ctrl == 1 || !(havexyz & 2)) {
672 wcs->obsgeo[1] = y;
673
674 if ((k = strlen(cp))) {
675 strcat(cp+k, ".\n");
676 cp += k + 2;
677 }
678
679 sprintf(cp, "%s OBSGEO-Y to %12.3f from OBSGEO-[LBH]",
680 (havexyz & 2) ? "Reset" : "Set", y);
681 }
682
683 if (ctrl == 1 || !(havexyz & 4)) {
684 wcs->obsgeo[2] = z;
685 if ((k = strlen(cp))) {
686 strcat(cp+k, ".\n");
687 cp += k + 2;
688 }
689
690 sprintf(cp, "%s OBSGEO-Z to %12.3f from OBSGEO-[LBH]",
691 (havexyz & 4) ? "Reset" : "Set", z);
692 }
693
694 wcserr_set(WCSERR_SET(FIXERR_OBSGEO_FIX), infomsg);
695
696 if (havexyz == 0) {
697 // Skip the consistency check.
698 return status;
699 }
700 }
701
702 } else if (havexyz == 7) {
703 // Compute (lng,lat,hgt) from (x,y,z).
704 x = wcs->obsgeo[0];
705 y = wcs->obsgeo[1];
706 z = wcs->obsgeo[2];
707 r2 = x*x + y*y;
708
709 // Iterate over the value of zeta.
710 zeta = z;
711 for (i = 0; i < 4; i++) {
712 rho = sqrt(r2 + zeta*zeta);
713 sinlat = zeta / rho;
714 n = a / sqrt(1.0 - e2*sinlat*sinlat);
715
716 zeta = z / (1.0 - n*e2/rho);
717 }
718
719 lng = atan2d(y, x);
720 lat = asind(sinlat);
721 hgt = rho - n;
722
723 if (havelbh < 7) {
724 // One or more of the Geodetic elements was undefined.
725 status = 0;
726 cp = infomsg;
727
728 if (ctrl == 1 || !(havelbh & 1)) {
729 wcs->obsgeo[3] = lng;
730 sprintf(cp, "%s OBSGEO-L to %12.6f from OBSGEO-[XYZ]",
731 (havelbh & 1) ? "Reset" : "Set", lng);
732 }
733
734 if (ctrl == 1 || !(havelbh & 2)) {
735 wcs->obsgeo[4] = lat;
736 if ((k = strlen(cp))) {
737 strcat(cp+k, ".\n");
738 cp += k + 2;
739 }
740
741 sprintf(cp, "%s OBSGEO-B to %12.6f from OBSGEO-[XYZ]",
742 (havelbh & 2) ? "Reset" : "Set", lat);
743 }
744
745 if (ctrl == 1 || !(havelbh & 4)) {
746 wcs->obsgeo[5] = hgt;
747 if ((k = strlen(cp))) {
748 strcat(cp+k, ".\n");
749 cp += k + 2;
750 }
751
752 sprintf(cp, "%s OBSGEO-H to %12.3f from OBSGEO-[XYZ]",
753 (havelbh & 4) ? "Reset" : "Set", hgt);
754 }
755
756 wcserr_set(WCSERR_SET(FIXERR_OBSGEO_FIX), infomsg);
757
758 if (havelbh == 0) {
759 // Skip the consistency check.
760 return status;
761 }
762 }
763
764 // Compute (x,y,z) from (lng,lat,hgt) for consistency checking.
765 sincosd(wcs->obsgeo[3], &sinlng, &coslng);
766 sincosd(wcs->obsgeo[4], &sinlat, &coslat);
767 n = a / sqrt(1.0 - e2*sinlat*sinlat);
768 rho = n + wcs->obsgeo[5];
769
770 x = rho*coslng*coslat;
771 y = rho*sinlng*coslat;
772 z = (rho - n*e2)*sinlat;
773
774 } else {
775 return wcserr_set(WCSERR_SET(FIXERR_BAD_PARAM),
776 "Observatory coordinates incomplete");
777 }
778
779
780 // Check consistency.
781 r2 = 0.0;
782 d = wcs->obsgeo[0] - x;
783 r2 += d*d;
784 d = wcs->obsgeo[1] - y;
785 r2 += d*d;
786 d = wcs->obsgeo[2] - z;
787 r2 += d*d;
788
789 if (1.0 < r2) {
790 d = sqrt(r2);
791 return wcserr_set(WCSERR_SET(FIXERR_BAD_PARAM),
792 "Observatory coordinates inconsistent by %.1f metres", d);
793 }
794
795 return status;
796 }
797
798
799 //----------------------------------------------------------------------------
800
unitfix(int ctrl,struct wcsprm * wcs)801 int unitfix(int ctrl, struct wcsprm *wcs)
802
803 {
804 const char *function = "unitfix";
805
806 char orig_unit[72], msg[512], msgtmp[192];
807 int i, result, status = FIXERR_NO_CHANGE;
808 size_t msglen;
809 struct wcserr **err;
810
811 if (wcs == 0x0) return FIXERR_NULL_POINTER;
812 err = &(wcs->err);
813
814 strncpy(msg, "Changed units:", 512);
815
816 for (i = 0; i < wcs->naxis; i++) {
817 strncpy(orig_unit, wcs->cunit[i], 71);
818 result = wcsutrne(ctrl, wcs->cunit[i], &(wcs->err));
819 if (result == 0 || result == 12) {
820 msglen = strlen(msg);
821 if (msglen < 511) {
822 wcsutil_null_fill(72, orig_unit);
823 sprintf(msgtmp, "\n '%s' -> '%s',", orig_unit, wcs->cunit[i]);
824 strncpy(msg+msglen, msgtmp, 511-msglen);
825 status = FIXERR_UNITS_ALIAS;
826 }
827 }
828 }
829
830 if (status == FIXERR_UNITS_ALIAS) {
831 // Chop off the trailing ", ".
832 msglen = strlen(msg) - 2;
833 msg[msglen] = '\0';
834 wcserr_set(WCSERR_SET(FIXERR_UNITS_ALIAS), msg);
835
836 status = 0;
837 }
838
839 return status;
840 }
841
842 //----------------------------------------------------------------------------
843
spcfix(struct wcsprm * wcs)844 int spcfix(struct wcsprm *wcs)
845
846 {
847 static const char *function = "spcfix";
848
849 char ctype[9], specsys[9];
850 int i, status;
851 struct wcserr **err;
852
853 if (wcs == 0x0) return FIXERR_NULL_POINTER;
854 err = &(wcs->err);
855
856 for (i = 0; i < wcs->naxis; i++) {
857 // Translate an AIPS-convention spectral type if present.
858 status = spcaips(wcs->ctype[i], wcs->velref, ctype, specsys);
859 if (status == 0) {
860 // An AIPS type was found but it may match what we already have.
861 status = FIXERR_NO_CHANGE;
862
863 // Was specsys translated?
864 if (wcs->specsys[0] == '\0' && *specsys) {
865 strncpy(wcs->specsys, specsys, 9);
866 wcserr_set(WCSERR_SET(FIXERR_SPC_UPDATE),
867 "Changed SPECSYS to '%s'", specsys);
868 status = 0;
869 }
870
871 // Was ctype translated? Have to null-fill for comparing them.
872 wcsutil_null_fill(9, wcs->ctype[i]);
873 if (strncmp(wcs->ctype[i], ctype, 9)) {
874 // ctype was translated...
875 if (status == 0) {
876 // ...and specsys was also.
877 wcserr_set(WCSERR_SET(FIXERR_SPC_UPDATE),
878 "Changed CTYPE%d from '%s' to '%s', and SPECSYS to '%s' "
879 "(VELREF=%d)", i+1, wcs->ctype[i], ctype, wcs->specsys,
880 wcs->velref);
881 } else {
882 wcserr_set(WCSERR_SET(FIXERR_SPC_UPDATE),
883 "Changed CTYPE%d from '%s' to '%s' (VELREF=%d)", i+1,
884 wcs->ctype[i], ctype, wcs->velref);
885 status = 0;
886 }
887
888 strncpy(wcs->ctype[i], ctype, 9);
889 }
890
891 // Tidy up.
892 if (status == 0) {
893 wcsutil_null_fill(72, wcs->ctype[i]);
894 wcsutil_null_fill(72, wcs->specsys);
895 }
896
897 // No need to check for others, wcsset() will fail if so.
898 return status;
899
900 } else if (status == SPCERR_BAD_SPEC_PARAMS) {
901 // An AIPS spectral type was found but with invalid velref.
902 return wcserr_set(WCSERR_SET(FIXERR_BAD_PARAM),
903 "Invalid parameter value: velref = %d", wcs->velref);
904 }
905 }
906
907 return FIXERR_NO_CHANGE;
908 }
909
910 //----------------------------------------------------------------------------
911
celfix(struct wcsprm * wcs)912 int celfix(struct wcsprm *wcs)
913
914 {
915 static const char *function = "celfix";
916
917 int k, status;
918 struct celprm *wcscel = &(wcs->cel);
919 struct prjprm *wcsprj = &(wcscel->prj);
920 struct wcserr **err;
921
922 if (wcs == 0x0) return FIXERR_NULL_POINTER;
923 err = &(wcs->err);
924
925 // Initialize if required.
926 if (wcs->flag != WCSSET) {
927 if ((status = wcsset(wcs))) return fix_wcserr[status];
928 }
929
930 // Was an NCP or GLS projection code translated?
931 if (wcs->lat >= 0) {
932 // Check ctype.
933 if (strcmp(wcs->ctype[wcs->lat]+5, "NCP") == 0) {
934 strcpy(wcs->ctype[wcs->lng]+5, "SIN");
935 strcpy(wcs->ctype[wcs->lat]+5, "SIN");
936
937 if (wcs->npvmax < wcs->npv + 2) {
938 // Allocate space for two more PVi_ma keyvalues.
939 if (wcs->m_flag == WCSSET && wcs->pv == wcs->m_pv) {
940 if (!(wcs->pv = calloc(wcs->npv+2, sizeof(struct pvcard)))) {
941 wcs->pv = wcs->m_pv;
942 return wcserr_set(WCSFIX_ERRMSG(FIXERR_MEMORY));
943 }
944
945 wcs->npvmax = wcs->npv + 2;
946 wcs->m_flag = WCSSET;
947
948 for (k = 0; k < wcs->npv; k++) {
949 wcs->pv[k] = wcs->m_pv[k];
950 }
951
952 if (wcs->m_pv) free(wcs->m_pv);
953 wcs->m_pv = wcs->pv;
954
955 } else {
956 return wcserr_set(WCSFIX_ERRMSG(FIXERR_MEMORY));
957 }
958 }
959
960 wcs->pv[wcs->npv].i = wcs->lat + 1;
961 wcs->pv[wcs->npv].m = 1;
962 wcs->pv[wcs->npv].value = wcsprj->pv[1];
963 (wcs->npv)++;
964
965 wcs->pv[wcs->npv].i = wcs->lat + 1;
966 wcs->pv[wcs->npv].m = 2;
967 wcs->pv[wcs->npv].value = wcsprj->pv[2];
968 (wcs->npv)++;
969
970 return 0;
971
972 } else if (strcmp(wcs->ctype[wcs->lat]+5, "GLS") == 0) {
973 strcpy(wcs->ctype[wcs->lng]+5, "SFL");
974 strcpy(wcs->ctype[wcs->lat]+5, "SFL");
975
976 if (wcs->crval[wcs->lng] != 0.0 || wcs->crval[wcs->lat] != 0.0) {
977 // In the AIPS convention, setting the reference longitude and
978 // latitude for GLS does not create an oblique graticule. A non-zero
979 // reference longitude introduces an offset in longitude in the normal
980 // way, whereas a non-zero reference latitude simply translates the
981 // reference point (i.e. the map as a whole) to that latitude. This
982 // might be effected by adjusting CRPIXja but that is complicated by
983 // the linear transformation and instead is accomplished here by
984 // setting theta_0.
985 if (wcs->npvmax < wcs->npv + 3) {
986 // Allocate space for three more PVi_ma keyvalues.
987 if (wcs->m_flag == WCSSET && wcs->pv == wcs->m_pv) {
988 if (!(wcs->pv = calloc(wcs->npv+3, sizeof(struct pvcard)))) {
989 wcs->pv = wcs->m_pv;
990 return wcserr_set(WCSFIX_ERRMSG(FIXERR_MEMORY));
991 }
992
993 wcs->npvmax = wcs->npv + 3;
994 wcs->m_flag = WCSSET;
995
996 for (k = 0; k < wcs->npv; k++) {
997 wcs->pv[k] = wcs->m_pv[k];
998 }
999
1000 if (wcs->m_pv) free(wcs->m_pv);
1001 wcs->m_pv = wcs->pv;
1002
1003 } else {
1004 return wcserr_set(WCSFIX_ERRMSG(FIXERR_MEMORY));
1005 }
1006 }
1007
1008 wcs->pv[wcs->npv].i = wcs->lng + 1;
1009 wcs->pv[wcs->npv].m = 0;
1010 wcs->pv[wcs->npv].value = 1.0;
1011 (wcs->npv)++;
1012
1013 // Note that the reference longitude is still zero.
1014 wcs->pv[wcs->npv].i = wcs->lng + 1;
1015 wcs->pv[wcs->npv].m = 1;
1016 wcs->pv[wcs->npv].value = 0.0;
1017 (wcs->npv)++;
1018
1019 wcs->pv[wcs->npv].i = wcs->lng + 1;
1020 wcs->pv[wcs->npv].m = 2;
1021 wcs->pv[wcs->npv].value = wcs->crval[wcs->lat];
1022 (wcs->npv)++;
1023 }
1024
1025 return 0;
1026 }
1027 }
1028
1029 return FIXERR_NO_CHANGE;
1030 }
1031
1032 //----------------------------------------------------------------------------
1033
cylfix(const int naxis[],struct wcsprm * wcs)1034 int cylfix(const int naxis[], struct wcsprm *wcs)
1035
1036 {
1037 static const char *function = "cylfix";
1038
1039 unsigned short icnr, indx[NMAX], ncnr;
1040 int j, k, stat[4], status;
1041 double img[4][NMAX], lat, lng, phi[4], phi0, phimax, phimin, pix[4][NMAX],
1042 *pixj, theta[4], theta0, world[4][NMAX], x, y;
1043 struct wcserr **err;
1044
1045 if (naxis == 0x0) return FIXERR_NO_CHANGE;
1046 if (wcs == 0x0) return FIXERR_NULL_POINTER;
1047 err = &(wcs->err);
1048
1049 // Initialize if required.
1050 if (wcs->flag != WCSSET) {
1051 if ((status = wcsset(wcs))) return fix_wcserr[status];
1052 }
1053
1054 // Check that we have a cylindrical projection.
1055 if (wcs->cel.prj.category != CYLINDRICAL) return FIXERR_NO_CHANGE;
1056 if (wcs->naxis < 2) return FIXERR_NO_CHANGE;
1057
1058
1059 // Compute the native longitude in each corner of the image.
1060 ncnr = 1 << wcs->naxis;
1061
1062 for (k = 0; k < NMAX; k++) {
1063 indx[k] = 1 << k;
1064 }
1065
1066 phimin = 1.0e99;
1067 phimax = -1.0e99;
1068 for (icnr = 0; icnr < ncnr;) {
1069 // Do four corners at a time.
1070 for (j = 0; j < 4; j++, icnr++) {
1071 pixj = pix[j];
1072
1073 for (k = 0; k < wcs->naxis; k++) {
1074 if (icnr & indx[k]) {
1075 *(pixj++) = naxis[k] + 0.5;
1076 } else {
1077 *(pixj++) = 0.5;
1078 }
1079 }
1080 }
1081
1082 if (!(status = wcsp2s(wcs, 4, NMAX, pix[0], img[0], phi, theta, world[0],
1083 stat))) {
1084 for (j = 0; j < 4; j++) {
1085 if (phi[j] < phimin) phimin = phi[j];
1086 if (phi[j] > phimax) phimax = phi[j];
1087 }
1088 }
1089 }
1090
1091 if (phimin > phimax) return fix_wcserr[status];
1092
1093 // Any changes needed?
1094 if (phimin >= -180.0 && phimax <= 180.0) return FIXERR_NO_CHANGE;
1095
1096
1097 // Compute the new reference pixel coordinates.
1098 phi0 = (phimin + phimax) / 2.0;
1099 theta0 = 0.0;
1100
1101 if ((status = prjs2x(&(wcs->cel.prj), 1, 1, 1, 1, &phi0, &theta0, &x, &y,
1102 stat))) {
1103 if (status == PRJERR_BAD_PARAM) {
1104 status = FIXERR_BAD_PARAM;
1105 } else {
1106 status = FIXERR_NO_REF_PIX_COORD;
1107 }
1108 return wcserr_set(WCSFIX_ERRMSG(status));
1109 }
1110
1111 for (k = 0; k < wcs->naxis; k++) {
1112 img[0][k] = 0.0;
1113 }
1114 img[0][wcs->lng] = x;
1115 img[0][wcs->lat] = y;
1116
1117 if ((status = linx2p(&(wcs->lin), 1, 0, img[0], pix[0]))) {
1118 return wcserr_set(WCSFIX_ERRMSG(fix_linerr[status]));
1119 }
1120
1121
1122 // Compute celestial coordinates at the new reference pixel.
1123 if ((status = wcsp2s(wcs, 1, 0, pix[0], img[0], phi, theta, world[0],
1124 stat))) {
1125 return fix_wcserr[status];
1126 }
1127
1128 // Compute native coordinates of the celestial pole.
1129 lng = 0.0;
1130 lat = 90.0;
1131 (void)sphs2x(wcs->cel.euler, 1, 1, 1, 1, &lng, &lat, phi, theta);
1132
1133 wcs->crpix[wcs->lng] = pix[0][wcs->lng];
1134 wcs->crpix[wcs->lat] = pix[0][wcs->lat];
1135 wcs->crval[wcs->lng] = world[0][wcs->lng];
1136 wcs->crval[wcs->lat] = world[0][wcs->lat];
1137 wcs->lonpole = phi[0] - phi0;
1138
1139 return wcsset(wcs);
1140 }
1141