1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5
6 #include <stdlib.h>
7 #include <stdio.h>
8 #include <unistd.h>
9 #include <string.h>
10 #include <math.h>
11
12 #include "os-features.h"
13 #include "sip.h"
14 #include "sip-utils.h"
15 #include "sip_qfits.h"
16 #include "starutil.h"
17 #include "mathutil.h"
18 #include "boilerplate.h"
19 #include "errors.h"
20
21 const char* OPTIONS = "he:W:H:";
22
printHelp(char * progname)23 void printHelp(char* progname) {
24 BOILERPLATE_HELP_HEADER(stderr);
25 fprintf(stderr, "\nUsage: %s [options] <wcs-file>\n"
26 " [-e <extension>] Read from given HDU (default 0 = primary)\n"
27 " [-W <image width>] Set/override IMAGEW\n"
28 " [-H <image height>] Set/override IMAGEH\n"
29 "\n", progname);
30 }
31
32
main(int argc,char ** args)33 int main(int argc, char** args) {
34 int argchar;
35 char* progname = args[0];
36 char** inputfiles = NULL;
37 int ninputfiles = 0;
38 int ext = 0;
39 sip_t wcs;
40 double imw=0, imh=0;
41 double rac, decc;
42 double det, parity, orient, orientc;
43 int rah, ram, decsign, decd, decm;
44 double ras, decs;
45 char* units;
46 double pixscale;
47 double fldw, fldh;
48 double ramin, ramax, decmin, decmax;
49 double mxlo, mxhi, mylo, myhi;
50 double dm;
51 int merczoom;
52 char rastr[32];
53 char decstr[32];
54
55 while ((argchar = getopt (argc, args, OPTIONS)) != -1) {
56 switch (argchar) {
57 case 'e':
58 ext = atoi(optarg);
59 break;
60 case 'W':
61 imw = atof(optarg);
62 break;
63 case 'H':
64 imh = atof(optarg);
65 break;
66 case 'h':
67 default:
68 printHelp(progname);
69 exit(-1);
70 }
71 }
72 if (optind < argc) {
73 ninputfiles = argc - optind;
74 inputfiles = args + optind;
75 }
76 if (ninputfiles != 1) {
77 printHelp(progname);
78 exit(-1);
79 }
80
81 if (!sip_read_header_file_ext(inputfiles[0], ext, &wcs)) {
82 ERROR("failed to read WCS header from file %s, extension %i", inputfiles[0], ext);
83 return -1;
84 }
85
86 if (imw == 0)
87 imw = wcs.wcstan.imagew;
88 if (imh == 0)
89 imh = wcs.wcstan.imageh;
90 if ((imw == 0.0) || (imh == 0.0)) {
91 ERROR("failed to find IMAGE{W,H} in WCS file");
92 return -1;
93 }
94 // If W,H were set on the cmdline...
95 if (wcs.wcstan.imagew == 0)
96 wcs.wcstan.imagew = imw;
97 if (wcs.wcstan.imageh == 0)
98 wcs.wcstan.imageh = imh;
99
100 printf("crpix0 %.12g\n", wcs.wcstan.crpix[0]);
101 printf("crpix1 %.12g\n", wcs.wcstan.crpix[1]);
102 printf("crval0 %.12g\n", wcs.wcstan.crval[0]);
103 printf("crval1 %.12g\n", wcs.wcstan.crval[1]);
104 printf("ra_tangent %.12g\n", wcs.wcstan.crval[0]);
105 printf("dec_tangent %.12g\n", wcs.wcstan.crval[1]);
106 printf("pixx_tangent %.12g\n", wcs.wcstan.crpix[0]);
107 printf("pixy_tangent %.12g\n", wcs.wcstan.crpix[1]);
108
109 printf("imagew %.12g\n", imw);
110 printf("imageh %.12g\n", imh);
111
112 printf("cd11 %.12g\n", wcs.wcstan.cd[0][0]);
113 printf("cd12 %.12g\n", wcs.wcstan.cd[0][1]);
114 printf("cd21 %.12g\n", wcs.wcstan.cd[1][0]);
115 printf("cd22 %.12g\n", wcs.wcstan.cd[1][1]);
116
117 det = sip_det_cd(&wcs);
118 parity = (det >= 0 ? 1.0 : -1.0);
119 pixscale = sip_pixel_scale(&wcs);
120 printf("det %.12g\n", det);
121 printf("parity %i\n", (int)parity);
122 printf("pixscale %.12g\n", pixscale);
123
124 orient = sip_get_orientation(&wcs);
125 printf("orientation %.8g\n", orient);
126
127 sip_get_radec_center(&wcs, &rac, &decc);
128 printf("ra_center %.12g\n", rac);
129 printf("dec_center %.12g\n", decc);
130
131 // contributed by Rob Johnson, user rob at the domain whim.org, Nov 13, 2009
132 orientc = orient + rad2deg(atan(tan(deg2rad(rac - wcs.wcstan.crval[0])) * sin(deg2rad(wcs.wcstan.crval[1]))));
133 printf("orientation_center %.8g\n", orientc);
134
135 sip_get_radec_center_hms(&wcs, &rah, &ram, &ras, &decsign, &decd, &decm, &decs);
136 printf("ra_center_h %i\n", rah);
137 printf("ra_center_m %i\n", ram);
138 printf("ra_center_s %.12g\n", ras);
139 printf("dec_center_sign %i\n", decsign);
140 printf("dec_center_d %i\n", decd);
141 printf("dec_center_m %i\n", decm);
142 printf("dec_center_s %.12g\n", decs);
143
144 sip_get_radec_center_hms_string(&wcs, rastr, decstr);
145 printf("ra_center_hms %s\n", rastr);
146 printf("dec_center_dms %s\n", decstr);
147
148 // mercator
149 printf("ra_center_merc %.8g\n", ra2mercx(rac));
150 printf("dec_center_merc %.8g\n", dec2mercy(decc));
151
152 fldw = imw * pixscale;
153 fldh = imh * pixscale;
154 // area of the field, in square degrees.
155 printf("fieldarea %g\n", (arcsec2deg(fldw) * arcsec2deg(fldh)));
156
157 sip_get_field_size(&wcs, &fldw, &fldh, &units);
158 printf("fieldw %.4g\n", fldw);
159 printf("fieldh %.4g\n", fldh);
160 printf("fieldunits %s\n", units);
161
162 sip_get_radec_bounds(&wcs, 10, &ramin, &ramax, &decmin, &decmax);
163 printf("decmin %g\n", decmin);
164 printf("decmax %g\n", decmax);
165 printf("ramin %g\n", ramin);
166 printf("ramax %g\n", ramax);
167
168 // merc zoom level
169 mxlo = ra2mercx(ramax);
170 mxhi = ra2mercx(ramin);
171 mylo = dec2mercy(decmax);
172 myhi = dec2mercy(decmin);
173 printf("ra_min_merc %g\n", mxlo);
174 printf("ra_max_merc %g\n", mxhi);
175 printf("dec_min_merc %g\n", mylo);
176 printf("dec_max_merc %g\n", myhi);
177
178 dm = MAX(fabs(mxlo - mxhi), fabs(mylo - myhi));
179 printf("merc_diff %g\n", dm);
180 merczoom = 0 - (int)floor(log(dm) / log(2.0));
181 printf("merczoom %i\n", merczoom);
182 return 0;
183 }
184