1 /*** File libwcs/wcsinit.c
2 *** July 24, 2016
3 *** By Jessica Mink, jmink@cfa.harvard.edu
4 *** Harvard-Smithsonian Center for Astrophysics
5 *** Copyright (C) 1998-2016
6 *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA
7
8 This library is free software; you can redistribute it and/or
9 modify it under the terms of the GNU Lesser General Public
10 License as published by the Free Software Foundation; either
11 version 2 of the License, or (at your option) any later version.
12
13 This library is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 Lesser General Public License for more details.
17
18 You should have received a copy of the GNU Lesser General Public
19 License along with this library; if not, write to the Free Software
20 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21
22 Correspondence concerning WCSTools should be addressed as follows:
23 Internet email: jmink@cfa.harvard.edu
24 Postal address: Jessica Mink
25 Smithsonian Astrophysical Observatory
26 60 Garden St.
27 Cambridge, MA 02138 USA
28
29 * Module: wcsinit.c (World Coordinate Systems)
30 * Purpose: Convert FITS WCS to pixels and vice versa:
31 * Subroutine: wcsinit (hstring) sets a WCS structure from an image header
32 * Subroutine: wcsninit (hstring,lh) sets a WCS structure from fixed-length header
33 * Subroutine: wcsinitn (hstring, name) sets a WCS structure for specified WCS
34 * Subroutine: wcsninitn (hstring,lh, name) sets a WCS structure for specified WCS
35 * Subroutine: wcsinitc (hstring, mchar) sets a WCS structure if multiple
36 * Subroutine: wcsninitc (hstring,lh,mchar) sets a WCS structure if multiple
37 * Subroutine: wcschar (hstring, name) returns suffix for specifed WCS
38 * Subroutine: wcseq (hstring, wcs) set radecsys and equinox from image header
39 * Subroutine: wcseqm (hstring, wcs, mchar) set radecsys and equinox if multiple
40 */
41
42 #include <string.h> /* strstr, NULL */
43 #include <stdio.h> /* stderr */
44 #include <math.h>
45 #include "wcs.h"
46 #ifndef VMS
47 #include <stdlib.h>
48 #endif
49
50 static void wcseq();
51 static void wcseqm();
52 static void wcsioset();
53 void wcsrotset();
54 char wcschar();
55
56 /* set up a WCS structure from a FITS image header lhstring bytes long
57 * for a specified WCS name */
58
59 struct WorldCoor *
wcsninitn(hstring,lhstring,name)60 wcsninitn (hstring, lhstring, name)
61
62 const char *hstring; /* character string containing FITS header information
63 in the format <keyword>= <value> [/ <comment>] */
64 int lhstring; /* Length of FITS header in bytes */
65 const char *name; /* character string with identifying name of WCS */
66 {
67 hlength (hstring, lhstring);
68 return (wcsinitn (hstring, name));
69 }
70
71
72 /* set up a WCS structure from a FITS image header for specified WCSNAME */
73
74 struct WorldCoor *
wcsinitn(hstring,name)75 wcsinitn (hstring, name)
76
77 const char *hstring; /* character string containing FITS header information
78 in the format <keyword>= <value> [/ <comment>] */
79 const char *name; /* character string with identifying name of WCS */
80 {
81 char mchar; /* Suffix character for one of multiple WCS */
82
83 mchar = wcschar (hstring, name);
84 if (mchar == '_') {
85 fprintf (stderr, "WCSINITN: WCS name %s not matched in FITS header\n",
86 name);
87 return (NULL);
88 }
89 return (wcsinitc (hstring, &mchar));
90 }
91
92
93 /* WCSCHAR -- Find the letter for a specific WCS conversion */
94
95 char
wcschar(hstring,name)96 wcschar (hstring, name)
97
98 const char *hstring; /* character string containing FITS header information
99 in the format <keyword>= <value> [/ <comment>] */
100 const char *name; /* Name of WCS conversion to be matched
101 (case-independent) */
102 {
103 char *upname;
104 char cwcs, charwcs;
105 int iwcs;
106 char keyword[12];
107 char *upval, value[72];
108
109 /* If no WCS character, return 0 */
110 if (name == NULL)
111 return ((char) 0);
112
113 /* Convert input name to upper case */
114 upname = uppercase (name);
115
116 /* If single character name, return that character */
117 if (strlen (upname) == 1)
118 return (upname[0]);
119
120 /* Try to match input name to available WCSNAME names in header */
121 strcpy (keyword, "WCSNAME");
122 keyword[8] = (char) 0;
123 charwcs = '_';
124 for (iwcs = 0; iwcs < 27; iwcs++) {
125 if (iwcs > 0)
126 cwcs = (char) (64 + iwcs);
127 else
128 cwcs = (char) 0;
129 keyword[7] = cwcs;
130 if (hgets (hstring, keyword, 72, value)) {
131 upval = uppercase (value);
132 if (!strcmp (upval, upname))
133 charwcs = cwcs;
134 free (upval);
135 }
136 }
137 free (upname);
138 return (charwcs);
139 }
140
141
142 /* Make string of arbitrary case all uppercase */
143
144 char *
uppercase(string)145 uppercase (string)
146 const char *string;
147 {
148 int lstring, i;
149 char *upstring;
150
151 lstring = strlen (string);
152 upstring = (char *) calloc (1,lstring+1);
153 for (i = 0; i < lstring; i++) {
154 if (string[i] > 96 && string[i] < 123)
155 upstring[i] = string[i] - 32;
156 else
157 upstring[i] = string[i];
158 }
159 upstring[lstring] = (char) 0;
160 return (upstring);
161 }
162
163
164 /* set up a WCS structure from a FITS image header lhstring bytes long */
165
166 struct WorldCoor *
wcsninit(hstring,lhstring)167 wcsninit (hstring, lhstring)
168
169 const char *hstring; /* character string containing FITS header information
170 in the format <keyword>= <value> [/ <comment>] */
171 int lhstring; /* Length of FITS header in bytes */
172 {
173 char mchar; /* Suffix character for one of multiple WCS */
174 mchar = (char) 0;
175 hlength (hstring, lhstring);
176 return (wcsinitc (hstring, &mchar));
177 }
178
179
180 /* set up a WCS structure from a FITS image header lhstring bytes long */
181
182 struct WorldCoor *
wcsninitc(hstring,lhstring,mchar)183 wcsninitc (hstring, lhstring, mchar)
184
185 const char *hstring; /* character string containing FITS header information
186 in the format <keyword>= <value> [/ <comment>] */
187 int lhstring; /* Length of FITS header in bytes */
188 char *mchar; /* Suffix character for one of multiple WCS */
189 {
190 hlength (hstring, lhstring);
191 if (mchar[0] == ' ')
192 mchar[0] = (char) 0;
193 return (wcsinitc (hstring, mchar));
194 }
195
196
197 /* set up a WCS structure from a FITS image header */
198
199 struct WorldCoor *
wcsinit(hstring)200 wcsinit (hstring)
201
202 const char *hstring; /* character string containing FITS header information
203 in the format <keyword>= <value> [/ <comment>] */
204 {
205 char mchar; /* Suffix character for one of multiple WCS */
206 mchar = (char) 0;
207 return (wcsinitc (hstring, &mchar));
208 }
209
210
211 /* set up a WCS structure from a FITS image header for specified suffix */
212
213 struct WorldCoor *
wcsinitc(hstring,wchar)214 wcsinitc (hstring, wchar)
215
216 const char *hstring; /* character string containing FITS header information
217 in the format <keyword>= <value> [/ <comment>] */
218 char *wchar; /* Suffix character for one of multiple WCS */
219 {
220 struct WorldCoor *wcs, *depwcs;
221 char ctype1[32], ctype2[32], tstring[32];
222 char pvkey1[8],pvkey2[8],pvkey3[8];
223 char *hcoeff; /* pointer to first coeff's in header */
224 char decsign;
225 double rah,ram,ras, dsign,decd,decm,decs;
226 double dec_deg,ra_hours, secpix, ra0, ra1, dec0, dec1, cvel;
227 double cdelt1, cdelt2, cd[4], pc[81];
228 char keyword[16];
229 int ieq, i, j, k, naxes, cd11p, cd12p, cd21p, cd22p;
230 int ilat; /* coordinate for latitude or declination */
231 /*
232 int ix1, ix2, iy1, iy2, idx1, idx2, idy1, idy2;
233 double dxrefpix, dyrefpix;
234 */
235 char temp[80];
236 char wcsname[64]; /* Name of WCS depended on by current WCS */
237 char mchar;
238 char cspace = (char) ' ';
239 char cnull = (char) 0;
240 double mjd;
241 double rot;
242 double ut;
243 int nax;
244 int twod;
245 extern int tnxinit();
246 extern int zpxinit();
247 extern int platepos();
248 extern int dsspos();
249 void invert_wcs();
250
251 wcs = (struct WorldCoor *) calloc (1, sizeof(struct WorldCoor));
252
253 /* Set WCS character and name in structure */
254 mchar = wchar[0];
255 if (mchar == ' ')
256 mchar = cnull;
257 wcs->wcschar = mchar;
258 if (hgetsc (hstring, "WCSNAME", &mchar, 63, wcsname)) {
259 wcs->wcsname = (char *) calloc (strlen (wcsname)+2, 1);
260 strcpy (wcs->wcsname, wcsname);
261 }
262
263
264 /* Set WCSLIB flags so that structures will be reinitialized */
265 wcs->cel.flag = 0;
266 wcs->lin.flag = 0;
267 wcs->wcsl.flag = 0;
268 wcs->wcsl.cubeface = -1;
269
270 /* Initialize to no plate fit */
271 wcs->ncoeff1 = 0;
272 wcs->ncoeff2 = 0;
273
274 /* Initialize to no CD matrix */
275 cdelt1 = 0.0;
276 cdelt2 = 0.0;
277 cd[0] = 0.0;
278 cd[1] = 0.0;
279 cd[2] = 0.0;
280 cd[3] = 0.0;
281 pc[0] = 0.0;
282 wcs->rotmat = 0;
283 wcs->rot = 0.0;
284
285 /* Header parameters independent of projection */
286 naxes = 0;
287 hgeti4c (hstring, "WCSAXES", &mchar, &naxes);
288 if (naxes == 0)
289 hgeti4 (hstring, "WCSAXES", &naxes);
290 if (naxes == 0)
291 hgeti4 (hstring, "NAXIS", &naxes);
292 if (naxes == 0)
293 hgeti4 (hstring, "WCSDIM", &naxes);
294 if (naxes < 1) {
295 setwcserr ("WCSINIT: No WCSAXES, NAXIS, or WCSDIM keyword");
296 wcsfree (wcs);
297 return (NULL);
298 }
299 if (naxes > 2)
300 naxes = 2;
301 wcs->naxis = naxes;
302 wcs->naxes = naxes;
303 wcs->lin.naxis = naxes;
304 wcs->nxpix = 0;
305 hgetr8 (hstring, "NAXIS1", &wcs->nxpix);
306 if (wcs->nxpix < 1)
307 hgetr8 (hstring, "IMAGEW", &wcs->nxpix);
308 if (wcs->nxpix < 1) {
309 setwcserr ("WCSINIT: No NAXIS1 or IMAGEW keyword");
310 wcsfree (wcs);
311 return (NULL);
312 }
313 wcs->nypix = 0;
314 hgetr8 (hstring, "NAXIS2", &wcs->nypix);
315 if (wcs->nypix < 1)
316 hgetr8 (hstring, "IMAGEH", &wcs->nypix);
317 if (naxes > 1 && wcs->nypix < 1) {
318 setwcserr ("WCSINIT: No NAXIS2 or IMAGEH keyword");
319 wcsfree (wcs);
320 return (NULL);
321 }
322
323 /* Reset number of axes to only those with dimension greater than one */
324 nax = 0;
325 for (i = 0; i < naxes; i++) {
326
327 /* Check for number of pixels in axis more than one */
328 strcpy (keyword, "NAXIS");
329 sprintf (temp, "%d", i+1);
330 strcat (keyword, temp);
331 if (!hgeti4 (hstring, keyword, &j)) {
332 if (i == 0 && wcs->nxpix > 1) {
333 /* fprintf (stderr,"WCSINIT: Missing keyword %s set to %.0f from IMAGEW\n",
334 keyword, wcs->nxpix); */
335 j = wcs->nxpix;
336 }
337 else if (i == 1 && wcs->nypix > 1) {
338 /* fprintf (stderr,"WCSINIT: Missing keyword %s set to %.0f from IMAGEH\n",
339 keyword, wcs->nypix); */
340 j = wcs->nypix;
341 }
342 else
343 fprintf (stderr,"WCSINIT: Missing keyword %s assumed 1\n",keyword);
344 }
345
346 /* Check for TAB WCS in axis */
347 strcpy (keyword, "CTYPE");
348 strcat (keyword, temp);
349 if (hgets (hstring, keyword, 16, temp)) {
350 if (strsrch (temp, "-TAB"))
351 j = 0;
352 }
353 if (j > 1) nax = nax + 1;
354 }
355 naxes = nax;
356 wcs->naxes = nax;
357 wcs->naxis = nax;
358
359 hgets (hstring, "INSTRUME", 16, wcs->instrument);
360 hgeti4 (hstring, "DETECTOR", &wcs->detector);
361 wcs->wcsproj = getdefwcs();
362 wcs->logwcs = 0;
363 hgeti4 (hstring, "DC-FLAG", &wcs->logwcs);
364
365 /* Initialize rotation matrices */
366 for (i = 0; i < 81; i++) wcs->pc[i] = 0.0;
367 for (i = 0; i < 81; i++) pc[i] = 0.0;
368 for (i = 0; i < naxes; i++) wcs->pc[(i*naxes)+i] = 1.0;
369 for (i = 0; i < naxes; i++) pc[(i*naxes)+i] = 1.0;
370 for (i = 0; i < 9; i++) wcs->cdelt[i] = 0.0;
371 for (i = 0; i < naxes; i++) wcs->cdelt[i] = 1.0;
372
373 /* If the current world coordinate system depends on another, set it now */
374 if (hgetsc (hstring, "WCSDEP",&mchar, 63, wcsname)) {
375 if ((wcs->wcs = wcsinitn (hstring, wcsname)) == NULL) {
376 setwcserr ("WCSINIT: depended on WCS could not be set");
377 wcsfree (wcs);
378 return (NULL);
379 }
380 depwcs = wcs->wcs;
381 depwcs->wcsdep = wcs;
382 }
383 else
384 wcs->wcs = NULL;
385
386 /* Read radial velocity from image header */
387 wcs->radvel = 0.0;
388 wcs->zvel = 0.0;
389 cvel = 299792.5;
390 if (hgetr8c (hstring, "VSOURCE", &mchar, &wcs->radvel))
391 wcs->zvel = wcs->radvel / cvel;
392 else if (hgetr8c (hstring, "ZSOURCE", &mchar, &wcs->zvel))
393 wcs->radvel = wcs->zvel * cvel;
394 else if (hgetr8 (hstring, "VELOCITY", &wcs->radvel))
395 wcs->zvel = wcs->radvel / cvel;
396
397 for (i = 0; i < 10; i++) {
398 wcs->prj.p[i] = 0.0;
399 }
400
401 /* World coordinate system reference coordinate information */
402 if (hgetsc (hstring, "CTYPE1", &mchar, 16, ctype1)) {
403
404 /* Read second coordinate type */
405 strcpy (ctype2, ctype1);
406 if (!hgetsc (hstring, "CTYPE2", &mchar, 16, ctype2))
407 twod = 0;
408 else
409 twod = 1;
410 strcpy (wcs->ctype[0], ctype1);
411 strcpy (wcs->ctype[1], ctype2);
412 if (strsrch (ctype2, "LAT") || strsrch (ctype2, "DEC"))
413 ilat = 2;
414 else
415 ilat = 1;
416
417 /* Read third and fourth coordinate types, if present */
418 strcpy (wcs->ctype[2], "");
419 hgetsc (hstring, "CTYPE3", &mchar, 9, wcs->ctype[2]);
420 strcpy (wcs->ctype[3], "");
421 hgetsc (hstring, "CTYPE4", &mchar, 9, wcs->ctype[3]);
422
423 /* Set projection type in WCS data structure */
424 if (wcstype (wcs, ctype1, ctype2)) {
425 wcsfree (wcs);
426 return (NULL);
427 }
428
429 /* Get units, if present, for linear coordinates */
430 if (wcs->prjcode == WCS_LIN) {
431 if (!hgetsc (hstring, "CUNIT1", &mchar, 16, wcs->units[0])) {
432 if (!mgetstr (hstring, "WAT1", "units", 16, wcs->units[0])) {
433 wcs->units[0][0] = 0;
434 }
435 }
436 if (!strcmp (wcs->units[0], "pixel"))
437 wcs->prjcode = WCS_PIX;
438 if (twod) {
439 if (!hgetsc (hstring, "CUNIT2", &mchar, 16, wcs->units[1])) {
440 if (!mgetstr (hstring, "WAT2", "units", 16, wcs->units[1])) {
441 wcs->units[1][0] = 0;
442 }
443 }
444 if (!strcmp (wcs->units[0], "pixel"))
445 wcs->prjcode = WCS_PIX;
446 }
447 }
448
449 /* Reference pixel coordinates and WCS value */
450 wcs->crpix[0] = 1.0;
451 hgetr8c (hstring, "CRPIX1", &mchar, &wcs->crpix[0]);
452 wcs->crpix[1] = 1.0;
453 hgetr8c (hstring, "CRPIX2", &mchar, &wcs->crpix[1]);
454 wcs->xrefpix = wcs->crpix[0];
455 wcs->yrefpix = wcs->crpix[1];
456 wcs->crval[0] = 0.0;
457 hgetr8c (hstring, "CRVAL1", &mchar, &wcs->crval[0]);
458 wcs->crval[1] = 0.0;
459 hgetr8c (hstring, "CRVAL2", &mchar, &wcs->crval[1]);
460 if (wcs->syswcs == WCS_NPOLE)
461 wcs->crval[1] = 90.0 - wcs->crval[1];
462 if (wcs->syswcs == WCS_SPA)
463 wcs->crval[1] = wcs->crval[1] - 90.0;
464 wcs->xref = wcs->crval[0];
465 wcs->yref = wcs->crval[1];
466 if (wcs->coorflip) {
467 wcs->cel.ref[0] = wcs->crval[1];
468 wcs->cel.ref[1] = wcs->crval[0];
469 }
470 else {
471 wcs->cel.ref[0] = wcs->crval[0];
472 wcs->cel.ref[1] = wcs->crval[1];
473 }
474 wcs->longpole = 999.0;
475 hgetr8c (hstring, "LONPOLE", &mchar, &wcs->longpole);
476 wcs->cel.ref[2] = wcs->longpole;
477 wcs->latpole = 999.0;
478 hgetr8c (hstring, "LATPOLE", &mchar, &wcs->latpole);
479 wcs->cel.ref[3] = wcs->latpole;
480 wcs->lin.crpix = wcs->crpix;
481 wcs->lin.cdelt = wcs->cdelt;
482 wcs->lin.pc = wcs->pc;
483
484 /* Projection constants (this should be projection-dependent */
485 wcs->prj.r0 = 0.0;
486 hgetr8c (hstring, "PROJR0", &mchar, &wcs->prj.r0);
487
488 /* FITS WCS interim proposal projection constants */
489 for (i = 0; i < 10; i++) {
490 sprintf (keyword,"PROJP%d",i);
491 hgetr8c (hstring, keyword, &mchar, &wcs->prj.p[i]);
492 }
493
494 sprintf (pvkey1, "PV%d_1", ilat);
495 sprintf (pvkey2, "PV%d_2", ilat);
496 sprintf (pvkey3, "PV%d_3", ilat);
497
498 /* FITS WCS standard projection constants (projection-dependent) */
499 if (wcs->prjcode == WCS_AZP || wcs->prjcode == WCS_SIN ||
500 wcs->prjcode == WCS_COP || wcs->prjcode == WCS_COE ||
501 wcs->prjcode == WCS_COD || wcs->prjcode == WCS_COO) {
502 hgetr8c (hstring, pvkey1, &mchar, &wcs->prj.p[1]);
503 hgetr8c (hstring, pvkey2, &mchar, &wcs->prj.p[2]);
504 }
505 else if (wcs->prjcode == WCS_SZP) {
506 hgetr8c (hstring, pvkey1, &mchar, &wcs->prj.p[1]);
507 hgetr8c (hstring, pvkey2, &mchar, &wcs->prj.p[2]);
508 if (wcs->prj.p[3] == 0.0)
509 wcs->prj.p[3] = 90.0;
510 hgetr8c (hstring, pvkey3, &mchar, &wcs->prj.p[3]);
511 }
512 else if (wcs->prjcode == WCS_CEA) {
513 if (wcs->prj.p[1] == 0.0)
514 wcs->prj.p[1] = 1.0;
515 hgetr8c (hstring, pvkey1, &mchar, &wcs->prj.p[1]);
516 }
517 else if (wcs->prjcode == WCS_CYP) {
518 if (wcs->prj.p[1] == 0.0)
519 wcs->prj.p[1] = 1.0;
520 hgetr8c (hstring, pvkey1, &mchar, &wcs->prj.p[1]);
521 if (wcs->prj.p[2] == 0.0)
522 wcs->prj.p[2] = 1.0;
523 hgetr8c (hstring, pvkey2, &mchar, &wcs->prj.p[2]);
524 }
525 else if (wcs->prjcode == WCS_AIR) {
526 if (wcs->prj.p[1] == 0.0)
527 wcs->prj.p[1] = 90.0;
528 hgetr8c (hstring, pvkey1, &mchar, &wcs->prj.p[1]);
529 }
530 else if (wcs->prjcode == WCS_BON) {
531 hgetr8c (hstring, pvkey1, &mchar, &wcs->prj.p[1]);
532 }
533 else if (wcs->prjcode == WCS_ZPN) {
534 for (i = 0; i < 10; i++) {
535 sprintf (keyword,"PV%d_%d", ilat, i);
536 hgetr8c (hstring, keyword, &mchar, &wcs->prj.p[i]);
537 }
538 }
539
540 /* Initialize TNX, defaulting to TAN if there is a problem */
541 if (wcs->prjcode == WCS_TNX) {
542 if (tnxinit (hstring, wcs)) {
543 wcs->ctype[0][6] = 'A';
544 wcs->ctype[0][7] = 'N';
545 wcs->ctype[1][6] = 'A';
546 wcs->ctype[1][7] = 'N';
547 wcs->prjcode = WCS_TAN;
548 }
549 }
550
551 /* Initialize ZPX, defaulting to ZPN if there is a problem */
552 if (wcs->prjcode == WCS_ZPX) {
553 if (zpxinit (hstring, wcs)) {
554 wcs->ctype[0][7] = 'N';
555 wcs->ctype[1][7] = 'N';
556 wcs->prjcode = WCS_ZPN;
557 }
558 }
559
560 /* Set TPV to TAN as SCAMP coefficients will be added below */
561 if (wcs->prjcode == WCS_TPV) {
562 wcs->ctype[0][6] = 'A';
563 wcs->ctype[0][7] = 'N';
564 wcs->ctype[1][6] = 'A';
565 wcs->ctype[1][7] = 'N';
566 wcs->prjcode = WCS_TAN;
567 }
568
569 /* Coordinate reference frame, equinox, and epoch */
570 if (wcs->wcsproj > 0)
571 wcseqm (hstring, wcs, &mchar);
572 wcsioset (wcs);
573
574 /* Read distortion coefficients, if present */
575 distortinit (wcs, hstring);
576
577 /* Use polynomial fit instead of projection, if present */
578 wcs->ncoeff1 = 0;
579 wcs->ncoeff2 = 0;
580 cd11p = hgetr8c (hstring, "CD1_1", &mchar, &cd[0]);
581 cd12p = hgetr8c (hstring, "CD1_2", &mchar, &cd[1]);
582 cd21p = hgetr8c (hstring, "CD2_1", &mchar, &cd[2]);
583 cd22p = hgetr8c (hstring, "CD2_2", &mchar, &cd[3]);
584 if (wcs->wcsproj != WCS_OLD &&
585 (hcoeff = ksearch (hstring,"CO1_1")) != NULL) {
586 wcs->prjcode = WCS_PLT;
587 (void)strcpy (wcs->ptype, "PLA");
588 for (i = 0; i < 20; i++) {
589 sprintf (keyword,"CO1_%d", i+1);
590 wcs->x_coeff[i] = 0.0;
591 if (hgetr8 (hcoeff, keyword, &wcs->x_coeff[i]))
592 wcs->ncoeff1 = i + 1;
593 }
594 hcoeff = ksearch (hstring,"CO2_1");
595 for (i = 0; i < 20; i++) {
596 sprintf (keyword,"CO2_%d",i+1);
597 wcs->y_coeff[i] = 0.0;
598 if (hgetr8 (hcoeff, keyword, &wcs->y_coeff[i]))
599 wcs->ncoeff2 = i + 1;
600 }
601
602 /* Compute a nominal scale factor */
603 platepos (wcs->crpix[0], wcs->crpix[1], wcs, &ra0, &dec0);
604 platepos (wcs->crpix[0], wcs->crpix[1]+1.0, wcs, &ra1, &dec1);
605 wcs->yinc = dec1 - dec0;
606 wcs->xinc = -wcs->yinc;
607
608 /* Compute image rotation angle */
609 wcs->wcson = 1;
610 wcsrotset (wcs);
611 rot = degrad (wcs->rot);
612
613 /* Compute scale at reference pixel */
614 platepos (wcs->crpix[0], wcs->crpix[1], wcs, &ra0, &dec0);
615 platepos (wcs->crpix[0]+cos(rot),
616 wcs->crpix[1]+sin(rot), wcs, &ra1, &dec1);
617 wcs->cdelt[0] = -wcsdist (ra0, dec0, ra1, dec1);
618 wcs->xinc = wcs->cdelt[0];
619 platepos (wcs->crpix[0]+sin(rot),
620 wcs->crpix[1]+cos(rot), wcs, &ra1, &dec1);
621 wcs->cdelt[1] = wcsdist (ra0, dec0, ra1, dec1);
622 wcs->yinc = wcs->cdelt[1];
623
624 /* Set CD matrix from header */
625 wcs->cd[0] = cd[0];
626 wcs->cd[1] = cd[1];
627 wcs->cd[2] = cd[2];
628 wcs->cd[3] = cd[3];
629 (void) matinv (2, wcs->cd, wcs->dc);
630 }
631
632 /* Else use CD matrix, if present */
633 else if (cd11p || cd12p || cd21p || cd22p) {
634 wcs->rotmat = 1;
635 wcscdset (wcs, cd);
636 }
637
638 /* Else get scaling from CDELT1 and CDELT2 */
639 else if (hgetr8c (hstring, "CDELT1", &mchar, &cdelt1) != 0) {
640 hgetr8c (hstring, "CDELT2", &mchar, &cdelt2);
641
642 /* If CDELT1 or CDELT2 is 0 or missing */
643 if (cdelt1 == 0.0 || (wcs->nypix > 1 && cdelt2 == 0.0)) {
644 if (ksearch (hstring,"SECPIX") != NULL ||
645 ksearch (hstring,"PIXSCALE") != NULL ||
646 ksearch (hstring,"PIXSCAL1") != NULL ||
647 ksearch (hstring,"XPIXSIZE") != NULL ||
648 ksearch (hstring,"SECPIX1") != NULL) {
649 secpix = 0.0;
650 hgetr8 (hstring,"SECPIX",&secpix);
651 if (secpix == 0.0)
652 hgetr8 (hstring,"PIXSCALE",&secpix);
653 if (secpix == 0.0) {
654 hgetr8 (hstring,"SECPIX1",&secpix);
655 if (secpix != 0.0) {
656 if (cdelt1 == 0.0)
657 cdelt1 = -secpix / 3600.0;
658 if (cdelt2 == 0.0) {
659 hgetr8 (hstring,"SECPIX2",&secpix);
660 cdelt2 = secpix / 3600.0;
661 }
662 }
663 else {
664 hgetr8 (hstring,"XPIXSIZE",&secpix);
665 if (secpix != 0.0) {
666 if (cdelt1 == 0.0)
667 cdelt1 = -secpix / 3600.0;
668 if (cdelt2 == 0.0) {
669 hgetr8 (hstring,"YPIXSIZE",&secpix);
670 cdelt2 = secpix / 3600.0;
671 }
672 }
673 else {
674 hgetr8 (hstring,"PIXSCAL1",&secpix);
675 if (secpix != 0.0 && cdelt1 == 0.0)
676 cdelt1 = -secpix / 3600.0;
677 if (cdelt2 == 0.0) {
678 hgetr8 (hstring,"PIXSCAL2",&secpix);
679 cdelt2 = secpix / 3600.0;
680 }
681 }
682 }
683 }
684 else {
685 if (cdelt1 == 0.0)
686 cdelt1 = -secpix / 3600.0;
687 if (cdelt2 == 0.0)
688 cdelt2 = secpix / 3600.0;
689 }
690 }
691 }
692 if (cdelt2 == 0.0 && wcs->nypix > 1)
693 cdelt2 = -cdelt1;
694 wcs->cdelt[2] = 1.0;
695 wcs->cdelt[3] = 1.0;
696
697 /* Initialize rotation matrix */
698 for (i = 0; i < 81; i++) {
699 pc[i] = 0.0;
700 wcs->pc[i] = 0.0;
701 }
702 for (i = 0; i < naxes; i++)
703 pc[(i*naxes)+i] = 1.0;
704
705 /* Read FITS WCS interim rotation matrix */
706 if (!mchar && hgetr8 (hstring,"PC001001",&pc[0]) != 0) {
707 k = 0;
708 for (i = 0; i < naxes; i++) {
709 for (j = 0; j < naxes; j++) {
710 if (i == j)
711 pc[k] = 1.0;
712 else
713 pc[k] = 0.0;
714 sprintf (keyword, "PC00%1d00%1d", i+1, j+1);
715 hgetr8 (hstring, keyword, &pc[k++]);
716 }
717 }
718 wcspcset (wcs, cdelt1, cdelt2, pc);
719 }
720
721 /* Read FITS WCS standard rotation matrix */
722 else if (hgetr8c (hstring, "PC1_1", &mchar, &pc[0]) != 0) {
723 k = 0;
724 for (i = 0; i < naxes; i++) {
725 for (j = 0; j < naxes; j++) {
726 if (i == j)
727 pc[k] = 1.0;
728 else
729 pc[k] = 0.0;
730 sprintf (keyword, "PC%1d_%1d", i+1, j+1);
731 hgetr8c (hstring, keyword, &mchar, &pc[k++]);
732 }
733 }
734 wcspcset (wcs, cdelt1, cdelt2, pc);
735 }
736
737 /* Otherwise, use CROTAn */
738 else {
739 rot = 0.0;
740 if (ilat == 2)
741 hgetr8c (hstring, "CROTA2", &mchar, &rot);
742 else
743 hgetr8c (hstring,"CROTA1", &mchar, &rot);
744 wcsdeltset (wcs, cdelt1, cdelt2, rot);
745 }
746 }
747
748 /* If no scaling is present, set to 1 per pixel, no rotation */
749 else {
750 wcs->xinc = 1.0;
751 wcs->yinc = 1.0;
752 wcs->cdelt[0] = 1.0;
753 wcs->cdelt[1] = 1.0;
754 wcs->rot = 0.0;
755 wcs->rotmat = 0;
756 setwcserr ("WCSINIT: setting CDELT to 1");
757 }
758
759 /* SCAMP convention */
760 if (wcs->prjcode == WCS_TAN && wcs->naxis == 2) {
761 int n = 0;
762 if (wcs->inv_x) {
763 poly_end(wcs->inv_x);
764 wcs->inv_x = NULL;
765 }
766 if (wcs->inv_y) {
767 poly_end(wcs->inv_y);
768 wcs->inv_y = NULL;
769 }
770 wcs->pvfail = 0;
771 for (i = 0; i < (2*MAXPV); i++) {
772 wcs->projppv[i] = 0.0;
773 wcs->prj.ppv[i] = 0.0;
774 }
775 for (k = 0; k < 2; k++) {
776 for (j = 0; j < MAXPV; j++) {
777 sprintf(keyword, "PV%d_%d", k+1, j);
778 if (hgetr8c(hstring, keyword,&mchar, &wcs->projppv[j+k*MAXPV]) == 0) {
779 wcs->projppv[j+k*MAXPV] = 0.0;
780 }
781 else
782 n++;
783 }
784 }
785
786 /* If any PVi_j are set, add them in the structure if no SIRTF distortion*/
787 if (n > 0 && wcs->distcode != DISTORT_SIRTF) {
788 n = 0;
789
790 for (k = MAXPV; k >= 0; k--) {
791 /* lat comes first for compatibility reasons */
792 wcs->prj.ppv[k] = wcs->projppv[k+wcs->wcsl.lat*MAXPV];
793 wcs->prj.ppv[k+MAXPV] = wcs->projppv[k+wcs->wcsl.lng*MAXPV];
794 if (!n && (wcs->prj.ppv[k] || wcs->prj.ppv[k+MAXPV])) {
795 n = k+1;
796 }
797 }
798 invert_wcs(wcs);
799
800 /* Need to call tanset again */
801 wcs->cel.flag = 0;
802 }
803 }
804
805 /* If linear or pixel WCS, print "degrees" */
806 if (!strncmp (wcs->ptype,"LIN",3) ||
807 !strncmp (wcs->ptype,"PIX",3)) {
808 wcs->degout = -1;
809 wcs->ndec = 5;
810 }
811
812 /* Epoch of image (from observation date, if possible) */
813 if (hgetr8 (hstring, "MJD-OBS", &mjd))
814 wcs->epoch = 1900.0 + (mjd - 15019.81352) / 365.242198781;
815 else if (!hgetdate (hstring,"DATE-OBS",&wcs->epoch)) {
816 if (!hgetdate (hstring,"DATE",&wcs->epoch)) {
817 if (!hgetr8 (hstring,"EPOCH",&wcs->epoch))
818 wcs->epoch = wcs->equinox;
819 }
820 }
821
822 /* Add time of day if not part of DATE-OBS string */
823 else {
824 hgets (hstring,"DATE-OBS",32,tstring);
825 if (!strchr (tstring,'T')) {
826 if (hgetr8 (hstring, "UT",&ut))
827 wcs->epoch = wcs->epoch + (ut / (24.0 * 365.242198781));
828 else if (hgetr8 (hstring, "UTMID",&ut))
829 wcs->epoch = wcs->epoch + (ut / (24.0 * 365.242198781));
830 }
831 }
832
833 wcs->wcson = 1;
834 }
835
836 else if (mchar != cnull && mchar != cspace) {
837 (void) sprintf (temp, "WCSINITC: No image scale for WCS %c", mchar);
838 setwcserr (temp);
839 wcsfree (wcs);
840 return (NULL);
841 }
842
843 /* Plate solution coefficients */
844 else if (ksearch (hstring,"PLTRAH") != NULL) {
845 wcs->prjcode = WCS_DSS;
846 hcoeff = ksearch (hstring,"PLTRAH");
847 hgetr8 (hcoeff,"PLTRAH",&rah);
848 hgetr8 (hcoeff,"PLTRAM",&ram);
849 hgetr8 (hcoeff,"PLTRAS",&ras);
850 ra_hours = rah + (ram / (double)60.0) + (ras / (double)3600.0);
851 wcs->plate_ra = hrrad (ra_hours);
852 decsign = '+';
853 hgets (hcoeff,"PLTDECSN", 1, &decsign);
854 if (decsign == '-')
855 dsign = -1.;
856 else
857 dsign = 1.;
858 hgetr8 (hcoeff,"PLTDECD",&decd);
859 hgetr8 (hcoeff,"PLTDECM",&decm);
860 hgetr8 (hcoeff,"PLTDECS",&decs);
861 dec_deg = dsign * (decd+(decm/(double)60.0)+(decs/(double)3600.0));
862 wcs->plate_dec = degrad (dec_deg);
863 hgetr8 (hstring,"EQUINOX",&wcs->equinox);
864 hgeti4 (hstring,"EQUINOX",&ieq);
865 if (ieq == 1950)
866 strcpy (wcs->radecsys,"FK4");
867 else
868 strcpy (wcs->radecsys,"FK5");
869 wcs->epoch = wcs->equinox;
870 hgetr8 (hstring,"EPOCH",&wcs->epoch);
871 (void)sprintf (wcs->center,"%2.0f:%2.0f:%5.3f %c%2.0f:%2.0f:%5.3f %s",
872 rah,ram,ras,decsign,decd,decm,decs,wcs->radecsys);
873 hgetr8 (hstring,"PLTSCALE",&wcs->plate_scale);
874 hgetr8 (hstring,"XPIXELSZ",&wcs->x_pixel_size);
875 hgetr8 (hstring,"YPIXELSZ",&wcs->y_pixel_size);
876 hgetr8 (hstring,"CNPIX1",&wcs->x_pixel_offset);
877 hgetr8 (hstring,"CNPIX2",&wcs->y_pixel_offset);
878 hcoeff = ksearch (hstring,"PPO1");
879 for (i = 0; i < 6; i++) {
880 sprintf (keyword,"PPO%d", i+1);
881 wcs->ppo_coeff[i] = 0.0;
882 hgetr8 (hcoeff,keyword,&wcs->ppo_coeff[i]);
883 }
884 hcoeff = ksearch (hstring,"AMDX1");
885 for (i = 0; i < 20; i++) {
886 sprintf (keyword,"AMDX%d", i+1);
887 wcs->x_coeff[i] = 0.0;
888 hgetr8 (hcoeff, keyword, &wcs->x_coeff[i]);
889 }
890 hcoeff = ksearch (hstring,"AMDY1");
891 for (i = 0; i < 20; i++) {
892 sprintf (keyword,"AMDY%d",i+1);
893 wcs->y_coeff[i] = 0.0;
894 hgetr8 (hcoeff, keyword, &wcs->y_coeff[i]);
895 }
896 wcs->wcson = 1;
897 (void)strcpy (wcs->c1type, "RA");
898 (void)strcpy (wcs->c2type, "DEC");
899 (void)strcpy (wcs->ptype, "DSS");
900 wcs->degout = 0;
901 wcs->ndec = 3;
902
903 /* Compute a nominal reference pixel at the image center */
904 strcpy (wcs->ctype[0], "RA---DSS");
905 strcpy (wcs->ctype[1], "DEC--DSS");
906 wcs->crpix[0] = 0.5 * wcs->nxpix;
907 wcs->crpix[1] = 0.5 * wcs->nypix;
908 wcs->xrefpix = wcs->crpix[0];
909 wcs->yrefpix = wcs->crpix[1];
910 dsspos (wcs->crpix[0], wcs->crpix[1], wcs, &ra0, &dec0);
911 wcs->crval[0] = ra0;
912 wcs->crval[1] = dec0;
913 wcs->xref = wcs->crval[0];
914 wcs->yref = wcs->crval[1];
915
916 /* Compute a nominal scale factor */
917 dsspos (wcs->crpix[0], wcs->crpix[1]+1.0, wcs, &ra1, &dec1);
918 wcs->yinc = dec1 - dec0;
919 wcs->xinc = -wcs->yinc;
920 wcsioset (wcs);
921
922 /* Compute image rotation angle */
923 wcs->wcson = 1;
924 wcsrotset (wcs);
925 rot = degrad (wcs->rot);
926
927 /* Compute image scale at center */
928 dsspos (wcs->crpix[0]+cos(rot),
929 wcs->crpix[1]+sin(rot), wcs, &ra1, &dec1);
930 wcs->cdelt[0] = -wcsdist (ra0, dec0, ra1, dec1);
931 dsspos (wcs->crpix[0]+sin(rot),
932 wcs->crpix[1]+cos(rot), wcs, &ra1, &dec1);
933 wcs->cdelt[1] = wcsdist (ra0, dec0, ra1, dec1);
934
935 /* Set all other image scale parameters */
936 wcsdeltset (wcs, wcs->cdelt[0], wcs->cdelt[1], wcs->rot);
937 }
938
939 /* Approximate world coordinate system if plate scale is known */
940 else if ((ksearch (hstring,"SECPIX") != NULL ||
941 ksearch (hstring,"PIXSCALE") != NULL ||
942 ksearch (hstring,"PIXSCAL1") != NULL ||
943 ksearch (hstring,"XPIXSIZE") != NULL ||
944 ksearch (hstring,"SECPIX1") != NULL)) {
945 secpix = 0.0;
946 hgetr8 (hstring,"SECPIX",&secpix);
947 if (secpix == 0.0)
948 hgetr8 (hstring,"PIXSCALE",&secpix);
949 if (secpix == 0.0) {
950 hgetr8 (hstring,"SECPIX1",&secpix);
951 if (secpix != 0.0) {
952 cdelt1 = -secpix / 3600.0;
953 hgetr8 (hstring,"SECPIX2",&secpix);
954 cdelt2 = secpix / 3600.0;
955 }
956 else {
957 hgetr8 (hstring,"XPIXSIZE",&secpix);
958 if (secpix != 0.0) {
959 cdelt1 = -secpix / 3600.0;
960 hgetr8 (hstring,"YPIXSIZE",&secpix);
961 cdelt2 = secpix / 3600.0;
962 }
963 else {
964 hgetr8 (hstring,"PIXSCAL1",&secpix);
965 cdelt1 = -secpix / 3600.0;
966 hgetr8 (hstring,"PIXSCAL2",&secpix);
967 cdelt2 = secpix / 3600.0;
968 }
969 }
970 }
971 else {
972 cdelt2 = secpix / 3600.0;
973 cdelt1 = -cdelt2;
974 }
975
976 /* Get rotation angle from the header, if it's there */
977 rot = 0.0;
978 hgetr8 (hstring,"CROTA1", &rot);
979 if (wcs->rot == 0.)
980 hgetr8 (hstring,"CROTA2", &rot);
981
982 /* Set CD and PC matrices */
983 wcsdeltset (wcs, cdelt1, cdelt2, rot);
984
985 /* By default, set reference pixel to center of image */
986 wcs->crpix[0] = 0.5 + (wcs->nxpix * 0.5);
987 wcs->crpix[1] = 0.5 + (wcs->nypix * 0.5);
988
989 /* Get reference pixel from the header, if it's there */
990 if (ksearch (hstring,"CRPIX1") != NULL) {
991 hgetr8 (hstring,"CRPIX1",&wcs->crpix[0]);
992 hgetr8 (hstring,"CRPIX2",&wcs->crpix[1]);
993 }
994
995 /* Use center of detector array as reference pixel
996 else if (ksearch (hstring,"DETSIZE") != NULL ||
997 ksearch (hstring,"DETSEC") != NULL) {
998 char *ic;
999 hgets (hstring, "DETSIZE", 32, temp);
1000 ic = strchr (temp, ':');
1001 if (ic != NULL)
1002 *ic = ' ';
1003 ic = strchr (temp, ',');
1004 if (ic != NULL)
1005 *ic = ' ';
1006 ic = strchr (temp, ':');
1007 if (ic != NULL)
1008 *ic = ' ';
1009 ic = strchr (temp, ']');
1010 if (ic != NULL)
1011 *ic = cnull;
1012 sscanf (temp, "%d %d %d %d", &idx1, &idx2, &idy1, &idy2);
1013 dxrefpix = 0.5 * (double) (idx1 + idx2 - 1);
1014 dyrefpix = 0.5 * (double) (idy1 + idy2 - 1);
1015 hgets (hstring, "DETSEC", 32, temp);
1016 ic = strchr (temp, ':');
1017 if (ic != NULL)
1018 *ic = ' ';
1019 ic = strchr (temp, ',');
1020 if (ic != NULL)
1021 *ic = ' ';
1022 ic = strchr (temp, ':');
1023 if (ic != NULL)
1024 *ic = ' ';
1025 ic = strchr (temp, ']');
1026 if (ic != NULL)
1027 *ic = cnull;
1028 sscanf (temp, "%d %d %d %d", &ix1, &ix2, &iy1, &iy2);
1029 wcs->crpix[0] = dxrefpix - (double) (ix1 - 1);
1030 wcs->crpix[1] = dyrefpix - (double) (iy1 - 1);
1031 } */
1032 wcs->xrefpix = wcs->crpix[0];
1033 wcs->yrefpix = wcs->crpix[1];
1034
1035 wcs->crval[0] = -999.0;
1036 if (!hgetra (hstring,"RA",&wcs->crval[0])) {
1037 setwcserr ("WCSINIT: No RA with SECPIX, no WCS");
1038 wcsfree (wcs);
1039 return (NULL);
1040 }
1041 wcs->crval[1] = -999.0;
1042 if (!hgetdec (hstring,"DEC",&wcs->crval[1])) {
1043 setwcserr ("WCSINIT No DEC with SECPIX, no WCS");
1044 wcsfree (wcs);
1045 return (NULL);
1046 }
1047 wcs->xref = wcs->crval[0];
1048 wcs->yref = wcs->crval[1];
1049 wcs->coorflip = 0;
1050
1051 wcs->cel.ref[0] = wcs->crval[0];
1052 wcs->cel.ref[1] = wcs->crval[1];
1053 wcs->cel.ref[2] = 999.0;
1054 if (!hgetr8 (hstring,"LONPOLE",&wcs->cel.ref[2]))
1055 hgetr8 (hstring,"LONGPOLE",&wcs->cel.ref[2]);
1056 wcs->cel.ref[3] = 999.0;
1057 hgetr8 (hstring,"LATPOLE",&wcs->cel.ref[3]);
1058
1059 /* Epoch of image (from observation date, if possible) */
1060 if (hgetr8 (hstring, "MJD-OBS", &mjd))
1061 wcs->epoch = 1900.0 + (mjd - 15019.81352) / 365.242198781;
1062 else if (!hgetdate (hstring,"DATE-OBS",&wcs->epoch)) {
1063 if (!hgetdate (hstring,"DATE",&wcs->epoch)) {
1064 if (!hgetr8 (hstring,"EPOCH",&wcs->epoch))
1065 wcs->epoch = wcs->equinox;
1066 }
1067 }
1068
1069 /* Add time of day if not part of DATE-OBS string */
1070 else {
1071 hgets (hstring,"DATE-OBS",32,tstring);
1072 if (!strchr (tstring,'T')) {
1073 if (hgetr8 (hstring, "UT",&ut))
1074 wcs->epoch = wcs->epoch + (ut / (24.0 * 365.242198781));
1075 else if (hgetr8 (hstring, "UTMID",&ut))
1076 wcs->epoch = wcs->epoch + (ut / (24.0 * 365.242198781));
1077 }
1078 }
1079
1080 /* Coordinate reference frame and equinox */
1081 (void) wcstype (wcs, "RA---TAN", "DEC--TAN");
1082 wcs->coorflip = 0;
1083 wcseq (hstring,wcs);
1084 wcsioset (wcs);
1085 wcs->degout = 0;
1086 wcs->ndec = 3;
1087 wcs->wcson = 1;
1088 }
1089
1090 else {
1091 setwcserr ("WCSINIT: No image scale");
1092 wcsfree (wcs);
1093 return (NULL);
1094 }
1095
1096 wcs->lin.crpix = wcs->crpix;
1097 wcs->lin.cdelt = wcs->cdelt;
1098 wcs->lin.pc = wcs->pc;
1099
1100 wcs->printsys = 1;
1101 wcs->tabsys = 0;
1102 wcs->linmode = 0;
1103
1104 /* Initialize special WCS commands */
1105 setwcscom (wcs);
1106
1107 return (wcs);
1108 }
1109
1110
1111 /******* invert_wcs ***********************************************************
1112 PROTO void invert_wcs(wcsstruct *wcs)
1113 PURPOSE Invert WCS projection mapping (using a polynomial).
1114 INPUT WCS structure.
1115 OUTPUT -.
1116 NOTES .
1117 AUTHOR E. Bertin (IAP)
1118 VERSION 06/11/2003
1119 ***/
1120
1121 void
invert_wcs(struct WorldCoor * wcs)1122 invert_wcs( struct WorldCoor *wcs)
1123
1124 {
1125 polystruct *poly;
1126 double pixin[NAXISPV],raw[NAXISPV],rawmin[NAXISPV];
1127 double *outpos,*outpost, *lngpos,*lngpost;
1128 double *latpos,*latpost,lngstep,latstep, rawsize, epsilon;
1129 int group[] = {1,1};
1130 /* Don't ask, this is needed by poly_init()! */
1131 int i,j,lng,lat,deg, maxflag;
1132 char errstr[80];
1133 double xmin;
1134 double ymin;
1135 double xmax;
1136 double ymax;
1137 double lngmin;
1138 double latmin;
1139
1140 /* Check first that inversion is not straightforward */
1141 lng = wcs->wcsl.lng;
1142 lat = wcs->wcsl.lat;
1143
1144 if (wcs->naxis != NAXISPV) {
1145 return;
1146 }
1147
1148 if (strcmp(wcs->wcsl.pcode, "TAN") != 0) {
1149 return;
1150 }
1151
1152 if ((wcs->projppv[1+lng*MAXPV] == 0) &&
1153 (wcs->projppv[1+lat*MAXPV] == 0)) {
1154 return;
1155 }
1156
1157 if (wcs->wcs != NULL) {
1158 pix2wcs(wcs->wcs,0,0,&xmin,&ymin);
1159 pix2wcs(wcs->wcs,wcs->nxpix,wcs->nypix,&xmax,&ymax);
1160 }
1161 else {
1162 xmin = 0;
1163 ymin = 0;
1164 xmax = wcs->nxpix;
1165 ymax = wcs->nypix;
1166 }
1167
1168 /* We define x as "longitude" and y as "latitude" projections */
1169 /* We assume that PCxx cross-terms with additional dimensions are small */
1170 /* Sample the whole image with a regular grid */
1171 if (lng == 0) {
1172 lngstep = (xmax-xmin)/(WCS_NGRIDPOINTS-1.0);
1173 lngmin = xmin;
1174 latstep = (ymax-ymin)/(WCS_NGRIDPOINTS-1.0);
1175 latmin = ymin;
1176 }
1177 else {
1178 lngstep = (ymax-ymin)/(WCS_NGRIDPOINTS-1.0);
1179 lngmin = ymin;
1180 latstep = (xmax-xmin)/(WCS_NGRIDPOINTS-1.0);
1181 latmin = xmin;
1182 }
1183
1184 outpos = (double *)calloc(2*WCS_NGRIDPOINTS2,sizeof(double));
1185 lngpos = (double *)calloc(WCS_NGRIDPOINTS2,sizeof(double));
1186 latpos = (double *)calloc(WCS_NGRIDPOINTS2,sizeof(double));
1187 raw[lat] = rawmin[lat] = 0.5+latmin;
1188 raw[lng] = rawmin[lng] = 0.5+lngmin;
1189 outpost = outpos;
1190 lngpost = lngpos;
1191 latpost = latpos;
1192 for (j=WCS_NGRIDPOINTS; j--; raw[lat]+=latstep) {
1193 raw[lng] = rawmin[lng];
1194 for (i=WCS_NGRIDPOINTS; i--; raw[lng]+=lngstep) {
1195 if (linrev(raw, &wcs->lin, pixin)) {
1196 sprintf (errstr,"*Error*: incorrect linear conversion in %s",
1197 wcs->wcsl.pcode);
1198 setwcserr (errstr);
1199 }
1200 *(lngpost++) = pixin[lng];
1201 *(latpost++) = pixin[lat];
1202 raw_to_pv (&wcs->prj,pixin[lng],pixin[lat], outpost, outpost+1);
1203 outpost += 2;
1204 }
1205 }
1206
1207 /* Invert "longitude" */
1208 /* Compute the extent of the pixel in reduced projected coordinates */
1209 linrev(rawmin, &wcs->lin, pixin);
1210 pixin[lng] += S2D;
1211 linfwd(pixin, &wcs->lin, raw);
1212 rawsize = sqrt((raw[lng]-rawmin[lng])*(raw[lng]-rawmin[lng])
1213 +(raw[lat]-rawmin[lat])*(raw[lat]-rawmin[lat]))*D2S;
1214 if (!rawsize) {
1215 sprintf (errstr,"*Error*: incorrect linear conversion in %s",
1216 wcs->wcsl.pcode);
1217 setwcserr (errstr);
1218 }
1219 epsilon = WCS_INVACCURACY/rawsize;
1220
1221 /* Find the lowest degree polynom */
1222 poly = NULL; /* to avoid gcc -Wall warnings */
1223 maxflag = 1;
1224 for (deg=1; deg<=WCS_INVMAXDEG && maxflag; deg++) {
1225 if (deg>1) {
1226 poly_end(poly);
1227 }
1228 poly = poly_init(group, 2, °, 1);
1229 poly_fit(poly, outpos, lngpos, NULL, WCS_NGRIDPOINTS2, NULL);
1230 maxflag = 0;
1231 outpost = outpos;
1232 lngpost = lngpos;
1233 for (i=WCS_NGRIDPOINTS2; i--; outpost+=2) {
1234 if (fabs(poly_func(poly, outpost)-*(lngpost++))>epsilon) {
1235 maxflag = 1;
1236 break;
1237 }
1238 }
1239 }
1240 if (maxflag) {
1241 setwcserr ("WARNING: Significant inaccuracy likely to occur in projection");
1242 wcs->pvfail = 1;
1243 }
1244
1245 /* Now link the created structure */
1246 wcs->prj.inv_x = wcs->inv_x = poly;
1247
1248 /* Invert "latitude" */
1249 /* Compute the extent of the pixel in reduced projected coordinates */
1250 linrev(rawmin, &wcs->lin, pixin);
1251 pixin[lat] += S2D;
1252 linfwd(pixin, &wcs->lin, raw);
1253 rawsize = sqrt((raw[lng]-rawmin[lng])*(raw[lng]-rawmin[lng])
1254 +(raw[lat]-rawmin[lat])*(raw[lat]-rawmin[lat]))*D2S;
1255 if (!rawsize) {
1256 sprintf (errstr,"*Error*: incorrect linear conversion in %s",
1257 wcs->wcsl.pcode);
1258 setwcserr (errstr);
1259 }
1260 epsilon = WCS_INVACCURACY/rawsize;
1261
1262 /* Find the lowest degree polynom */
1263 maxflag = 1;
1264 for (deg=1; deg<=WCS_INVMAXDEG && maxflag; deg++) {
1265 if (deg>1)
1266 poly_end(poly);
1267 poly = poly_init(group, 2, °, 1);
1268 poly_fit(poly, outpos, latpos, NULL, WCS_NGRIDPOINTS2, NULL);
1269 maxflag = 0;
1270 outpost = outpos;
1271 latpost = latpos;
1272 for (i=WCS_NGRIDPOINTS2; i--; outpost+=2) {
1273 if (fabs(poly_func(poly, outpost)-*(latpost++))>epsilon) {
1274 maxflag = 1;
1275 break;
1276 }
1277 }
1278 }
1279 if (maxflag) {
1280 setwcserr ("WARNING: Significant inaccuracy likely to occur in projection");
1281 wcs->pvfail = 1;
1282 }
1283
1284 /* Now link the created structure */
1285 wcs->prj.inv_y = wcs->inv_y = poly;
1286
1287 /* Free memory */
1288 free(outpos);
1289 free(lngpos);
1290 free(latpos);
1291
1292 return;
1293 }
1294
1295
1296 /* Set coordinate system of image, input, and output */
1297
1298 static void
wcsioset(wcs)1299 wcsioset (wcs)
1300
1301 struct WorldCoor *wcs;
1302 {
1303 if (strlen (wcs->radecsys) == 0 || wcs->prjcode == WCS_LIN)
1304 strcpy (wcs->radecsys, "LINEAR");
1305 if (wcs->prjcode == WCS_PIX)
1306 strcpy (wcs->radecsys, "PIXEL");
1307 wcs->syswcs = wcscsys (wcs->radecsys);
1308
1309 if (wcs->syswcs == WCS_B1950)
1310 strcpy (wcs->radecout, "FK4");
1311 else if (wcs->syswcs == WCS_J2000)
1312 strcpy (wcs->radecout, "FK5");
1313 else
1314 strcpy (wcs->radecout, wcs->radecsys);
1315 wcs->sysout = wcscsys (wcs->radecout);
1316 wcs->eqout = wcs->equinox;
1317 strcpy (wcs->radecin, wcs->radecsys);
1318 wcs->sysin = wcscsys (wcs->radecin);
1319 wcs->eqin = wcs->equinox;
1320 return;
1321 }
1322
1323
1324 static void
wcseq(hstring,wcs)1325 wcseq (hstring, wcs)
1326
1327 char *hstring; /* character string containing FITS header information
1328 in the format <keyword>= <value> [/ <comment>] */
1329 struct WorldCoor *wcs; /* World coordinate system data structure */
1330 {
1331 char mchar; /* Suffix character for one of multiple WCS */
1332 mchar = (char) 0;
1333 wcseqm (hstring, wcs, &mchar);
1334 return;
1335 }
1336
1337
1338 static void
wcseqm(hstring,wcs,mchar)1339 wcseqm (hstring, wcs, mchar)
1340
1341 char *hstring; /* character string containing FITS header information
1342 in the format <keyword>= <value> [/ <comment>] */
1343 struct WorldCoor *wcs; /* World coordinate system data structure */
1344 char *mchar; /* Suffix character for one of multiple WCS */
1345 {
1346 int ieq = 0;
1347 int eqhead = 0;
1348 char systring[32], eqstring[32];
1349 char radeckey[16], eqkey[16];
1350 char tstring[32];
1351 double ut;
1352
1353 /* Set equinox from EQUINOX, EPOCH, or RADECSYS; default to 2000 */
1354 systring[0] = 0;
1355 eqstring[0] = 0;
1356 if (mchar[0]) {
1357 sprintf (eqkey, "EQUINOX%c", mchar[0]);
1358 sprintf (radeckey, "RADECSYS%c", mchar[0]);
1359 }
1360 else {
1361 strcpy (eqkey, "EQUINOX");
1362 sprintf (radeckey, "RADECSYS");
1363 }
1364 if (!hgets (hstring, eqkey, 31, eqstring)) {
1365 if (hgets (hstring, "EQUINOX", 31, eqstring))
1366 strcpy (eqkey, "EQUINOX");
1367 }
1368 if (!hgets (hstring, radeckey, 31, systring)) {
1369 if (hgets (hstring, "RADECSYS", 31, systring))
1370 sprintf (radeckey, "RADECSYS");
1371 }
1372
1373 if (eqstring[0] == 'J') {
1374 wcs->equinox = atof (eqstring+1);
1375 ieq = atoi (eqstring+1);
1376 strcpy (systring, "FK5");
1377 }
1378 else if (eqstring[0] == 'B') {
1379 wcs->equinox = atof (eqstring+1);
1380 ieq = (int) atof (eqstring+1);
1381 strcpy (systring, "FK4");
1382 }
1383 else if (hgeti4 (hstring, eqkey, &ieq)) {
1384 hgetr8 (hstring, eqkey, &wcs->equinox);
1385 eqhead = 1;
1386 }
1387
1388 else if (hgeti4 (hstring,"EPOCH",&ieq)) {
1389 if (ieq == 0) {
1390 ieq = 1950;
1391 wcs->equinox = 1950.0;
1392 }
1393 else {
1394 hgetr8 (hstring,"EPOCH",&wcs->equinox);
1395 eqhead = 1;
1396 }
1397 }
1398
1399 else if (systring[0] != (char)0) {
1400 if (!strncmp (systring,"FK4",3)) {
1401 wcs->equinox = 1950.0;
1402 ieq = 1950;
1403 }
1404 else if (!strncmp (systring,"ICRS",4)) {
1405 wcs->equinox = 2000.0;
1406 ieq = 2000;
1407 }
1408 else if (!strncmp (systring,"FK5",3)) {
1409 wcs->equinox = 2000.0;
1410 ieq = 2000;
1411 }
1412 else if (!strncmp (systring,"GAL",3)) {
1413 wcs->equinox = 2000.0;
1414 ieq = 2000;
1415 }
1416 else if (!strncmp (systring,"ECL",3)) {
1417 wcs->equinox = 2000.0;
1418 ieq = 2000;
1419 }
1420 }
1421
1422 if (ieq == 0) {
1423 wcs->equinox = 2000.0;
1424 ieq = 2000;
1425 if (!strncmp (wcs->c1type, "RA",2) || !strncmp (wcs->c1type,"DEC",3))
1426 strcpy (systring,"FK5");
1427 }
1428
1429 /* Epoch of image (from observation date, if possible) */
1430 if (!hgetdate (hstring,"DATE-OBS",&wcs->epoch)) {
1431 if (!hgetdate (hstring,"DATE",&wcs->epoch)) {
1432 if (!hgetr8 (hstring,"EPOCH",&wcs->epoch))
1433 wcs->epoch = wcs->equinox;
1434 }
1435 }
1436
1437 /* Add time of day if not part of DATE-OBS string */
1438 else {
1439 hgets (hstring,"DATE-OBS",32,tstring);
1440 if (!strchr (tstring,'T')) {
1441 if (hgetr8 (hstring, "UT",&ut))
1442 wcs->epoch = wcs->epoch + (ut / (24.0 * 365.242198781));
1443 else if (hgetr8 (hstring, "UTMID",&ut))
1444 wcs->epoch = wcs->epoch + (ut / (24.0 * 365.242198781));
1445 }
1446 }
1447 if (wcs->epoch == 0.0)
1448 wcs->epoch = wcs->equinox;
1449
1450 /* Set coordinate system from keyword, if it is present */
1451 if (systring[0] == (char) 0)
1452 hgets (hstring, radeckey, 31, systring);
1453 if (systring[0] != (char) 0) {
1454 strcpy (wcs->radecsys,systring);
1455 if (!eqhead) {
1456 if (!strncmp (wcs->radecsys,"FK4",3))
1457 wcs->equinox = 1950.0;
1458 else if (!strncmp (wcs->radecsys,"FK5",3))
1459 wcs->equinox = 2000.0;
1460 else if (!strncmp (wcs->radecsys,"ICRS",4))
1461 wcs->equinox = 2000.0;
1462 else if (!strncmp (wcs->radecsys,"GAL",3) && ieq == 0)
1463 wcs->equinox = 2000.0;
1464 }
1465 }
1466
1467 /* Otherwise set coordinate system from equinox */
1468 /* Systemless coordinates cannot be translated using b, j, or g commands */
1469 else if (wcs->syswcs != WCS_NPOLE) {
1470 if (ieq > 1980)
1471 strcpy (wcs->radecsys,"FK5");
1472 else
1473 strcpy (wcs->radecsys,"FK4");
1474 }
1475
1476 /* Set galactic coordinates if GLON or GLAT are in C1TYPE */
1477 if (wcs->c1type[0] == 'G')
1478 strcpy (wcs->radecsys,"GALACTIC");
1479 else if (wcs->c1type[0] == 'E')
1480 strcpy (wcs->radecsys,"ECLIPTIC");
1481 else if (wcs->c1type[0] == 'S')
1482 strcpy (wcs->radecsys,"SGALACTC");
1483 else if (wcs->c1type[0] == 'H')
1484 strcpy (wcs->radecsys,"HELIOECL");
1485 else if (wcs->c1type[0] == 'A')
1486 strcpy (wcs->radecsys,"ALTAZ");
1487 else if (wcs->c1type[0] == 'L')
1488 strcpy (wcs->radecsys,"LINEAR");
1489
1490 wcs->syswcs = wcscsys (wcs->radecsys);
1491
1492 return;
1493 }
1494
1495 /* Jun 11 1998 Split off header-dependent WCS initialization from other subs
1496 * Jun 15 1998 Fix major bug in wcsinit() when synthesizing WCS from header
1497 * Jun 18 1998 Fix bug in CD initialization; split PC initialization off
1498 * Jun 18 1998 Split PC initialization off into subroutine wcspcset()
1499 * Jun 24 1998 Set equinox from RADECSYS only if EQUINOX and EPOCH not present
1500 * Jul 6 1998 Read third and fourth axis CTYPEs
1501 * Jul 7 1998 Initialize eqin and eqout to equinox,
1502 * Jul 9 1998 Initialize rotation matrices correctly
1503 * Jul 13 1998 Initialize rotation, scale for polynomial and DSS projections
1504 * Aug 6 1998 Fix CROTA computation for DSS projection
1505 * Sep 4 1998 Fix CROTA, CDELT computation for DSS and polynomial projections
1506 * Sep 14 1998 If DATE-OBS not found, check for DATE
1507 * Sep 14 1998 If B or J present in EQUINOX, use that info to set system
1508 * Sep 29 1998 Initialize additional WCS commands from the environment
1509 * Sep 29 1998 Fix bug which read DATE as number rather than formatted date
1510 * Dec 2 1998 Read projection constants from header (bug fix)
1511 *
1512 * Feb 9 1999 Set rotation angle correctly when using DSS projection
1513 * Feb 19 1999 Fill in CDELTs from scale keyword if absent or zero
1514 * Feb 19 1999 Add PIXSCALE as possible default arcseconds per pixel
1515 * Apr 7 1999 Add error checking for NAXIS and NAXIS1 keywords
1516 * Apr 7 1999 Do not set systring if epoch is 0 and not RA/Dec
1517 * Jul 8 1999 In RADECSYS, use FK5 and FK4 instead of J2000 and B1950
1518 * Oct 15 1999 Free wcs using wcsfree()
1519 * Oct 20 1999 Add multiple WCS support using new subroutine names
1520 * Oct 21 1999 Delete unused variables after lint; declare dsspos()
1521 * Nov 9 1999 Add wcschar() to check WCSNAME keywords for desired WCS
1522 * Nov 9 1999 Check WCSPREx keyword to find out if chained WCS's
1523 *
1524 * Jan 6 1999 Add wcsinitn() to initialize from specific WCSNAME
1525 * Jan 24 2000 Set CD matrix from header even if using polynomial
1526 * Jan 27 2000 Fix MJD to epoch conversion for when MJD-OBS is the only date
1527 * Jan 28 2000 Set CD matrix for DSS projection, too
1528 * Jan 28 2000 Use wcsproj instead of oldwcs
1529 * Dec 18 2000 Fix error in hgets() call in wcschar()
1530 * Dec 29 2000 Compute inverse CD matrix even if polynomial solution
1531 * Dec 29 2000 Add PROJR0 keyword for WCSLIB projections
1532 * Dec 29 2000 Use CDi_j matrix if any elements are present
1533 *
1534 * Jan 31 2001 Fix to allow 1D WCS
1535 * Jan 31 2001 Treat single character WCS name as WCS character
1536 * Feb 20 2001 Implement WCSDEPx nested WCS's
1537 * Feb 23 2001 Initialize all 4 terms of CD matrix
1538 * Feb 28 2001 Fix bug which read CRPIX1 into CRPIX2
1539 * Mar 20 2001 Compare mchar to (char)0, not null
1540 * Mar 21 2001 Move ic declaration into commented out code
1541 * Jul 12 2001 Read PROJPn constants into proj.p array instead of PVn
1542 * Sep 7 2001 Set system to galactic or ecliptic based on CTYPE, not RADECSYS
1543 * Oct 11 2001 Set ctype[0] as well as ctype[1] to TAN for TNX projections
1544 * Oct 19 2001 WCSDIM keyword overrides zero value of NAXIS
1545 *
1546 * Feb 19 2002 Add XPIXSIZE/YPIXSIZE (KPNO) as default image scale keywords
1547 * Mar 12 2002 Add LONPOLE as well as LONGPOLE for WCSLIB 2.8
1548 * Apr 3 2002 Implement hget8c() and hgetsc() to simplify code
1549 * Apr 3 2002 Add PVj_n projection constants in addition to PROJPn
1550 * Apr 19 2002 Increase numeric keyword value length from 16 to 31
1551 * Apr 19 2002 Fix bug which didn't set radecsys keyword name
1552 * Apr 24 2002 If no WCS present for specified letter, return null
1553 * Apr 26 2002 Implement WCSAXESa keyword as first choice for number of axes
1554 * Apr 26 2002 Add wcschar and wcsname to WCS structure
1555 * May 9 2002 Add radvel and zvel to WCS structure
1556 * May 13 2002 Free everything which is allocated
1557 * May 28 2002 Read 10 prj.p instead of maximum of 100
1558 * May 31 2002 Fix bugs with PV reading
1559 * May 31 2002 Initialize syswcs, sysin, sysout in wcsioset()
1560 * Sep 25 2002 Fix subroutine calls for radvel and latpole
1561 * Dec 6 2002 Correctly compute pixel at center of image for default CRPIX
1562 *
1563 * Jan 2 2002 Do not reinitialize projection vector for PV input
1564 * Jan 3 2002 For ZPN, read PVi_0 to PVi_9, not PVi_1 to PVi_10
1565 * Mar 27 2003 Clean up default center computation
1566 * Apr 3 2003 Add input for SIRTF distortion coefficients
1567 * May 8 2003 Change PROJP reading to start with 0 instead of 1
1568 * May 22 2003 Add ZPX approximation, reading projpn from WATi
1569 * May 28 2003 Avoid reinitializing coefficients set by PROJP
1570 * Jun 26 2003 Initialize xref and yref to -999.0
1571 * Sep 23 2003 Change mgets() to mgetstr() to avoid name collision at UCO Lick
1572 * Oct 1 2003 Rename wcs->naxes to wcs->naxis to match WCSLIB 3.2
1573 * Nov 3 2003 Initialize distortion coefficients in distortinit() in distort.c
1574 * Dec 1 2003 Change p[0,1,2] initializations to p[1,2,3]
1575 * Dec 3 2003 Add back wcs->naxes for backward compatibility
1576 * Dec 3 2003 Remove unused variables j,m in wcsinitc()
1577 * Dec 12 2003 Fix call to setwcserr() with format in it
1578 *
1579 * Feb 26 2004 Add parameters for ZPX projection
1580 *
1581 * Jun 22 2005 Drop declaration of variable wcserrmsg which is not used
1582 * Nov 9 2005 Use CROTA1 if CTYPE1 is LAT/DEC, CROTA2 if CTYPE2 is LAT/DEC
1583 *
1584 * Mar 9 2006 Get Epoch of observation from MJD-OBS or DATE-OBS/UT unless DSS
1585 * Apr 24 2006 Initialize rotation matrices
1586 * Apr 25 2006 Ignore axes with dimension of one
1587 * May 19 2006 Initialize all of 9x9 PC matrix; read in loops
1588 * Aug 21 2006 Limit naxes to 2 everywhere; RA and DEC should always be 1st
1589 * Oct 6 2006 If units are pixels, projection type is PIXEL
1590 * Oct 30 2006 Initialize cube face to -1, not a cube projection
1591 *
1592 * Jan 4 2007 Drop declarations of wcsinitc() and wcsinitn() already in wcs.h
1593 * Jan 8 2007 Change WCS letter from char to char*
1594 * Feb 1 2007 Read IRAF log wavelength flag DC-FLAG to wcs.logwcs
1595 * Feb 15 2007 Check for wcs->wcsproj > 0 instead of CTYPEi != LINEAR or PIXEL
1596 * Mar 13 2007 Try for RA, DEC, SECPIX if WCS character is space or null
1597 * Apr 27 2007 Ignore axes with TAB WCS for now
1598 * Oct 17 2007 Fix bug testing &mchar instead of mchar in if statement
1599 *
1600 * May 9 2008 Initialize TNX projection when projection types first set
1601 * Jun 27 2008 If NAXIS1 and NAXIS2 not present, check for IMAGEW and IMAGEH
1602 *
1603 * Mar 24 2009 Fix dimension bug if NAXISi not present (fix from John Burns)
1604 *
1605 * Mar 11 2011 Add NOAO ZPX projection (Frank Valdes)
1606 * Mar 18 2011 Add invert_wcs() by Emmanuel Bertin for SCAMP
1607 * Mar 18 2011 Change Bertin's ARCSEC/DEG to S2D and DEG/ARCSEC to D2S
1608 * Sep 1 2011 Add TPV as TAN with SCAMP PVs
1609 *
1610 * Oct 19 2012 Drop unused variable iszpx; fix bug in latmin assignment
1611 *
1612 * Feb 1 2013 Externalize uppercase()
1613 * Feb 07 2013 Avoid SWARP distortion if SIRTF distortion is present
1614 * Jul 25 2013 Initialize n=0 when checking for SCAMP distortion terms (fix from Martin Kuemmel)
1615 *
1616 * Jun 24 2016 wcs->ptype contains only 3-letter projection code
1617 */
1618