1 #include <stdlib.h>
2 #include <math.h>
3 #include <string.h>
4 #include "fitsio2.h"
5
6 /*--------------------------------------------------------------------------*/
fits_read_wcstab(fitsfile * fptr,int nwtb,wtbarr * wtb,int * status)7 int fits_read_wcstab(
8 fitsfile *fptr, /* I - FITS file pointer */
9 int nwtb, /* Number of arrays to be read from the binary table(s) */
10 wtbarr *wtb, /* Address of the first element of an array of wtbarr
11 typedefs. This wtbarr typedef is defined below to
12 match the wtbarr struct defined in WCSLIB. An array
13 of such structs returned by the WCSLIB function
14 wcstab(). */
15 int *status)
16
17 /*
18 * Author: Mark Calabretta, Australia Telescope National Facility
19 * http://www.atnf.csiro.au/~mcalabre/index.html
20 *
21 * fits_read_wcstab() extracts arrays from a binary table required in
22 * constructing -TAB coordinates. This helper routine is intended for
23 * use by routines in the WCSLIB library when dealing with the -TAB table
24 * look up WCS convention.
25 */
26
27 {
28 int anynul, colnum, hdunum, iwtb, m, naxis, nostat;
29 long *naxes = 0, nelem;
30 wtbarr *wtbp;
31
32
33 if (*status) return *status;
34
35 if (fptr == 0) {
36 return (*status = NULL_INPUT_PTR);
37 }
38
39 if (nwtb == 0) return 0;
40
41 /* Zero the array pointers. */
42 wtbp = wtb;
43 for (iwtb = 0; iwtb < nwtb; iwtb++, wtbp++) {
44 *wtbp->arrayp = 0x0;
45 }
46
47 /* Save HDU number so that we can move back to it later. */
48 fits_get_hdu_num(fptr, &hdunum);
49
50 wtbp = wtb;
51 for (iwtb = 0; iwtb < nwtb; iwtb++, wtbp++) {
52 /* Move to the required binary table extension. */
53 if (fits_movnam_hdu(fptr, BINARY_TBL, (char *)(wtbp->extnam),
54 wtbp->extver, status)) {
55 goto cleanup;
56 }
57
58 /* Locate the table column. */
59 if (fits_get_colnum(fptr, CASEINSEN, (char *)(wtbp->ttype), &colnum,
60 status)) {
61 goto cleanup;
62 }
63
64 /* Get the array dimensions and check for consistency. */
65 if (wtbp->ndim < 1) {
66 *status = NEG_AXIS;
67 goto cleanup;
68 }
69
70 if (!(naxes = calloc(wtbp->ndim, sizeof(long)))) {
71 *status = MEMORY_ALLOCATION;
72 goto cleanup;
73 }
74
75 if (fits_read_tdim(fptr, colnum, wtbp->ndim, &naxis, naxes, status)) {
76 goto cleanup;
77 }
78
79 if (naxis != wtbp->ndim) {
80 if (wtbp->kind == 'c' && wtbp->ndim == 2) {
81 /* Allow TDIMn to be omitted for degenerate coordinate arrays. */
82 naxis = 2;
83 naxes[1] = naxes[0];
84 naxes[0] = 1;
85 } else {
86 *status = BAD_TDIM;
87 goto cleanup;
88 }
89 }
90
91 if (wtbp->kind == 'c') {
92 /* Coordinate array; calculate the array size. */
93 nelem = naxes[0];
94 for (m = 0; m < naxis-1; m++) {
95 *(wtbp->dimlen + m) = naxes[m+1];
96 nelem *= naxes[m+1];
97 }
98 } else {
99 /* Index vector; check length. */
100 if ((nelem = naxes[0]) != *(wtbp->dimlen)) {
101 /* N.B. coordinate array precedes the index vectors. */
102 *status = BAD_TDIM;
103 goto cleanup;
104 }
105 }
106
107 free(naxes);
108 naxes = 0;
109
110 /* Allocate memory for the array. */
111 if (!(*wtbp->arrayp = calloc((size_t)nelem, sizeof(double)))) {
112 *status = MEMORY_ALLOCATION;
113 goto cleanup;
114 }
115
116 /* Read the array from the table. */
117 if (fits_read_col_dbl(fptr, colnum, wtbp->row, 1L, nelem, 0.0,
118 *wtbp->arrayp, &anynul, status)) {
119 goto cleanup;
120 }
121 }
122
123 cleanup:
124 /* Move back to the starting HDU. */
125 nostat = 0;
126 fits_movabs_hdu(fptr, hdunum, 0, &nostat);
127
128 /* Release allocated memory. */
129 if (naxes) free(naxes);
130 if (*status) {
131 wtbp = wtb;
132 for (iwtb = 0; iwtb < nwtb; iwtb++, wtbp++) {
133 if (*wtbp->arrayp) free(*wtbp->arrayp);
134 }
135 }
136
137 return *status;
138 }
139 /*--------------------------------------------------------------------------*/
ffgiwcs(fitsfile * fptr,char ** header,int * status)140 int ffgiwcs(fitsfile *fptr, /* I - FITS file pointer */
141 char **header, /* O - pointer to the WCS related keywords */
142 int *status) /* IO - error status */
143 /*
144 int fits_get_image_wcs_keys
145 return a string containing all the image WCS header keywords.
146 This string is then used as input to the wcsinit WCSlib routine.
147
148 THIS ROUTINE IS DEPRECATED. USE fits_hdr2str INSTEAD
149 */
150 {
151 int hdutype;
152
153 if (*status > 0)
154 return(*status);
155
156 fits_get_hdu_type(fptr, &hdutype, status);
157 if (hdutype != IMAGE_HDU)
158 {
159 ffpmsg(
160 "Error in ffgiwcs. This HDU is not an image. Can't read WCS keywords");
161 return(*status = NOT_IMAGE);
162 }
163
164 /* read header keywords into a long string of chars */
165 if (ffh2st(fptr, header, status) > 0)
166 {
167 ffpmsg("error creating string of image WCS keywords (ffgiwcs)");
168 return(*status);
169 }
170
171 return(*status);
172 }
173
174 /*--------------------------------------------------------------------------*/
ffgics(fitsfile * fptr,double * xrval,double * yrval,double * xrpix,double * yrpix,double * xinc,double * yinc,double * rot,char * type,int * status)175 int ffgics(fitsfile *fptr, /* I - FITS file pointer */
176 double *xrval, /* O - X reference value */
177 double *yrval, /* O - Y reference value */
178 double *xrpix, /* O - X reference pixel */
179 double *yrpix, /* O - Y reference pixel */
180 double *xinc, /* O - X increment per pixel */
181 double *yinc, /* O - Y increment per pixel */
182 double *rot, /* O - rotation angle (degrees) */
183 char *type, /* O - type of projection ('-tan') */
184 int *status) /* IO - error status */
185 /*
186 read the values of the celestial coordinate system keywords.
187 These values may be used as input to the subroutines that
188 calculate celestial coordinates. (ffxypx, ffwldp)
189
190 Modified in Nov 1999 to convert the CD matrix keywords back
191 to the old CDELTn form, and to swap the axes if the dec-like
192 axis is given first, and to assume default values if any of the
193 keywords are not present.
194 */
195 {
196 int tstat = 0, cd_exists = 0, pc_exists = 0;
197 char ctype[FLEN_VALUE];
198 double cd11 = 0.0, cd21 = 0.0, cd22 = 0.0, cd12 = 0.0;
199 double pc11 = 1.0, pc21 = 0.0, pc22 = 1.0, pc12 = 0.0;
200 double pi = 3.1415926535897932;
201 double phia, phib, temp;
202 double toler = .0002; /* tolerance for angles to agree (radians) */
203 /* (= approximately 0.01 degrees) */
204
205 if (*status > 0)
206 return(*status);
207
208 tstat = 0;
209 if (ffgkyd(fptr, "CRVAL1", xrval, NULL, &tstat))
210 *xrval = 0.;
211
212 tstat = 0;
213 if (ffgkyd(fptr, "CRVAL2", yrval, NULL, &tstat))
214 *yrval = 0.;
215
216 tstat = 0;
217 if (ffgkyd(fptr, "CRPIX1", xrpix, NULL, &tstat))
218 *xrpix = 0.;
219
220 tstat = 0;
221 if (ffgkyd(fptr, "CRPIX2", yrpix, NULL, &tstat))
222 *yrpix = 0.;
223
224 /* look for CDELTn first, then CDi_j keywords */
225 tstat = 0;
226 if (ffgkyd(fptr, "CDELT1", xinc, NULL, &tstat))
227 {
228 /* CASE 1: no CDELTn keyword, so look for the CD matrix */
229 tstat = 0;
230 if (ffgkyd(fptr, "CD1_1", &cd11, NULL, &tstat))
231 tstat = 0; /* reset keyword not found error */
232 else
233 cd_exists = 1; /* found at least 1 CD_ keyword */
234
235 if (ffgkyd(fptr, "CD2_1", &cd21, NULL, &tstat))
236 tstat = 0; /* reset keyword not found error */
237 else
238 cd_exists = 1; /* found at least 1 CD_ keyword */
239
240 if (ffgkyd(fptr, "CD1_2", &cd12, NULL, &tstat))
241 tstat = 0; /* reset keyword not found error */
242 else
243 cd_exists = 1; /* found at least 1 CD_ keyword */
244
245 if (ffgkyd(fptr, "CD2_2", &cd22, NULL, &tstat))
246 tstat = 0; /* reset keyword not found error */
247 else
248 cd_exists = 1; /* found at least 1 CD_ keyword */
249
250 if (cd_exists) /* convert CDi_j back to CDELTn */
251 {
252 /* there are 2 ways to compute the angle: */
253 phia = atan2( cd21, cd11);
254 phib = atan2(-cd12, cd22);
255
256 /* ensure that phia <= phib */
257 temp = minvalue(phia, phib);
258 phib = maxvalue(phia, phib);
259 phia = temp;
260
261 /* there is a possible 180 degree ambiguity in the angles */
262 /* so add 180 degress to the smaller value if the values */
263 /* differ by more than 90 degrees = pi/2 radians. */
264 /* (Later, we may decide to take the other solution by */
265 /* subtracting 180 degrees from the larger value). */
266
267 if ((phib - phia) > (pi / 2.))
268 phia += pi;
269
270 if (fabs(phia - phib) > toler)
271 {
272 /* angles don't agree, so looks like there is some skewness */
273 /* between the axes. Return with an error to be safe. */
274 *status = APPROX_WCS_KEY;
275 }
276
277 phia = (phia + phib) /2.; /* use the average of the 2 values */
278 *xinc = cd11 / cos(phia);
279 *yinc = cd22 / cos(phia);
280 *rot = phia * 180. / pi;
281
282 /* common usage is to have a positive yinc value. If it is */
283 /* negative, then subtract 180 degrees from rot and negate */
284 /* both xinc and yinc. */
285
286 if (*yinc < 0)
287 {
288 *xinc = -(*xinc);
289 *yinc = -(*yinc);
290 *rot = *rot - 180.;
291 }
292 }
293 else /* no CD matrix keywords either */
294 {
295 *xinc = 1.;
296
297 /* there was no CDELT1 keyword, but check for CDELT2 just in case */
298 tstat = 0;
299 if (ffgkyd(fptr, "CDELT2", yinc, NULL, &tstat))
300 *yinc = 1.;
301
302 tstat = 0;
303 if (ffgkyd(fptr, "CROTA2", rot, NULL, &tstat))
304 *rot=0.;
305 }
306 }
307 else /* Case 2: CDELTn + optional PC matrix */
308 {
309 if (ffgkyd(fptr, "CDELT2", yinc, NULL, &tstat))
310 *yinc = 1.;
311
312 tstat = 0;
313 if (ffgkyd(fptr, "CROTA2", rot, NULL, &tstat))
314 {
315 *rot=0.;
316
317 /* no CROTA2 keyword, so look for the PC matrix */
318 tstat = 0;
319 if (ffgkyd(fptr, "PC1_1", &pc11, NULL, &tstat))
320 tstat = 0; /* reset keyword not found error */
321 else
322 pc_exists = 1; /* found at least 1 PC_ keyword */
323
324 if (ffgkyd(fptr, "PC2_1", &pc21, NULL, &tstat))
325 tstat = 0; /* reset keyword not found error */
326 else
327 pc_exists = 1; /* found at least 1 PC_ keyword */
328
329 if (ffgkyd(fptr, "PC1_2", &pc12, NULL, &tstat))
330 tstat = 0; /* reset keyword not found error */
331 else
332 pc_exists = 1; /* found at least 1 PC_ keyword */
333
334 if (ffgkyd(fptr, "PC2_2", &pc22, NULL, &tstat))
335 tstat = 0; /* reset keyword not found error */
336 else
337 pc_exists = 1; /* found at least 1 PC_ keyword */
338
339 if (pc_exists) /* convert PCi_j back to CDELTn */
340 {
341 /* there are 2 ways to compute the angle: */
342 phia = atan2( pc21, pc11);
343 phib = atan2(-pc12, pc22);
344
345 /* ensure that phia <= phib */
346 temp = minvalue(phia, phib);
347 phib = maxvalue(phia, phib);
348 phia = temp;
349
350 /* there is a possible 180 degree ambiguity in the angles */
351 /* so add 180 degress to the smaller value if the values */
352 /* differ by more than 90 degrees = pi/2 radians. */
353 /* (Later, we may decide to take the other solution by */
354 /* subtracting 180 degrees from the larger value). */
355
356 if ((phib - phia) > (pi / 2.))
357 phia += pi;
358
359 if (fabs(phia - phib) > toler)
360 {
361 /* angles don't agree, so looks like there is some skewness */
362 /* between the axes. Return with an error to be safe. */
363 *status = APPROX_WCS_KEY;
364 }
365
366 phia = (phia + phib) /2.; /* use the average of the 2 values */
367 *rot = phia * 180. / pi;
368 }
369 }
370 }
371
372 /* get the type of projection, if any */
373 tstat = 0;
374 if (ffgkys(fptr, "CTYPE1", ctype, NULL, &tstat))
375 type[0] = '\0';
376 else
377 {
378 /* copy the projection type string */
379 strncpy(type, &ctype[4], 4);
380 type[4] = '\0';
381
382 /* check if RA and DEC are inverted */
383 if (!strncmp(ctype, "DEC-", 4) || !strncmp(ctype+1, "LAT", 3))
384 {
385 /* the latitudinal axis is given first, so swap them */
386
387 /*
388 this case was removed on 12/9. Apparently not correct.
389
390 if ((*xinc / *yinc) < 0. )
391 *rot = -90. - (*rot);
392 else
393 */
394 *rot = 90. - (*rot);
395
396 /* Empirical tests with ds9 show the y-axis sign must be negated */
397 /* and the xinc and yinc values must NOT be swapped. */
398 *yinc = -(*yinc);
399
400 temp = *xrval;
401 *xrval = *yrval;
402 *yrval = temp;
403 }
404 }
405
406 return(*status);
407 }
408 /*--------------------------------------------------------------------------*/
ffgicsa(fitsfile * fptr,char version,double * xrval,double * yrval,double * xrpix,double * yrpix,double * xinc,double * yinc,double * rot,char * type,int * status)409 int ffgicsa(fitsfile *fptr, /* I - FITS file pointer */
410 char version, /* I - character code of desired version */
411 /* A - Z or blank */
412 double *xrval, /* O - X reference value */
413 double *yrval, /* O - Y reference value */
414 double *xrpix, /* O - X reference pixel */
415 double *yrpix, /* O - Y reference pixel */
416 double *xinc, /* O - X increment per pixel */
417 double *yinc, /* O - Y increment per pixel */
418 double *rot, /* O - rotation angle (degrees) */
419 char *type, /* O - type of projection ('-tan') */
420 int *status) /* IO - error status */
421 /*
422 read the values of the celestial coordinate system keywords.
423 These values may be used as input to the subroutines that
424 calculate celestial coordinates. (ffxypx, ffwldp)
425
426 Modified in Nov 1999 to convert the CD matrix keywords back
427 to the old CDELTn form, and to swap the axes if the dec-like
428 axis is given first, and to assume default values if any of the
429 keywords are not present.
430 */
431 {
432 int tstat = 0, cd_exists = 0, pc_exists = 0;
433 char ctype[FLEN_VALUE], keyname[FLEN_VALUE], alt[2];
434 double cd11 = 0.0, cd21 = 0.0, cd22 = 0.0, cd12 = 0.0;
435 double pc11 = 1.0, pc21 = 0.0, pc22 = 1.0, pc12 = 0.0;
436 double pi = 3.1415926535897932;
437 double phia, phib, temp;
438 double toler = .0002; /* tolerance for angles to agree (radians) */
439 /* (= approximately 0.01 degrees) */
440
441 if (*status > 0)
442 return(*status);
443
444 if (version == ' ') {
445 ffgics(fptr, xrval, yrval, xrpix, yrpix, xinc, yinc, rot, type, status);
446 return (*status);
447 }
448
449 if (version > 'Z' || version < 'A') {
450 ffpmsg("ffgicsa: illegal WCS version code (must be A - Z or blank)");
451 return(*status = WCS_ERROR);
452 }
453
454 alt[0] = version;
455 alt[1] = '\0';
456
457 tstat = 0;
458 strcpy(keyname, "CRVAL1");
459 strcat(keyname, alt);
460 if (ffgkyd(fptr, keyname, xrval, NULL, &tstat))
461 *xrval = 0.;
462
463 tstat = 0;
464 strcpy(keyname, "CRVAL2");
465 strcat(keyname, alt);
466 if (ffgkyd(fptr, keyname, yrval, NULL, &tstat))
467 *yrval = 0.;
468
469 tstat = 0;
470 strcpy(keyname, "CRPIX1");
471 strcat(keyname, alt);
472 if (ffgkyd(fptr, keyname, xrpix, NULL, &tstat))
473 *xrpix = 0.;
474
475 tstat = 0;
476 strcpy(keyname, "CRPIX2");
477 strcat(keyname, alt);
478 if (ffgkyd(fptr, keyname, yrpix, NULL, &tstat))
479 *yrpix = 0.;
480
481 /* look for CDELTn first, then CDi_j keywords */
482 tstat = 0;
483 strcpy(keyname, "CDELT1");
484 strcat(keyname, alt);
485 if (ffgkyd(fptr, keyname, xinc, NULL, &tstat))
486 {
487 /* CASE 1: no CDELTn keyword, so look for the CD matrix */
488 tstat = 0;
489 strcpy(keyname, "CD1_1");
490 strcat(keyname, alt);
491 if (ffgkyd(fptr, keyname, &cd11, NULL, &tstat))
492 tstat = 0; /* reset keyword not found error */
493 else
494 cd_exists = 1; /* found at least 1 CD_ keyword */
495
496 strcpy(keyname, "CD2_1");
497 strcat(keyname, alt);
498 if (ffgkyd(fptr, keyname, &cd21, NULL, &tstat))
499 tstat = 0; /* reset keyword not found error */
500 else
501 cd_exists = 1; /* found at least 1 CD_ keyword */
502
503 strcpy(keyname, "CD1_2");
504 strcat(keyname, alt);
505 if (ffgkyd(fptr, keyname, &cd12, NULL, &tstat))
506 tstat = 0; /* reset keyword not found error */
507 else
508 cd_exists = 1; /* found at least 1 CD_ keyword */
509
510 strcpy(keyname, "CD2_2");
511 strcat(keyname, alt);
512 if (ffgkyd(fptr, keyname, &cd22, NULL, &tstat))
513 tstat = 0; /* reset keyword not found error */
514 else
515 cd_exists = 1; /* found at least 1 CD_ keyword */
516
517 if (cd_exists) /* convert CDi_j back to CDELTn */
518 {
519 /* there are 2 ways to compute the angle: */
520 phia = atan2( cd21, cd11);
521 phib = atan2(-cd12, cd22);
522
523 /* ensure that phia <= phib */
524 temp = minvalue(phia, phib);
525 phib = maxvalue(phia, phib);
526 phia = temp;
527
528 /* there is a possible 180 degree ambiguity in the angles */
529 /* so add 180 degress to the smaller value if the values */
530 /* differ by more than 90 degrees = pi/2 radians. */
531 /* (Later, we may decide to take the other solution by */
532 /* subtracting 180 degrees from the larger value). */
533
534 if ((phib - phia) > (pi / 2.))
535 phia += pi;
536
537 if (fabs(phia - phib) > toler)
538 {
539 /* angles don't agree, so looks like there is some skewness */
540 /* between the axes. Return with an error to be safe. */
541 *status = APPROX_WCS_KEY;
542 }
543
544 phia = (phia + phib) /2.; /* use the average of the 2 values */
545 *xinc = cd11 / cos(phia);
546 *yinc = cd22 / cos(phia);
547 *rot = phia * 180. / pi;
548
549 /* common usage is to have a positive yinc value. If it is */
550 /* negative, then subtract 180 degrees from rot and negate */
551 /* both xinc and yinc. */
552
553 if (*yinc < 0)
554 {
555 *xinc = -(*xinc);
556 *yinc = -(*yinc);
557 *rot = *rot - 180.;
558 }
559 }
560 else /* no CD matrix keywords either */
561 {
562 *xinc = 1.;
563
564 /* there was no CDELT1 keyword, but check for CDELT2 just in case */
565 tstat = 0;
566 strcpy(keyname, "CDELT2");
567 strcat(keyname, alt);
568 if (ffgkyd(fptr, keyname, yinc, NULL, &tstat))
569 *yinc = 1.;
570
571 tstat = 0;
572 strcpy(keyname, "CROTA2");
573 strcat(keyname, alt);
574 if (ffgkyd(fptr, keyname, rot, NULL, &tstat))
575 *rot=0.;
576 }
577 }
578 else /* Case 2: CDELTn + optional PC matrix */
579 {
580 strcpy(keyname, "CDELT2");
581 strcat(keyname, alt);
582 if (ffgkyd(fptr, keyname, yinc, NULL, &tstat))
583 *yinc = 1.;
584
585 tstat = 0;
586 strcpy(keyname, "CROTA2");
587 strcat(keyname, alt);
588 if (ffgkyd(fptr, keyname, rot, NULL, &tstat))
589 {
590 *rot=0.;
591
592 /* no CROTA2 keyword, so look for the PC matrix */
593 tstat = 0;
594 strcpy(keyname, "PC1_1");
595 strcat(keyname, alt);
596 if (ffgkyd(fptr, keyname, &pc11, NULL, &tstat))
597 tstat = 0; /* reset keyword not found error */
598 else
599 pc_exists = 1; /* found at least 1 PC_ keyword */
600
601 strcpy(keyname, "PC2_1");
602 strcat(keyname, alt);
603 if (ffgkyd(fptr, keyname, &pc21, NULL, &tstat))
604 tstat = 0; /* reset keyword not found error */
605 else
606 pc_exists = 1; /* found at least 1 PC_ keyword */
607
608 strcpy(keyname, "PC1_2");
609 strcat(keyname, alt);
610 if (ffgkyd(fptr, keyname, &pc12, NULL, &tstat))
611 tstat = 0; /* reset keyword not found error */
612 else
613 pc_exists = 1; /* found at least 1 PC_ keyword */
614
615 strcpy(keyname, "PC2_2");
616 strcat(keyname, alt);
617 if (ffgkyd(fptr, keyname, &pc22, NULL, &tstat))
618 tstat = 0; /* reset keyword not found error */
619 else
620 pc_exists = 1; /* found at least 1 PC_ keyword */
621
622 if (pc_exists) /* convert PCi_j back to CDELTn */
623 {
624 /* there are 2 ways to compute the angle: */
625 phia = atan2( pc21, pc11);
626 phib = atan2(-pc12, pc22);
627
628 /* ensure that phia <= phib */
629 temp = minvalue(phia, phib);
630 phib = maxvalue(phia, phib);
631 phia = temp;
632
633 /* there is a possible 180 degree ambiguity in the angles */
634 /* so add 180 degress to the smaller value if the values */
635 /* differ by more than 90 degrees = pi/2 radians. */
636 /* (Later, we may decide to take the other solution by */
637 /* subtracting 180 degrees from the larger value). */
638
639 if ((phib - phia) > (pi / 2.))
640 phia += pi;
641
642 if (fabs(phia - phib) > toler)
643 {
644 /* angles don't agree, so looks like there is some skewness */
645 /* between the axes. Return with an error to be safe. */
646 *status = APPROX_WCS_KEY;
647 }
648
649 phia = (phia + phib) /2.; /* use the average of the 2 values */
650 *rot = phia * 180. / pi;
651 }
652 }
653 }
654
655 /* get the type of projection, if any */
656 tstat = 0;
657 strcpy(keyname, "CTYPE1");
658 strcat(keyname, alt);
659 if (ffgkys(fptr, keyname, ctype, NULL, &tstat))
660 type[0] = '\0';
661 else
662 {
663 /* copy the projection type string */
664 strncpy(type, &ctype[4], 4);
665 type[4] = '\0';
666
667 /* check if RA and DEC are inverted */
668 if (!strncmp(ctype, "DEC-", 4) || !strncmp(ctype+1, "LAT", 3))
669 {
670 /* the latitudinal axis is given first, so swap them */
671
672 *rot = 90. - (*rot);
673
674 /* Empirical tests with ds9 show the y-axis sign must be negated */
675 /* and the xinc and yinc values must NOT be swapped. */
676 *yinc = -(*yinc);
677
678 temp = *xrval;
679 *xrval = *yrval;
680 *yrval = temp;
681 }
682 }
683
684 return(*status);
685 }
686 /*--------------------------------------------------------------------------*/
ffgtcs(fitsfile * fptr,int xcol,int ycol,double * xrval,double * yrval,double * xrpix,double * yrpix,double * xinc,double * yinc,double * rot,char * type,int * status)687 int ffgtcs(fitsfile *fptr, /* I - FITS file pointer */
688 int xcol, /* I - column containing the RA coordinate */
689 int ycol, /* I - column containing the DEC coordinate */
690 double *xrval, /* O - X reference value */
691 double *yrval, /* O - Y reference value */
692 double *xrpix, /* O - X reference pixel */
693 double *yrpix, /* O - Y reference pixel */
694 double *xinc, /* O - X increment per pixel */
695 double *yinc, /* O - Y increment per pixel */
696 double *rot, /* O - rotation angle (degrees) */
697 char *type, /* O - type of projection ('-sin') */
698 int *status) /* IO - error status */
699 /*
700 read the values of the celestial coordinate system keywords
701 from a FITS table where the X and Y or RA and DEC coordinates
702 are stored in separate column. Do this by converting the
703 table to a temporary FITS image, then reading the keywords
704 from the image file.
705 These values may be used as input to the subroutines that
706 calculate celestial coordinates. (ffxypx, ffwldp)
707 */
708 {
709 int colnum[2];
710 long naxes[2];
711 fitsfile *tptr;
712
713 if (*status > 0)
714 return(*status);
715
716 colnum[0] = xcol;
717 colnum[1] = ycol;
718 naxes[0] = 10;
719 naxes[1] = 10;
720
721 /* create temporary FITS file, in memory */
722 ffinit(&tptr, "mem://", status);
723
724 /* create a temporary image; the datatype and size are not important */
725 ffcrim(tptr, 32, 2, naxes, status);
726
727 /* now copy the relevant keywords from the table to the image */
728 fits_copy_pixlist2image(fptr, tptr, 9, 2, colnum, status);
729
730 /* write default WCS keywords, if they are not present */
731 fits_write_keys_histo(fptr, tptr, 2, colnum, status);
732
733 if (*status > 0)
734 return(*status);
735
736 /* read the WCS keyword values from the temporary image */
737 ffgics(tptr, xrval, yrval, xrpix, yrpix, xinc, yinc, rot, type, status);
738
739 if (*status > 0)
740 {
741 ffpmsg
742 ("ffgtcs could not find all the celestial coordinate keywords");
743 return(*status = NO_WCS_KEY);
744 }
745
746 /* delete the temporary file */
747 fits_delete_file(tptr, status);
748
749 return(*status);
750 }
751 /*--------------------------------------------------------------------------*/
ffgtwcs(fitsfile * fptr,int xcol,int ycol,char ** header,int * status)752 int ffgtwcs(fitsfile *fptr, /* I - FITS file pointer */
753 int xcol, /* I - column number for the X column */
754 int ycol, /* I - column number for the Y column */
755 char **header, /* O - string of all the WCS keywords */
756 int *status) /* IO - error status */
757 /*
758 int fits_get_table_wcs_keys
759 Return string containing all the WCS keywords appropriate for the
760 pair of X and Y columns containing the coordinate
761 of each event in an event list table. This string may then be passed
762 to Doug Mink's WCS library wcsinit routine, to create and initialize the
763 WCS structure. The calling routine must free the header character string
764 when it is no longer needed.
765
766 THIS ROUTINE IS DEPRECATED. USE fits_hdr2str INSTEAD
767 */
768 {
769 int hdutype, ncols, tstatus, length;
770 int naxis1 = 1, naxis2 = 1;
771 long tlmin, tlmax;
772 char keyname[FLEN_KEYWORD];
773 char valstring[FLEN_VALUE];
774 char comm[2];
775 char *cptr;
776 /* construct a string of 80 blanks, for adding fill to the keywords */
777 /* 12345678901234567890123456789012345678901234567890123456789012345678901234567890 */
778 char blanks[] = " ";
779
780 if (*status > 0)
781 return(*status);
782
783 fits_get_hdu_type(fptr, &hdutype, status);
784 if (hdutype == IMAGE_HDU)
785 {
786 ffpmsg("Can't read table WSC keywords. This HDU is not a table");
787 return(*status = NOT_TABLE);
788 }
789
790 fits_get_num_cols(fptr, &ncols, status);
791
792 if (xcol < 1 || xcol > ncols)
793 {
794 ffpmsg("illegal X axis column number in fftwcs");
795 return(*status = BAD_COL_NUM);
796 }
797
798 if (ycol < 1 || ycol > ncols)
799 {
800 ffpmsg("illegal Y axis column number in fftwcs");
801 return(*status = BAD_COL_NUM);
802 }
803
804 /* allocate character string for all the WCS keywords */
805 *header = calloc(1, 2401); /* room for up to 30 keywords */
806 if (*header == 0)
807 {
808 ffpmsg("error allocating memory for WCS header keywords (fftwcs)");
809 return(*status = MEMORY_ALLOCATION);
810 }
811
812 cptr = *header;
813 comm[0] = '\0';
814
815 tstatus = 0;
816 ffkeyn("TLMIN",xcol,keyname,status);
817 ffgkyj(fptr,keyname, &tlmin,NULL,&tstatus);
818
819 if (!tstatus)
820 {
821 ffkeyn("TLMAX",xcol,keyname,status);
822 ffgkyj(fptr,keyname, &tlmax,NULL,&tstatus);
823 }
824
825 if (!tstatus)
826 {
827 naxis1 = tlmax - tlmin + 1;
828 }
829
830 tstatus = 0;
831 ffkeyn("TLMIN",ycol,keyname,status);
832 ffgkyj(fptr,keyname, &tlmin,NULL,&tstatus);
833
834 if (!tstatus)
835 {
836 ffkeyn("TLMAX",ycol,keyname,status);
837 ffgkyj(fptr,keyname, &tlmax,NULL,&tstatus);
838 }
839
840 if (!tstatus)
841 {
842 naxis2 = tlmax - tlmin + 1;
843 }
844
845 /* 123456789012345678901234567890 */
846 strcat(cptr, "NAXIS = 2");
847 strncat(cptr, blanks, 50);
848 cptr += 80;
849
850 ffi2c(naxis1, valstring, status); /* convert to formatted string */
851 ffmkky("NAXIS1", valstring, comm, cptr, status); /* construct the keyword*/
852 strncat(cptr, blanks, 50); /* pad with blanks */
853 cptr += 80;
854
855 strcpy(keyname, "NAXIS2");
856 ffi2c(naxis2, valstring, status); /* convert to formatted string */
857 ffmkky(keyname, valstring, comm, cptr, status); /* construct the keyword*/
858 strncat(cptr, blanks, 50); /* pad with blanks */
859 cptr += 80;
860
861 /* read the required header keywords (use defaults if not found) */
862
863 /* CTYPE1 keyword */
864 tstatus = 0;
865 ffkeyn("TCTYP",xcol,keyname,status);
866 if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
867 valstring[0] = '\0';
868 ffmkky("CTYPE1", valstring, comm, cptr, status); /* construct the keyword*/
869 length = strlen(cptr);
870 strncat(cptr, blanks, 80 - length); /* pad with blanks */
871 cptr += 80;
872
873 /* CTYPE2 keyword */
874 tstatus = 0;
875 ffkeyn("TCTYP",ycol,keyname,status);
876 if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
877 valstring[0] = '\0';
878 ffmkky("CTYPE2", valstring, comm, cptr, status); /* construct the keyword*/
879 length = strlen(cptr);
880 strncat(cptr, blanks, 80 - length); /* pad with blanks */
881 cptr += 80;
882
883 /* CRPIX1 keyword */
884 tstatus = 0;
885 ffkeyn("TCRPX",xcol,keyname,status);
886 if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
887 strcpy(valstring, "1");
888 ffmkky("CRPIX1", valstring, comm, cptr, status); /* construct the keyword*/
889 strncat(cptr, blanks, 50); /* pad with blanks */
890 cptr += 80;
891
892 /* CRPIX2 keyword */
893 tstatus = 0;
894 ffkeyn("TCRPX",ycol,keyname,status);
895 if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
896 strcpy(valstring, "1");
897 ffmkky("CRPIX2", valstring, comm, cptr, status); /* construct the keyword*/
898 strncat(cptr, blanks, 50); /* pad with blanks */
899 cptr += 80;
900
901 /* CRVAL1 keyword */
902 tstatus = 0;
903 ffkeyn("TCRVL",xcol,keyname,status);
904 if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
905 strcpy(valstring, "1");
906 ffmkky("CRVAL1", valstring, comm, cptr, status); /* construct the keyword*/
907 strncat(cptr, blanks, 50); /* pad with blanks */
908 cptr += 80;
909
910 /* CRVAL2 keyword */
911 tstatus = 0;
912 ffkeyn("TCRVL",ycol,keyname,status);
913 if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
914 strcpy(valstring, "1");
915 ffmkky("CRVAL2", valstring, comm, cptr, status); /* construct the keyword*/
916 strncat(cptr, blanks, 50); /* pad with blanks */
917 cptr += 80;
918
919 /* CDELT1 keyword */
920 tstatus = 0;
921 ffkeyn("TCDLT",xcol,keyname,status);
922 if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
923 strcpy(valstring, "1");
924 ffmkky("CDELT1", valstring, comm, cptr, status); /* construct the keyword*/
925 strncat(cptr, blanks, 50); /* pad with blanks */
926 cptr += 80;
927
928 /* CDELT2 keyword */
929 tstatus = 0;
930 ffkeyn("TCDLT",ycol,keyname,status);
931 if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) )
932 strcpy(valstring, "1");
933 ffmkky("CDELT2", valstring, comm, cptr, status); /* construct the keyword*/
934 strncat(cptr, blanks, 50); /* pad with blanks */
935 cptr += 80;
936
937 /* the following keywords may not exist */
938
939 /* CROTA2 keyword */
940 tstatus = 0;
941 ffkeyn("TCROT",ycol,keyname,status);
942 if (ffgkey(fptr, keyname, valstring, NULL, &tstatus) == 0 )
943 {
944 ffmkky("CROTA2", valstring, comm, cptr, status); /* construct keyword*/
945 strncat(cptr, blanks, 50); /* pad with blanks */
946 cptr += 80;
947 }
948
949 /* EPOCH keyword */
950 tstatus = 0;
951 if (ffgkey(fptr, "EPOCH", valstring, NULL, &tstatus) == 0 )
952 {
953 ffmkky("EPOCH", valstring, comm, cptr, status); /* construct keyword*/
954 length = strlen(cptr);
955 strncat(cptr, blanks, 80 - length); /* pad with blanks */
956 cptr += 80;
957 }
958
959 /* EQUINOX keyword */
960 tstatus = 0;
961 if (ffgkey(fptr, "EQUINOX", valstring, NULL, &tstatus) == 0 )
962 {
963 ffmkky("EQUINOX", valstring, comm, cptr, status); /* construct keyword*/
964 length = strlen(cptr);
965 strncat(cptr, blanks, 80 - length); /* pad with blanks */
966 cptr += 80;
967 }
968
969 /* RADECSYS keyword */
970 tstatus = 0;
971 if (ffgkey(fptr, "RADECSYS", valstring, NULL, &tstatus) == 0 )
972 {
973 ffmkky("RADECSYS", valstring, comm, cptr, status); /*construct keyword*/
974 length = strlen(cptr);
975 strncat(cptr, blanks, 80 - length); /* pad with blanks */
976 cptr += 80;
977 }
978
979 /* TELESCOPE keyword */
980 tstatus = 0;
981 if (ffgkey(fptr, "TELESCOP", valstring, NULL, &tstatus) == 0 )
982 {
983 ffmkky("TELESCOP", valstring, comm, cptr, status);
984 length = strlen(cptr);
985 strncat(cptr, blanks, 80 - length); /* pad with blanks */
986 cptr += 80;
987 }
988
989 /* INSTRUME keyword */
990 tstatus = 0;
991 if (ffgkey(fptr, "INSTRUME", valstring, NULL, &tstatus) == 0 )
992 {
993 ffmkky("INSTRUME", valstring, comm, cptr, status);
994 length = strlen(cptr);
995 strncat(cptr, blanks, 80 - length); /* pad with blanks */
996 cptr += 80;
997 }
998
999 /* DETECTOR keyword */
1000 tstatus = 0;
1001 if (ffgkey(fptr, "DETECTOR", valstring, NULL, &tstatus) == 0 )
1002 {
1003 ffmkky("DETECTOR", valstring, comm, cptr, status);
1004 length = strlen(cptr);
1005 strncat(cptr, blanks, 80 - length); /* pad with blanks */
1006 cptr += 80;
1007 }
1008
1009 /* MJD-OBS keyword */
1010 tstatus = 0;
1011 if (ffgkey(fptr, "MJD-OBS", valstring, NULL, &tstatus) == 0 )
1012 {
1013 ffmkky("MJD-OBS", valstring, comm, cptr, status);
1014 length = strlen(cptr);
1015 strncat(cptr, blanks, 80 - length); /* pad with blanks */
1016 cptr += 80;
1017 }
1018
1019 /* DATE-OBS keyword */
1020 tstatus = 0;
1021 if (ffgkey(fptr, "DATE-OBS", valstring, NULL, &tstatus) == 0 )
1022 {
1023 ffmkky("DATE-OBS", valstring, comm, cptr, status);
1024 length = strlen(cptr);
1025 strncat(cptr, blanks, 80 - length); /* pad with blanks */
1026 cptr += 80;
1027 }
1028
1029 /* DATE keyword */
1030 tstatus = 0;
1031 if (ffgkey(fptr, "DATE", valstring, NULL, &tstatus) == 0 )
1032 {
1033 ffmkky("DATE", valstring, comm, cptr, status);
1034 length = strlen(cptr);
1035 strncat(cptr, blanks, 80 - length); /* pad with blanks */
1036 cptr += 80;
1037 }
1038
1039 strcat(cptr, "END");
1040 strncat(cptr, blanks, 77);
1041
1042 return(*status);
1043 }
1044