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 <math.h>
7 #include <stdio.h>
8 #include <assert.h>
9
10 #include "os-features.h"
11 #include "cutest.h"
12 #include "starutil.h"
13 #include "healpix.h"
14 #include "healpix-utils.h"
15 #include "bl.h"
16
square(double x)17 static double square(double x) {
18 return x*x;
19 }
20
test_hp_rangesearch(CuTest * ct)21 void test_hp_rangesearch(CuTest* ct) {
22 double ra = 120.12998471831297;
23 double dec = 37.72245880730694;
24 double rad = 0.2242483527520963;
25 int nside = 8;
26 il* hps = il_new(4);
27 int i;
28 healpix_rangesearch_radec(ra, dec, rad, nside, hps);
29 for (i=0; i<il_size(hps); i++) {
30 int hp = il_get(hps, i);
31 printf("Healpix %i\n", hp);
32 }
33 CuAssertIntEquals(ct, il_size(hps), 2);
34 CuAssertIntEquals(ct, il_get(hps,0), 84);
35 CuAssertIntEquals(ct, il_get(hps,1), 85);
36 }
37
test_side_length(CuTest * ct)38 void test_side_length(CuTest* ct) {
39 double hp;
40 double len = healpix_side_length_arcmin(1);
41 CuAssertDblEquals(ct, 3517.9, len, 0.1);
42 hp = healpix_nside_for_side_length_arcmin(len);
43 CuAssertDblEquals(ct, 1.0, hp, 0.001);
44
45 len = healpix_side_length_arcmin(2);
46 CuAssertDblEquals(ct, 1758.969, len, 0.001);
47 hp = healpix_nside_for_side_length_arcmin(len);
48 CuAssertDblEquals(ct, 2.0, hp, 0.001);
49 }
50
add_plot_xyz_point(double * xyz)51 static void add_plot_xyz_point(double* xyz) {
52 double ra,dec;
53 xyzarr2radecdeg(xyz, &ra, &dec);
54 fprintf(stderr, "xp.append(%g)\n", ra);
55 fprintf(stderr, "yp.append(%g)\n", dec);
56 }
57
add_plot_point(int hp,int nside,double dx,double dy)58 static void add_plot_point(int hp, int nside, double dx, double dy) {
59 double xyz[3];
60 healpix_to_xyzarr(hp, nside, dx, dy, xyz);
61 add_plot_xyz_point(xyz);
62 }
63
plot_point(int hp,int nside,double dx,double dy,char * style)64 void plot_point(int hp, int nside, double dx, double dy, char* style) {
65 fprintf(stderr, "xp=[]\n");
66 fprintf(stderr, "yp=[]\n");
67 add_plot_point(hp, nside, dx, dy);
68 fprintf(stderr, "plot(xp, yp, '%s')\n", style);
69 }
70
plot_xyz_point(double * xyz,char * style)71 void plot_xyz_point(double* xyz, char* style) {
72 fprintf(stderr, "xp=[]\n");
73 fprintf(stderr, "yp=[]\n");
74 add_plot_xyz_point(xyz);
75 fprintf(stderr, "plot(xp, yp, '%s')\n", style);
76 }
77
plot_hp_boundary(int hp,int nside,double start,double step,char * style)78 static void plot_hp_boundary(int hp, int nside, double start, double step, char* style) {
79 double dx, dy;
80 fprintf(stderr, "xp=[]\n");
81 fprintf(stderr, "yp=[]\n");
82 dy = 0.0;
83 for (dx=start; dx<=1.0; dx+=step)
84 add_plot_point(hp, nside, dx, dy);
85 dx = 1.0;
86 for (dy=start; dy<=1.0; dy+=step)
87 add_plot_point(hp, nside, dx, dy);
88 dy = 1.0;
89 for (dx=1.0-start; dx>=0.0; dx-=step)
90 add_plot_point(hp, nside, dx, dy);
91 dx = 0.0;
92 for (dy=1.0-start; dy>=0.0; dy-=step)
93 add_plot_point(hp, nside, dx, dy);
94 dy = 0.0;
95 add_plot_point(hp, nside, dx, dy);
96 fprintf(stderr, "xp,yp = wrapxy(xp,yp)\nplot(xp, yp, '%s')\n", style);
97 }
98
hpmap(int nside,const char * fn)99 static void hpmap(int nside, const char* fn) {
100 #if 0
101 int nhp;
102 #endif
103 double xyz[3];
104 double range;
105 int hps[9];
106 int i;
107 int hp;
108 double dx, dy;
109
110 // pick a point on the edge.
111 //hp = 8;
112 hp = 9;
113 dx = 0.95;
114 dy = 0.0;
115 range = 0.1;
116 /*
117 hp = 6;
118 dx = 0.05;
119 dy = 0.95;
120 range = 0.1;
121 */
122 healpix_to_xyzarr(hp, nside, dx, dy, xyz);
123
124 for (i=0; i<12*nside*nside; i++) {
125 plot_hp_boundary(i, nside, 0.005, 0.01, "b-");
126 }
127
128 #if 0
129 nhp = healpix_get_neighbours_within_range(xyz, range, hps, nside);
130 assert(nhp >= 1);
131 assert(nhp <= 9);
132 #else
133 (void)healpix_get_neighbours_within_range(xyz, range, hps, nside);
134 #endif
135
136 /*
137 for (i=0; i<nhp; i++) {
138 printf("in range: %i\n", hps[i]);
139 plot_hp_boundary(hps[i], nside, 0.005, 0.01, "b-");
140 }
141 plot_hp_boundary(hp, nside, 0.005, 0.01, "k-");
142 plot_point(hp, nside, dx, dy, "r.");
143 */
144
145 for (i=0; i<12*nside*nside; i++) {
146 fprintf(stderr, "xp=[]; yp=[]\n");
147 add_plot_point(i, nside, 0.5, 0.5);
148 //fprintf(stderr, "text(xp[0], yp[0], '%i', color='b', horizontalalignment='center', verticalalignment='center', bbox=dict(facecolor='w', edgecolor='w'))\n", i);
149 fprintf(stderr, "text(xp[0], yp[0], '%i', color='b', horizontalalignment='center', verticalalignment='center')\n", i);
150 }
151
152 fprintf(stderr, "axis((360, 0, -90, 90))\n");
153 fprintf(stderr, "xlabel('RA (deg)')\n");
154 fprintf(stderr, "ylabel('Dec (deg)')\n");
155 fprintf(stderr, "title('healpixels: nside=%i')\n", nside);
156 fprintf(stderr, "savefig('%s')\n", fn);
157 fprintf(stderr, "clf()\n");
158 }
159
testX_make_map(CuTest * ct)160 void testX_make_map(CuTest* ct) {
161 fprintf(stderr, "%s", "from pylab import *\n");
162 fprintf(stderr, "clf()\n"
163 "def wrapxy(x,y):\n"
164 " lastx = x[0]\n"
165 " lasty = y[0]\n"
166 " outx = [lastx]\n"
167 " outy = [lasty]\n"
168 " for xx,yy in zip(x[1:],y[1:]):\n"
169 " if (xx-lastx)**2 + (yy - lasty)**2 > 1.:\n"
170 " if xx < 180:\n"
171 " xx += 360\n"
172 " else:\n"
173 " xx -= 360\n"
174 " outx.append(xx)\n"
175 " outy.append(yy)\n"
176 " lastx = xx\n"
177 " lasty = yy\n"
178 " return (array(outx),array(outy))\n"
179 );
180
181 hpmap(1, "hp.png");
182 hpmap(2, "hp2.png");
183 }
184
185
tst_xyztohpf(CuTest * ct,int hp,int nside,double dx,double dy)186 int tst_xyztohpf(CuTest* ct,
187 int hp, int nside, double dx, double dy) {
188 double x,y,z;
189 double outdx, outdy;
190 int outhp;
191 double outx,outy,outz;
192 double dist;
193 healpix_to_xyz(hp, nside, dx, dy, &x, &y, &z);
194 outhp = xyztohealpixf(x, y, z, nside, &outdx, &outdy);
195 healpix_to_xyz(outhp, nside, outdx, outdy, &outx, &outy, &outz);
196 dist = sqrt(MAX(0, square(x-outx) + square(y-outy) + square(z-outz)));
197 printf("true/computed:\n"
198 "hp: %i / %i\n"
199 "dx: %.20g / %.20g\n"
200 "dy: %.20g / %.20g\n"
201 "x: %g / %g\n"
202 "y: %g / %g\n"
203 "z: %g / %g\n"
204 "dist: %g\n\n",
205 hp, outhp, dx, outdx, dy, outdy,
206 x, outx, y, outy, z, outz, dist);
207
208 if (dist > 1e-6) {
209 double a, b;
210 double outa, outb;
211 a = xy2ra(x,y) / (2.0 * M_PI);
212 b = z2dec(z) / (M_PI);
213 outa = xy2ra(outx, outy) / (2.0 * M_PI);
214 outb = z2dec(outz) / (M_PI);
215 fprintf(stderr,
216 "plot([%g, %g],[%g, %g],'r.-')\n", a, outa, b, outb);
217 fprintf(stderr,
218 "text(%g, %g, \"(%g,%g)\")\n",
219 a, b, dx, dy);
220 }
221
222 CuAssertIntEquals(ct, 1, (dist < 1e-6)?1:0);
223 return (dist > 1e-6);
224 }
225
tEst_xyztohpf(CuTest * ct)226 void tEst_xyztohpf(CuTest* ct) {
227 double dx, dy;
228 int hp;
229 int nside;
230 double step = 0.1;
231 double a, b;
232 nside = 1;
233
234 fprintf(stderr, "%s", "from pylab import plot,text,savefig,clf\n");
235 fprintf(stderr, "clf()\n");
236
237 /*
238 Plot the grid of healpixes with dx,dy=step steps.
239 */
240 step = 0.25;
241 //for (hp=0; hp<12*nside*nside; hp++) {
242 for (hp=0; hp<1*nside*nside; hp++) {
243 double x,y,z;
244 for (dx=0.0; dx<=1.05; dx+=step) {
245 fprintf(stderr, "xp=[]\n");
246 fprintf(stderr, "yp=[]\n");
247 for (dy=0.0; dy<=1.05; dy+=step) {
248 healpix_to_xyz(hp, nside, dx, dy, &x, &y, &z);
249 a = xy2ra(x,y) / (2.0 * M_PI);
250 b = z2dec(z) / (M_PI);
251 fprintf(stderr, "xp.append(%g)\n", a);
252 fprintf(stderr, "yp.append(%g)\n", b);
253 }
254 fprintf(stderr, "plot(xp, yp, 'k-')\n");
255 }
256 for (dy=0.0; dy<=1.05; dy+=step) {
257 fprintf(stderr, "xp=[]\n");
258 fprintf(stderr, "yp=[]\n");
259 for (dx=0.0; dx<=1.0; dx+=step) {
260 healpix_to_xyz(hp, nside, dx, dy, &x, &y, &z);
261 a = xy2ra(x,y) / (2.0 * M_PI);
262 b = z2dec(z) / (M_PI);
263 fprintf(stderr, "xp.append(%g)\n", a);
264 fprintf(stderr, "yp.append(%g)\n", b);
265 }
266 fprintf(stderr, "plot(xp, yp, 'k-')\n");
267 }
268 }
269
270 step = 0.5;
271 /*
272 Plot places where the conversion screws up.
273 */
274 for (hp=0; hp<12*nside*nside; hp++) {
275 for (dx=0.0; dx<=1.01; dx+=step) {
276 for (dy=0.0; dy<=1.01; dy+=step) {
277 tst_xyztohpf(ct, hp, nside, dx, dy);
278 }
279 }
280 }
281 fprintf(stderr, "savefig('plot.png')\n");
282
283 }
284
tst_neighbours(CuTest * ct,int pix,int * true_neigh,int true_nn,int Nside)285 static void tst_neighbours(CuTest* ct, int pix, int* true_neigh, int true_nn,
286 int Nside) {
287 int neigh[8];
288 int nn;
289 int i;
290 for (i=0; i<8; i++)
291 neigh[i] = -1;
292 nn = healpix_get_neighbours(pix, neigh, Nside);
293 /*
294 printf("true(%i) : [ ", pix);
295 for (i=0; i<true_nn; i++)
296 printf("%u, ", true_neigh[i]);
297 printf("]\n");
298 printf("got (%i) : [ ", pix);
299 for (i=0; i<nn; i++)
300 printf("%u, ", neigh[i]);
301 printf("]\n");
302 */
303 CuAssertIntEquals(ct, nn, true_nn);
304
305 for (i=0; i<true_nn; i++)
306 CuAssertIntEquals(ct, true_neigh[i], neigh[i]);
307 }
308
tst_nested(CuTest * ct,int pix,int * true_neigh,int true_nn,int Nside)309 static void tst_nested(CuTest* ct, int pix, int* true_neigh, int true_nn,
310 int Nside) {
311 int i;
312 int truexy[8];
313 int xypix;
314
315 /*
316 printf("nested true(%i) : [ ", pix);
317 for (i=0; i<true_nn; i++)
318 printf("%u ", true_neigh[i]);
319 printf("]\n");
320 */
321
322 CuAssert(ct, "true_nn <= 8", true_nn <= 8);
323 for (i=0; i<true_nn; i++) {
324 truexy[i] = healpix_nested_to_xy(true_neigh[i], Nside);
325 CuAssertIntEquals(ct, true_neigh[i], healpix_xy_to_nested(truexy[i], Nside));
326 }
327 xypix = healpix_nested_to_xy(pix, Nside);
328 CuAssertIntEquals(ct, pix, healpix_xy_to_nested(xypix, Nside));
329
330 tst_neighbours(ct, xypix, truexy, true_nn, Nside);
331 }
332
print_node(double z,double phi,int Nside)333 void print_node(double z, double phi, int Nside) {
334 double ra, dec;
335 int hp;
336 int nn;
337 int neigh[8];
338 int k;
339
340 double scale = 10.0;
341
342 ra = phi;
343 dec = asin(z);
344 while (ra < 0.0)
345 ra += 2.0 * M_PI;
346 while (ra > 2.0 * M_PI)
347 ra -= 2.0 * M_PI;
348
349 // find its healpix.
350 hp = radectohealpix(ra, dec, Nside);
351 // find its neighbourhood.
352 nn = healpix_get_neighbours(hp, neigh, Nside);
353 fprintf(stderr, " N%i [ label=\"%i\", pos=\"%g,%g!\" ];\n", hp, hp,
354 scale * ra/M_PI, scale * z);
355 for (k=0; k<nn; k++) {
356 fprintf(stderr, " N%i -- N%i\n", hp, neigh[k]);
357 }
358 }
359
test_healpix_distance_to_radec(CuTest * ct)360 void test_healpix_distance_to_radec(CuTest *ct) {
361 double d;
362 double rd[2];
363
364 d = healpix_distance_to_radec(4, 1, 0, 0, NULL);
365 CuAssertDblEquals(ct, 0, d, 0);
366 d = healpix_distance_to_radec(4, 1, 45, 0, NULL);
367 CuAssertDblEquals(ct, 0, d, 0);
368 d = healpix_distance_to_radec(4, 1, 45+1, 0, NULL);
369 CuAssertDblEquals(ct, 1, d, 1e-9);
370 d = healpix_distance_to_radec(4, 1, 45+1, 0+1, NULL);
371 CuAssertDblEquals(ct, 1.414, d, 1e-3);
372
373 d = healpix_distance_to_radec(4, 1, 45+10, 0, NULL);
374 CuAssertDblEquals(ct, 10, d, 1e-9);
375
376 // top corner
377 d = healpix_distance_to_radec(4, 1, 0, rad2deg(asin(2.0/3.0)), NULL);
378 CuAssertDblEquals(ct, 0, d, 1e-9);
379
380 d = healpix_distance_to_radec(4, 1, 0, 1 + rad2deg(asin(2.0/3.0)), NULL);
381 CuAssertDblEquals(ct, 1, d, 1e-9);
382
383 d = healpix_distance_to_radec(4, 1, -45-10, -10, NULL);
384 CuAssertDblEquals(ct, 14.106044, d, 1e-6);
385
386 d = healpix_distance_to_radec(10, 1, 225, 5, NULL);
387 CuAssertDblEquals(ct, 5, d, 1e-6);
388
389 d = healpix_distance_to_radec(44, 2, 300, -50, NULL);
390 CuAssertDblEquals(ct, 3.007643, d, 1e-6);
391
392 d = healpix_distance_to_radec(45, 2, 310, -50, NULL);
393 CuAssertDblEquals(ct, 1.873942, d, 1e-6);
394
395 // south-polar hp, north pole.
396 d = healpix_distance_to_radec(36, 2, 180, 90, NULL);
397 // The hp corner is -41.8 deg; add 90.
398 CuAssertDblEquals(ct, 131.810, d, 1e-3);
399
400 // just south of equator to nearly across the sphere
401 d = healpix_distance_to_radec(35, 2, 225, 20, NULL);
402 // this one actually has the midpoint further than A and B.
403 CuAssertDblEquals(ct, 158.189685, d, 1e-6);
404
405 /*
406 xyz[0] = -100.0;
407 ra = dec = -1.0;
408 d = healpix_distance_to_xyz(4, 1, 0, 0, xyz);
409 xyzarr2radecdeg(xyz, &ra, &dec);
410 CuAssertDblEquals(ct, 0, ra, 0);
411 CuAssertDblEquals(ct, 0, dec, 0);
412 */
413 rd[0] = rd[1] = -1.0;
414 d = healpix_distance_to_radec(4, 1, 0, 0, rd);
415 CuAssertDblEquals(ct, 0, rd[0], 0);
416 CuAssertDblEquals(ct, 0, rd[1], 0);
417
418
419 /*
420 xyz[0] = -100.0;
421 ra = dec = -1.0;
422 d = healpix_distance_to_xyz(4, 1, 45, 0, xyz);
423 xyzarr2radecdeg(xyz, &ra, &dec);
424 CuAssertDblEquals(ct, 45, ra, 0);
425 CuAssertDblEquals(ct, 0, dec, 0);
426 */
427
428 rd[0] = rd[1] = -1.0;
429 d = healpix_distance_to_radec(4, 1, 45+1, 0, rd);
430 CuAssertDblEquals(ct, 45, rd[0], 1e-12);
431 CuAssertDblEquals(ct, 0, rd[1], 1e-8);
432
433 d = healpix_distance_to_radec(4, 1, 45+1, 0+1, rd);
434 CuAssertDblEquals(ct, 45, rd[0], 1e-12);
435 CuAssertDblEquals(ct, 0, rd[1], 0);
436 // really??
437
438 d = healpix_distance_to_radec(4, 1, 20, 25, rd);
439 CuAssertDblEquals(ct, d, 2.297298, 1e-6);
440 CuAssertDblEquals(ct, 18.200995, rd[0], 1e-6);
441 CuAssertDblEquals(ct, 23.392159, rd[1], 1e-6);
442
443 }
444
test_healpix_neighbours(CuTest * ct)445 void test_healpix_neighbours(CuTest *ct) {
446 int n0[] = { 1,3,2,71,69,143,90,91 };
447 int n5[] = { 26,27,7,6,4,94,95 };
448 int n13[] = { 30,31,15,14,12,6,7,27 };
449 int n15[] = { 31,47,63,61,14,12,13,30 };
450 int n30[] = { 31,15,13,7,27,25,28,29 };
451 int n101[] = { 32,34,103,102,100,174,175,122 };
452 int n127[] = { 58,37,36,126,124,125,56 };
453 int n64[] = { 65,67,66,183,181,138,139 };
454 int n133[] = { 80,82,135,134,132,152,154 };
455 int n148[] = { 149,151,150,147,145,162,168,170 };
456 int n160[] = { 161,163,162,145,144,128,176,178 };
457 int n24[] = { 25,27,26,95,93,87,18,19 };
458 int n42[] = { 43,23,21,111,109,40,41 };
459 int n59[] = { 62,45,39,37,58,56,57,60 };
460 int n191[] = { 74,48,117,116,190,188,189,72 };
461 int n190[] = { 191,117,116,113,187,185,188,189 };
462 int n186[] = { 187,113,112,165,164,184,185 };
463 int n184[] = { 185,187,186,165,164,161,178,179 };
464
465 // These were taken (IIRC) from the Healpix paper, so the healpix
466 // numbers are all in the NESTED scheme.
467
468 tst_nested(ct, 0, n0, sizeof(n0) /sizeof(int), 4);
469 tst_nested(ct, 5, n5, sizeof(n5) /sizeof(int), 4);
470 tst_nested(ct, 13, n13, sizeof(n13) /sizeof(int), 4);
471 tst_nested(ct, 15, n15, sizeof(n15) /sizeof(int), 4);
472 tst_nested(ct, 30, n30, sizeof(n30) /sizeof(int), 4);
473 tst_nested(ct, 101, n101, sizeof(n101)/sizeof(int), 4);
474 tst_nested(ct, 127, n127, sizeof(n127)/sizeof(int), 4);
475 tst_nested(ct, 64, n64, sizeof(n64) /sizeof(int), 4);
476 tst_nested(ct, 133, n133, sizeof(n133)/sizeof(int), 4);
477 tst_nested(ct, 148, n148, sizeof(n148)/sizeof(int), 4);
478 tst_nested(ct, 160, n160, sizeof(n160)/sizeof(int), 4);
479 tst_nested(ct, 24, n24, sizeof(n24) /sizeof(int), 4);
480 tst_nested(ct, 42, n42, sizeof(n42) /sizeof(int), 4);
481 tst_nested(ct, 59, n59, sizeof(n59) /sizeof(int), 4);
482 tst_nested(ct, 191, n191, sizeof(n191)/sizeof(int), 4);
483 tst_nested(ct, 190, n190, sizeof(n190)/sizeof(int), 4);
484 tst_nested(ct, 186, n186, sizeof(n186)/sizeof(int), 4);
485 tst_nested(ct, 184, n184, sizeof(n184)/sizeof(int), 4);
486 }
487
488 /*
489 void pnprime_to_xy(int, int*, int*, int);
490 int xy_to_pnprime(int, int, int);
491
492 void tst_healpix_pnprime_to_xy(CuTest *ct) {
493 int px,py;
494 pnprime_to_xy(6, &px, &py, 3);
495 CuAssertIntEquals(ct, px, 2);
496 CuAssertIntEquals(ct, py, 0);
497 pnprime_to_xy(8, &px, &py, 3);
498 CuAssertIntEquals(ct, px, 2);
499 CuAssertIntEquals(ct, py, 2);
500 pnprime_to_xy(0, &px, &py, 3);
501 CuAssertIntEquals(ct, px, 0);
502 CuAssertIntEquals(ct, py, 0);
503 pnprime_to_xy(2, &px, &py, 3);
504 CuAssertIntEquals(ct, px, 0);
505 CuAssertIntEquals(ct, py, 2);
506 pnprime_to_xy(4, &px, &py, 3);
507 CuAssertIntEquals(ct, px, 1);
508 CuAssertIntEquals(ct, py, 1);
509 }
510
511 void tst_healpix_xy_to_pnprime(CuTest *ct) {
512 CuAssertIntEquals(ct, xy_to_pnprime(0,0,3), 0);
513 CuAssertIntEquals(ct, xy_to_pnprime(1,0,3), 3);
514 CuAssertIntEquals(ct, xy_to_pnprime(2,0,3), 6);
515 CuAssertIntEquals(ct, xy_to_pnprime(0,1,3), 1);
516 CuAssertIntEquals(ct, xy_to_pnprime(1,1,3), 4);
517 CuAssertIntEquals(ct, xy_to_pnprime(2,1,3), 7);
518 CuAssertIntEquals(ct, xy_to_pnprime(0,2,3), 2);
519 CuAssertIntEquals(ct, xy_to_pnprime(1,2,3), 5);
520 CuAssertIntEquals(ct, xy_to_pnprime(2,2,3), 8);
521 }
522 */
print_test_healpix_output(int Nside)523 void print_test_healpix_output(int Nside) {
524
525 int i, j;
526 double z;
527 double phi;
528 fprintf(stderr, "graph Nside4 {\n");
529
530 // north polar
531 for (i=1; i<=Nside; i++) {
532 for (j=1; j<=(4*i); j++) {
533 // find the center of the pixel in ring i
534 // and longitude j.
535 z = 1.0 - square((double)i / (double)Nside)/3.0;
536 phi = M_PI / (2.0 * i) * ((double)j - 0.5);
537 fprintf(stderr, " // North polar, i=%i, j=%i. z=%g, phi=%g\n", i, j, z, phi);
538 print_node(z, phi, Nside);
539 }
540 }
541 // south polar
542 for (i=1; i<=Nside; i++) {
543 for (j=1; j<=(4*i); j++) {
544 z = 1.0 - square((double)i / (double)Nside)/3.0;
545 z *= -1.0;
546 phi = M_PI / (2.0 * i) * ((double)j - 0.5);
547 fprintf(stderr, " // South polar, i=%i, j=%i. z=%g, phi=%g\n", i, j, z, phi);
548 print_node(z, phi, Nside);
549 }
550 }
551 // north equatorial
552 for (i=Nside+1; i<=2*Nside; i++) {
553 for (j=1; j<=(4*Nside); j++) {
554 int s;
555 z = 4.0/3.0 - 2.0 * i / (3.0 * Nside);
556 s = (i - Nside + 1) % 2;
557 s = (s + 2) % 2;
558 phi = M_PI / (2.0 * Nside) * ((double)j - (double)s / 2.0);
559 fprintf(stderr, " // North equatorial, i=%i, j=%i. z=%g, phi=%g, s=%i\n", i, j, z, phi, s);
560 print_node(z, phi, Nside);
561 }
562 }
563 // south equatorial
564 for (i=Nside+1; i<2*Nside; i++) {
565 for (j=1; j<=(4*Nside); j++) {
566 int s;
567 z = 4.0/3.0 - 2.0 * i / (3.0 * Nside);
568 z *= -1.0;
569 s = (i - Nside + 1) % 2;
570 s = (s + 2) % 2;
571 phi = M_PI / (2.0 * Nside) * ((double)j - s / 2.0);
572 fprintf(stderr, " // South equatorial, i=%i, j=%i. z=%g, phi=%g, s=%i\n", i, j, z, phi, s);
573 print_node(z, phi, Nside);
574 }
575 }
576
577 fprintf(stderr, " node [ shape=point ]\n");
578 fprintf(stderr, " C0 [ pos=\"0,-10!\" ];\n");
579 fprintf(stderr, " C1 [ pos=\"20,-10!\" ];\n");
580 fprintf(stderr, " C2 [ pos=\"20,10!\" ];\n");
581 fprintf(stderr, " C3 [ pos=\"0,10!\" ];\n");
582 fprintf(stderr, " C0 -- C1 -- C2 -- C3 -- C0\n");
583 fprintf(stderr, "}\n");
584 }
585
print_healpix_grid(int Nside)586 void print_healpix_grid(int Nside) {
587 int i;
588 int j;
589 int N = 500;
590
591 fprintf(stderr, "x%i=[", Nside);
592 for (i=0; i<N; i++) {
593 for (j=0; j<N; j++) {
594 fprintf(stderr, "%i ", radectohealpix(i*2*M_PI/N, M_PI*(j-N/2)/N, Nside));
595 }
596 fprintf(stderr, ";");
597 }
598 fprintf(stderr, "];\n\n");
599 fflush(stderr);
600 }
601
print_healpix_borders(int Nside)602 void print_healpix_borders(int Nside) {
603 int i;
604 int j;
605 int N = 1;
606
607 fprintf(stderr, "x%i=[", Nside);
608 for (i=0; i<N; i++) {
609 for (j=0; j<N; j++) {
610 fprintf(stderr, "%i ", radectohealpix(i*2*M_PI/N, M_PI*(j-N/2)/N, Nside));
611 }
612 fprintf(stderr, ";");
613 }
614 fprintf(stderr, "];\n\n");
615 fflush(stderr);
616 }
617
test_big_nside(CuTest * ct)618 void test_big_nside(CuTest* ct) {
619 double ra1, dec1, ra2, dec2;
620 int Nside = 2097152;
621 int64_t hp;
622 double dx, dy;
623
624 // just a random val...
625 ra1 = 43.7;
626 dec1 = -38.4;
627
628 hp = radecdegtohealpixl(ra1, dec1, Nside);
629 healpixl_to_radecdeg(hp, Nside, 0, 0, &ra2, &dec2);
630
631 CuAssertDblEquals(ct, ra1, ra2, arcsec2deg(0.1));
632 CuAssertDblEquals(ct, dec1, dec2, arcsec2deg(0.1));
633
634 // another random val...
635 ra1 = 0.0003;
636 dec1 = 75.3;
637
638 dx = dy = -1.0;
639 hp = radecdegtohealpixlf(ra1, dec1, Nside, &dx, &dy);
640 CuAssert(ct, "dx", dx >= 0.0);
641 CuAssert(ct, "dx", dx <= 1.0);
642 CuAssert(ct, "dy", dy >= 0.0);
643 CuAssert(ct, "dy", dy <= 1.0);
644 healpixl_to_radecdeg(hp, Nside, dx, dy, &ra2, &dec2);
645
646 CuAssertDblEquals(ct, ra1, ra2, arcsec2deg(1e-10));
647 CuAssertDblEquals(ct, dec1, dec2, arcsec2deg(1e-10));
648
649 printf("RA,Dec difference: %g, %g arcsec\n", deg2arcsec(ra2-ra1), deg2arcsec(dec2-dec1));
650 }
651
test_distortion_at_pole(CuTest * ct)652 void test_distortion_at_pole(CuTest* ct) {
653 // not really a test of the code, more of healpix itself...
654 double ra1, dec1, ra2, dec2, ra3, dec3, ra4, dec4;
655 int Nside = 2097152;
656 int64_t hp;
657 double d1, d2, d3, d4, d5, d6;
658 double testras [] = { 0.0, 45.0, 0.0, 0.0 };
659 double testdecs[] = { 90.0, 50.0, 40.0, 0.0 };
660 char* testnames[] = { "north pole", "mid-polar", "mid-equatorial",
661 "equator" };
662 double ra, dec;
663 int i;
664
665 for (i=0; i<sizeof(testras)/sizeof(double); i++) {
666 ra = testras[i];
667 dec = testdecs[i];
668
669 hp = radecdegtohealpixl(ra, dec, Nside);
670
671 healpixl_to_radecdeg(hp, Nside, 0, 0, &ra1, &dec1);
672 healpixl_to_radecdeg(hp, Nside, 0, 1, &ra2, &dec2);
673 healpixl_to_radecdeg(hp, Nside, 1, 1, &ra3, &dec3);
674 healpixl_to_radecdeg(hp, Nside, 1, 0, &ra4, &dec4);
675
676 // sides
677 d1 = arcsec_between_radecdeg(ra1, dec1, ra2, dec2);
678 d2 = arcsec_between_radecdeg(ra2, dec2, ra3, dec3);
679 d3 = arcsec_between_radecdeg(ra3, dec3, ra4, dec4);
680 d4 = arcsec_between_radecdeg(ra4, dec4, ra1, dec1);
681 // diagonals
682 d5 = arcsec_between_radecdeg(ra1, dec1, ra3, dec3);
683 d6 = arcsec_between_radecdeg(ra2, dec2, ra4, dec4);
684
685 printf("%-15s (%4.1f, %4.1f): %-5.3f, %-5.3f, %-5.3f, %-5.3f / %-5.3f, %-5.3f\n", testnames[i], ra, dec, d1, d2, d3, d4, d5, d6);
686 }
687 }
688
689
690 #if defined(TEST_HEALPIX_MAIN)
main(int argc,char ** args)691 int main(int argc, char** args) {
692
693 /* Run all tests */
694 CuString *output = CuStringNew();
695 CuSuite* suite = CuSuiteNew();
696
697 /* Add new tests here */
698 SUITE_ADD_TEST(suite, test_healpix_neighbours);
699 SUITE_ADD_TEST(suite, test_healpix_pnprime_to_xy);
700 SUITE_ADD_TEST(suite, test_healpix_xy_to_pnprime);
701 SUITE_ADD_TEST(suite, test_healpix_distance_to_radec);
702
703 /* Run the suite, collect results and display */
704 CuSuiteRun(suite);
705 CuSuiteSummary(suite, output);
706 CuSuiteDetails(suite, output);
707 printf("%s\n", output->buffer);
708
709 /*
710 print_healpix_grid(1);
711 print_healpix_grid(2);
712 print_healpix_grid(3);
713 print_healpix_grid(4);
714 print_healpix_grid(5);
715 */
716
717 //print_test_healpix_output();
718
719 /*
720 int rastep, decstep;
721 int Nra = 100;
722 int Ndec = 100;
723 double ra, dec;
724 int healpix;
725 printf("radechealpix=zeros(%i,3);\n", Nra*Ndec);
726 for (rastep=0; rastep<Nra; rastep++) {
727 ra = ((double)rastep / (double)(Nra-1)) * 2.0 * M_PI;
728 for (decstep=0; decstep<Ndec; decstep++) {
729 dec = (((double)decstep / (double)(Ndec-1)) * M_PI) - M_PI/2.0;
730 healpix = radectohealpix(ra, dec);
731 printf("radechealpix(%i,:)=[%g,%g,%i];\n",
732 (rastep*Ndec) + decstep + 1, ra, dec, healpix);
733 }
734 }
735 printf("ra=radechealpix(:,1);\n");
736 printf("dec=radechealpix(:,2);\n");
737 printf("healpix=radechealpix(:,3);\n");
738 return 0;
739 */
740 return 0;
741 }
742 #endif
743