1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
|
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "fitsio.h"
int voc_imgQParams (char *name, double *ra, double *dec, double *rad);
int main (int argc, char *argv[])
{
double ra, dec, rad;
int status = 0;
status = voc_imgQParams (argv[1], &ra, &dec, &rad);
printf ("RA=%f DEC=%f SZ=%f\n", ra, dec, rad);
}
int voc_imgQParams (char *name, double *ra, double *dec, double *rad)
{
fitsfile *fptr;
double xpos=0.0, ypos=0.0, xpix=0.0, ypix=0.0;
double xrval=0.0, yrval=0.0, xrpix=0.0, yrpix=0.0;
double xinc=0.0, yinc=0.0, rot=0.0;
double ll_x=0.0, ll_y=0.0, ur_x=0.0, ur_y=0.0, c_x=0.0, c_y=0.0;
long naxes[3], pcount, gcount;
int status=0, extend, simple, naxis, bitpix;
char ctype[5];
/* Open the FITS file.
*/
if (ffopen (&fptr, name, READWRITE, &status) > 0) {
fprintf (stderr, "Error: open status = %d\n", status);
exit (0);
}
/* Get the primary header keywords.
*/
ffghpr (fptr, 99, &simple, &bitpix, &naxis, naxes, &pcount,
&gcount, &extend, &status);
fprintf (stderr, "gcount = %d pcount = %d extend = %d \n", gcount, pcount, extend);
/* Get the header WCS keywords.
*/
ffgics (fptr, &xrval, &yrval, &xrpix,
&yrpix, &xinc, &yinc, &rot, ctype, &status);
if (status != 506)
fprintf (stderr, "Read WCS keywords with ffgics status = %d\n",status);
xpix = 0.5;
ypix = 0.5;
status = 0;
if (ffwldp (xpix, ypix, xrval, yrval, xrpix, yrpix, xinc, yinc, rot, ctype,
&ll_x, &ll_y, &status) > 0) {
fprintf (stderr, "Calculated sky coordinate with ffwldp status = %d\n",
status);
} else {
status = 0;
xpix = (double) naxes[0];
ypix = (double) naxes[1];
ffwldp (xpix, ypix, xrval, yrval, xrpix, yrpix, xinc, yinc, rot, ctype,
&ur_x, &ur_y, &status);
status = 0;
xpix = (double) naxes[0] / 2.;
ypix = (double) naxes[1] / 2.;
ffwldp (xpix, ypix, xrval, yrval, xrpix, yrpix, xinc, yinc, rot, ctype,
&c_x, &c_y, &status);
}
*ra = c_x;
*dec = c_y;
*rad = sqrt ((c_x - ll_x)*(c_x - ll_x) + (c_y - ll_y)*(c_y - ll_y));
#ifdef DBG_WCS
printf (" CRVAL1, CRVAL2 = %16.12f, %16.12f\n", xrval, yrval);
printf (" CRPIX1, CRPIX2 = %16.12f, %16.12f\n", xrpix, yrpix);
printf (" CDELT1, CDELT2 = %16.12f, %16.12f\n", xinc, yinc);
printf (" Rotation = %10.3f, CTYPE = %s\n", rot, ctype);
printf (" Pixels (%8.4f,%8.4f) --> (%11.6f, %11.6f) Sky\n",
xpix, ypix, ll_x, ll_y);
fprintf (stderr, "RA=%f DEC=%f SZ=%f\n", c_x, c_y,
sqrt ((c_x - ll_x)*(c_x - ll_x) + (c_y - ll_y)*(c_y - ll_y)));
ffxypx (xpos, ypos, xrval, yrval, xrpix, yrpix, xinc, yinc, rot, ctype,
&xpix, &ypix, &status);
printf ("Calculated pixel coordinate with ffxypx status = %d\n", status);
printf (" Sky (%11.6f, %11.6f) --> (%8.4f,%8.4f) Pixels\n",
xpos, ypos, xpix, ypix);
#endif
if (ffclos (fptr, &status) > 0) {
fprintf (stderr, "Error: close status = %d\n", status);
return (-1);
}
return (0);
}
|
