1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5 #include "sip.h"
6 #include "sip_qfits.h"
7 #include "anwcs.h"
8 #include "fitsioutils.h"
9
10 #include <stdio.h>
11 #include <unistd.h>
12 #include <string.h>
13
14 #include "cutest.h"
15
16 const char* tan1 = "SIMPLE = T / Standard FITS file BITPIX = 8 / ASCII or bytes array NAXIS = 0 / Minimal header EXTEND = T / There may be FITS ext CTYPE1 = 'RA---TAN' / TAN projection CTYPE2 = 'DEC--TAN' / TAN projection WCSAXES = 2 / no comment EQUINOX = 2000.0 / Equatorial coordinates definition (yr) LONPOLE = 180.0 / no comment LATPOLE = 0.0 / no comment CRVAL1 = 6.71529334958 / RA of reference point CRVAL2 = -72.7547910844 / DEC of reference point CRPIX1 = 477.760482899 / X reference pixel CRPIX2 = 361.955063329 / Y reference pixel CUNIT1 = 'deg ' / X pixel scale units CUNIT2 = 'deg ' / Y pixel scale units CD1_1 = -3.77685217356E-05 / Transformation matrix CD1_2 = -0.00273782095286 / no comment CD2_1 = 0.00273550077464 / no comment CD2_2 = -3.654231597E-05 / no comment IMAGEW = 1024 / Image width, in pixels. IMAGEH = 683 / Image height, in pixels. ";
17
18 const char* tan2 = "SIMPLE = T / Standard FITS file BITPIX = 8 / ASCII or bytes array NAXIS = 0 / Minimal header EXTEND = T / There may be FITS ext CTYPE1 = 'RA---TAN-SIP' / TAN (gnomic) projection + SIP distortions CTYPE2 = 'DEC--TAN-SIP' / TAN (gnomic) projection + SIP distortions WCSAXES = 2 / no comment EQUINOX = 2000.0 / Equatorial coordinates definition (yr) LONPOLE = 180.0 / no comment LATPOLE = 0.0 / no comment CRVAL1 = 6.71529334958 / RA of reference point CRVAL2 = -72.7547910844 / DEC of reference point CRPIX1 = 477.760482899 / X reference pixel CRPIX2 = 361.955063329 / Y reference pixel CUNIT1 = 'deg ' / X pixel scale units CUNIT2 = 'deg ' / Y pixel scale units CD1_1 = -3.77685217356E-05 / Transformation matrix CD1_2 = -0.00273782095286 / no comment CD2_1 = 0.00273550077464 / no comment CD2_2 = -3.654231597E-05 / no comment IMAGEW = 1024 / Image width, in pixels. IMAGEH = 683 / Image height, in pixels. A_ORDER = 2 / Polynomial order, axis 1 A_0_2 = 2.25528052918E-06 / no comment A_1_1 = -4.81359073875E-07 / no comment A_2_0 = 4.30780269448E-07 / no comment B_ORDER = 2 / Polynomial order, axis 2 B_0_2 = -1.03727867632E-06 / no comment B_1_1 = 6.9796860706E-07 / no comment B_2_0 = -3.7809448368E-07 / no comment AP_ORDER= 2 / Inv polynomial order, axis 1 AP_0_1 = -6.65805505351E-07 / no comment AP_0_2 = -2.2549438026E-06 / no comment AP_1_0 = 3.30484183954E-07 / no comment AP_1_1 = 4.80936510428E-07 / no comment AP_2_0 = -4.30764936375E-07 / no comment BP_ORDER= 2 / Inv polynomial order, axis 2 BP_0_1 = 4.50013020053E-07 / no comment BP_0_2 = 1.03706596388E-06 / no comment BP_1_0 = -2.70330536785E-07 / no comment BP_1_1 = -6.97662208285E-07 / no comment BP_2_0 = 3.78065361127E-07 / no comment ";
19
20 const char* sin1 = "SIMPLE = T / Standard FITS file BITPIX = 8 / ASCII or bytes array NAXIS = 0 / Minimal header EXTEND = T / There may be FITS ext CTYPE1 = 'RA---SIN' / SIN projection CTYPE2 = 'DEC--SIN' / SIN projection WCSAXES = 2 / no comment EQUINOX = 2000.0 / Equatorial coordinates definition (yr) LONPOLE = 180.0 / no comment LATPOLE = 0.0 / no comment CRVAL1 = 6.71529334958 / RA of reference point CRVAL2 = -72.7547910844 / DEC of reference point CRPIX1 = 477.760482899 / X reference pixel CRPIX2 = 361.955063329 / Y reference pixel CUNIT1 = 'deg ' / X pixel scale units CUNIT2 = 'deg ' / Y pixel scale units CD1_1 = -3.77685217356E-05 / Transformation matrix CD1_2 = -0.00273782095286 / no comment CD2_1 = 0.00273550077464 / no comment CD2_2 = -3.654231597E-05 / no comment IMAGEW = 1024 / Image width, in pixels. IMAGEH = 683 / Image height, in pixels. ";
21
22 const char* sin2 = "SIMPLE = T / Standard FITS file BITPIX = 8 / ASCII or bytes array NAXIS = 0 / Minimal header EXTEND = T / There may be FITS ext CTYPE1 = 'RA---SIN-SIP' / TAN (gnomic) projection + SIP distortions CTYPE2 = 'DEC--SIN-SIP' / TAN (gnomic) projection + SIP distortions WCSAXES = 2 / no comment EQUINOX = 2000.0 / Equatorial coordinates definition (yr) LONPOLE = 180.0 / no comment LATPOLE = 0.0 / no comment CRVAL1 = 6.71529334958 / RA of reference point CRVAL2 = -72.7547910844 / DEC of reference point CRPIX1 = 477.760482899 / X reference pixel CRPIX2 = 361.955063329 / Y reference pixel CUNIT1 = 'deg ' / X pixel scale units CUNIT2 = 'deg ' / Y pixel scale units CD1_1 = -3.77685217356E-05 / Transformation matrix CD1_2 = -0.00273782095286 / no comment CD2_1 = 0.00273550077464 / no comment CD2_2 = -3.654231597E-05 / no comment IMAGEW = 1024 / Image width, in pixels. IMAGEH = 683 / Image height, in pixels. A_ORDER = 2 / Polynomial order, axis 1 A_0_2 = 2.25528052918E-06 / no comment A_1_1 = -4.81359073875E-07 / no comment A_2_0 = 4.30780269448E-07 / no comment B_ORDER = 2 / Polynomial order, axis 2 B_0_2 = -1.03727867632E-06 / no comment B_1_1 = 6.9796860706E-07 / no comment B_2_0 = -3.7809448368E-07 / no comment AP_ORDER= 2 / Inv polynomial order, axis 1 AP_0_1 = -6.65805505351E-07 / no comment AP_0_2 = -2.2549438026E-06 / no comment AP_1_0 = 3.30484183954E-07 / no comment AP_1_1 = 4.80936510428E-07 / no comment AP_2_0 = -4.30764936375E-07 / no comment BP_ORDER= 2 / Inv polynomial order, axis 2 BP_0_1 = 4.50013020053E-07 / no comment BP_0_2 = 1.03706596388E-06 / no comment BP_1_0 = -2.70330536785E-07 / no comment BP_1_1 = -6.97662208285E-07 / no comment BP_2_0 = 3.78065361127E-07 / no comment ";
23
parse(const char * hdr,sip_t ** sip,anwcs_t ** anwcs,CuTest * tc)24 static void parse(const char* hdr, sip_t** sip, anwcs_t** anwcs, CuTest* tc) {
25 int len;
26 int hlen;
27 char* str;
28
29 fits_use_error_system();
30
31 len = strlen(hdr);
32 hlen = fits_bytes_needed(len + 80);
33 str = malloc(hlen + 1);
34 memcpy(str, hdr, len);
35 memset(str + len, ' ', hlen - len);
36 memcpy(str + hlen - 80, "END", 3);
37 str[hlen] = '\0';
38
39 //printf("Passing FITS header string: len %i, >>\n%s\n<<\n", hlen, str);
40
41 *sip = sip_from_string(str, hlen, NULL);
42 CuAssertPtrNotNull(tc, *sip);
43
44 *anwcs = anwcs_wcslib_from_string(str, hlen);
45 CuAssertPtrNotNull(tc, *anwcs);
46
47 free(str);
48 }
49
tst_rd2xy(double ra,double dec,double xtrue,double ytrue,anwcs_t * an1,sip_t * sip1,tan_t * tan1,CuTest * tc)50 static void tst_rd2xy(double ra, double dec, double xtrue, double ytrue,
51 anwcs_t* an1, sip_t* sip1, tan_t* tan1, CuTest* tc) {
52 double x, y;
53 anbool ok;
54 int err;
55
56 if (an1) {
57 x = y = 0.0;
58 err = anwcs_radec2pixelxy(an1, ra, dec, &x, &y);
59 printf("RA,Dec (%g, %g) -> x,y (%g, %g)\n", ra, dec, x, y);
60 CuAssertIntEquals(tc, 0, err);
61 CuAssertDblEquals(tc, xtrue, x, 1e-8);
62 CuAssertDblEquals(tc, ytrue, y, 1e-8);
63 }
64
65 if (sip1) {
66 x = y = 0.0;
67 ok = sip_radec2pixelxy(sip1, ra, dec, &x, &y);
68 printf("RA,Dec (%g, %g) -> x,y (%g, %g)\n", ra, dec, x, y);
69 CuAssertIntEquals(tc, TRUE, ok);
70 CuAssertDblEquals(tc, xtrue, x, 1e-8);
71 CuAssertDblEquals(tc, ytrue, y, 1e-8);
72 }
73
74 if (tan1) {
75 x = y = 0.0;
76 ok = tan_radec2pixelxy(tan1, ra, dec, &x, &y);
77 CuAssertIntEquals(tc, TRUE, ok);
78 CuAssertDblEquals(tc, xtrue, x, 1e-8);
79 CuAssertDblEquals(tc, ytrue, y, 1e-8);
80 }
81 }
82
tst_xy2rd(double x,double y,double ratrue,double dectrue,anwcs_t * an1,sip_t * sip1,tan_t * tan1,CuTest * tc)83 static void tst_xy2rd(double x, double y, double ratrue, double dectrue,
84 anwcs_t* an1, sip_t* sip1, tan_t* tan1, CuTest* tc) {
85 double ra, dec;
86 int err;
87
88 if (an1) {
89 ra = dec = 0.0;
90 err = anwcs_pixelxy2radec(an1, x, y, &ra, &dec);
91 printf("x,y (%g, %g) -> RA,Dec (%g, %g)\n", x, y, ra, dec);
92 CuAssertIntEquals(tc, 0, err);
93 CuAssertDblEquals(tc, ratrue, ra, 1e-8);
94 CuAssertDblEquals(tc, dectrue, dec, 1e-8);
95 }
96
97 if (sip1) {
98 ra = dec = 0.0;
99 sip_pixelxy2radec(sip1, x, y, &ra, &dec);
100 printf("x,y (%g, %g) -> RA,Dec (%g, %g)\n", x, y, ra, dec);
101 CuAssertDblEquals(tc, ratrue, ra, 1e-8);
102 CuAssertDblEquals(tc, dectrue, dec, 1e-8);
103 }
104
105 if (tan1) {
106 ra = dec = 0.0;
107 tan_pixelxy2radec(tan1, x, y, &ra, &dec);
108 printf("x,y (%g, %g) -> RA,Dec (%g, %g)\n", x, y, ra, dec);
109 CuAssertDblEquals(tc, ratrue, ra, 1e-8);
110 CuAssertDblEquals(tc, dectrue, dec, 1e-8);
111 }
112 }
113
test_tan1(CuTest * tc)114 void test_tan1(CuTest* tc) {
115 sip_t *sip1, *sip2, *sip3, *sip4;
116 anwcs_t * an1, * an2, * an3, * an4;
117
118 parse(tan1, &sip1, &an1, tc);
119 parse(tan2, &sip2, &an2, tc);
120 parse(sin1, &sip3, &an3, tc);
121 parse(sin2, &sip4, &an4, tc);
122
123 CuAssertIntEquals(tc, FALSE, sip1->wcstan.sin);
124 CuAssertIntEquals(tc, TRUE, sip3->wcstan.sin);
125
126 //anwcs_print(an1, stdout);
127 //sip_print(sip1);
128
129 /*
130 > wcs-rd2xy -w tan1.wcs -r 10.4 -d -74.1
131 RA,Dec (10.4000000000, -74.1000000000) -> pixel (-30.3363662787, 0.3458823213)
132 */
133
134 tst_rd2xy(10.4, -74.1, -30.3363662787, 0.3458823213,
135 an1, sip1, &(sip1->wcstan), tc);
136
137 /*
138 > wcs-xy2rd -w ../tan1.wcs -x 1 -y 1
139 Pixel (1.0000000000, 1.0000000000) -> RA,Dec (10.3701345009, -74.0147061244)
140 */
141 tst_xy2rd(1.0, 1.0, 10.3701345009, -74.0147061244,
142 an1, sip1, &(sip1->wcstan), tc);
143
144
145 /*
146 > wcs-rd2xy -w ../sin1.wcs -r 10.4 -d -74.1 -L
147 RA,Dec (10.4000000000, -74.1000000000) -> pixel (-30.1110266970, 0.5062550173)
148 */
149
150 tst_rd2xy(10.4, -74.1, -30.1110266970, 0.5062550173,
151 an3, sip3, &(sip3->wcstan), tc);
152
153 /*
154 > wcs-xy2rd -w ../sin1.wcs -x 1 -y 1 -L
155 Pixel (1.0000000000, 1.0000000000) -> RA,Dec (10.3717392951, -74.0152067481)
156 */
157
158 tst_xy2rd(1.0, 1.0, 10.3717392951, -74.0152067481,
159 an3, sip3, &(sip3->wcstan), tc);
160
161 // TAN-SIP / SIN-SIP
162
163 /*
164 > wcs-xy2rd -w ../tan2.wcs -x 1 -y 1
165 Pixel (1.0000000000, 1.0000000000) -> RA,Dec (10.3709060366, -74.0138433058)
166 */
167
168 tst_xy2rd(1.0, 1.0, 10.3709060366, -74.0138433058,
169 NULL, sip2, NULL, tc);
170
171
172 /*
173 > wcs-rd2xy -w tan2.wcs -r 10.4 -d -74.1
174 RA,Dec (10.4000000000, -74.1000000000) -> pixel (-30.6539962452, 0.4508839887)
175 */
176 tst_rd2xy(10.4, -74.1, -30.6539962452, 0.4508839887,
177 NULL, sip2, NULL, tc);
178
179
180 /*
181 > wcs-xy2rd -w sin2.wcs -x 1 -y 1
182 Pixel (1.0000000000, 1.0000000000) -> RA,Dec (10.3725100913, -74.0143432622)
183 */
184
185 tst_xy2rd(1.0, 1.0, 10.3725100913, -74.0143432622,
186 NULL, sip4, NULL, tc);
187
188 /*
189 > wcs-rd2xy -w sin2.wcs -r 10.4 -d -74.1
190 RA,Dec (10.4000000000, -74.1000000000) -> pixel (-30.4283749577, 0.6111635583)
191 */
192
193 tst_rd2xy(10.4, -74.1, -30.4283749577, 0.6111635583,
194 NULL, sip4, NULL, tc);
195
196
197 // SIP round-trips
198
199 double x, y, ra, dec, x2, y2;
200 anbool ok;
201
202 x = y = 100.0;
203
204 ra = dec = 0.0;
205 sip_pixelxy2radec(sip2, x, y, &ra, &dec);
206 printf("x,y (%g, %g) -> RA,Dec (%g, %g)\n", x, y, ra, dec);
207
208 x2 = y2 = 0.0;
209 ok = sip_radec2pixelxy(sip2, ra, dec, &x2, &y2);
210 printf("RA,Dec (%g, %g) -> x,y (%g, %g)\n", ra, dec, x2, y2);
211 CuAssertIntEquals(tc, TRUE, ok);
212
213 CuAssertDblEquals(tc, x, x2, 1e-3);
214 CuAssertDblEquals(tc, y, y2, 1e-3);
215
216
217 ra = dec = 0.0;
218 sip_pixelxy2radec(sip4, x, y, &ra, &dec);
219 printf("x,y (%g, %g) -> RA,Dec (%g, %g)\n", x, y, ra, dec);
220 x2 = y2 = 0.0;
221 ok = sip_radec2pixelxy(sip4, ra, dec, &x2, &y2);
222 printf("RA,Dec (%g, %g) -> x,y (%g, %g)\n", ra, dec, x2, y2);
223 CuAssertIntEquals(tc, TRUE, ok);
224 CuAssertDblEquals(tc, x, x2, 1e-3);
225 CuAssertDblEquals(tc, y, y2, 1e-3);
226
227 }
228
test_northpole(CuTest * tc)229 void test_northpole(CuTest* tc) {
230 tan_t thewcs;
231 tan_t *wcs;
232 double x, y, ra, dec;
233
234 //double dec_test = 89.7;
235 double dec_test = -89.7;
236
237 wcs = &thewcs;
238 thewcs.crval[0] = 180.;
239 //thewcs.crval[1] = +90.;
240 thewcs.crval[1] = -90.;
241 thewcs.crpix[0] = 1024.5;
242 thewcs.crpix[1] = 1024.5;
243 thewcs.cd[0][0] = -0.00152778;
244 thewcs.cd[0][1] = 0.;
245 thewcs.cd[1][0] = 0.;
246 thewcs.cd[1][1] = 0.00152778;
247 thewcs.imagew = 2048;
248 thewcs.imageh = 2048;
249 thewcs.sin = FALSE;
250
251 printf("\npixelxy2iwc\n");
252 ra = dec = 0.0;
253 tan_pixelxy2iwc(wcs, 800., 1024., &ra, &dec);
254 printf("800,1024 -> %.3f, %.3f\n", ra, dec);
255 ra = dec = 0.0;
256 tan_pixelxy2iwc(wcs, 1024., 800., &ra, &dec);
257 printf("1024,800 -> %.3f, %.3f\n", ra, dec);
258 ra = dec = 0.0;
259 tan_pixelxy2iwc(wcs, 1200., 1024., &ra, &dec);
260 printf("1200,1024 -> %.3f, %.3f\n", ra, dec);
261 ra = dec = 0.0;
262 tan_pixelxy2iwc(wcs, 1024., 1200., &ra, &dec);
263 printf("1024,1200 -> %.3f, %.3f\n", ra, dec);
264
265 printf("\niwc2pixelxy\n");
266 x = y = 0.0;
267 tan_iwc2pixelxy(wcs, +0.3, 0., &x, &y);
268 printf("+,0 -> %.1f, %.1f\n", x, y);
269 x = y = 0.0;
270 tan_iwc2pixelxy(wcs, 0., -0.3, &x, &y);
271 printf("0,- -> %.1f, %.1f\n", x, y);
272 x = y = 0.0;
273 tan_iwc2pixelxy(wcs, -0.3, 0., &x, &y);
274 printf("-,0 -> %.1f, %.1f\n", x, y);
275 x = y = 0.0;
276 tan_iwc2pixelxy(wcs, 0., +0.3, &x, &y);
277 printf("0,+ -> %.1f, %.1f\n", x, y);
278
279 printf("\niwc2radec\n");
280 x = y = 0.0;
281 tan_iwc2radec(wcs, +0.3, 0., &ra, &dec);
282 printf("+,0 -> %.1f, %.1f\n", ra, dec);
283 x = y = 0.0;
284 tan_iwc2radec(wcs, 0., -0.3, &ra, &dec);
285 printf("0,- -> %.1f, %.1f\n", ra, dec);
286 x = y = 0.0;
287 tan_iwc2radec(wcs, -0.3, 0., &ra, &dec);
288 printf("-,0 -> %.1f, %.1f\n", ra, dec);
289 x = y = 0.0;
290 tan_iwc2radec(wcs, 0., +0.3, &ra, &dec);
291 printf("0,+ -> %.1f, %.1f\n", ra, dec);
292
293 printf("\nradec2iwc\n");
294 x = y = 0.0;
295 tan_radec2iwc(wcs, 0., dec_test, &x, &y);
296 printf("0, %.1f -> %.3f, %.3f\n", dec_test, x, y);
297 x = y = 0.0;
298 tan_radec2iwc(wcs, 90., dec_test, &x, &y);
299 printf("90, %.1f -> %.3f, %.3f\n", dec_test, x, y);
300 x = y = 0.0;
301 tan_radec2iwc(wcs, 180., dec_test, &x, &y);
302 printf("180, %.1f -> %.3f, %.3f\n", dec_test, x, y);
303 x = y = 0.0;
304 tan_radec2iwc(wcs, 270., dec_test, &x, &y);
305 printf("270, %.1f -> %.3f, %.3f\n", dec_test, x, y);
306
307
308 printf("\npixelxy2radec\n");
309 ra = dec = 0.0;
310 tan_pixelxy2radec(wcs, 800., 1024., &ra, &dec);
311 printf("800,1024 -> %.3f, %.3f\n", ra, dec);
312 ra = dec = 0.0;
313 tan_pixelxy2radec(wcs, 1024., 800., &ra, &dec);
314 printf("1024,800 -> %.3f, %.3f\n", ra, dec);
315 ra = dec = 0.0;
316 tan_pixelxy2radec(wcs, 1200., 1024., &ra, &dec);
317 printf("1200,1024 -> %.3f, %.3f\n", ra, dec);
318 ra = dec = 0.0;
319 tan_pixelxy2radec(wcs, 1024., 1200., &ra, &dec);
320 printf("1024,1200 -> %.3f, %.3f\n", ra, dec);
321
322 printf("\nradec2pixelxy\n");
323 x = y = 0.0;
324 tan_radec2pixelxy(wcs, 0., dec_test, &x, &y);
325 printf("0, %.1f -> %.1f, %.1f\n", dec_test, x, y);
326 x = y = 0.0;
327 tan_radec2pixelxy(wcs, 90., dec_test, &x, &y);
328 printf("90, %.1f -> %.1f, %.1f\n", dec_test, x, y);
329 x = y = 0.0;
330 tan_radec2pixelxy(wcs, 180., dec_test, &x, &y);
331 printf("180, %.1f -> %.1f, %.1f\n", dec_test, x, y);
332 x = y = 0.0;
333 tan_radec2pixelxy(wcs, 270., dec_test, &x, &y);
334 printf("270, %.1f -> %.1f, %.1f\n", dec_test, x, y);
335
336
337 /*
338 printf("\n\n\ncrval[0] = 0.\n");
339 thewcs.crval[0] = 0.;
340 printf("\niwc2radec\n");
341 x = y = 0.0;
342 tan_iwc2radec(wcs, +0.3, 0., &ra, &dec);
343 printf("+,0 -> %.1f, %.1f\n", ra, dec);
344 x = y = 0.0;
345 tan_iwc2radec(wcs, 0., -0.3, &ra, &dec);
346 printf("0,- -> %.1f, %.1f\n", ra, dec);
347 x = y = 0.0;
348 tan_iwc2radec(wcs, -0.3, 0., &ra, &dec);
349 printf("-,0 -> %.1f, %.1f\n", ra, dec);
350 x = y = 0.0;
351 tan_iwc2radec(wcs, 0., +0.3, &ra, &dec);
352 printf("0,+ -> %.1f, %.1f\n", ra, dec);
353 */
354
355 }
356
357