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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
|
#include <math.h>
/* HEALPix functions. The algorithm to subdivide the sky in equal area
tiles is described by Gorsky et al. 2005
(http://adsabs.harvard.edu/abs/2005ApJ...622..759G). All that is implemented
here is the conversion between RA/dec and HEALPix cell numbers, which is all
that's needed to figure out which lines in the bias data should be read. */
#define ASIN_TWO_THIRDS 0.7297276562269663634547966598133206953965059140477136907089494914618189
#define pi 3.1415926535897932384626433832795028841971693993751058209749445923
void ra_dec_to_xy( const double ra, const double dec, double *x, double *y);
unsigned xy_to_healpix( const double x, const double y, const unsigned N);
/* The HEALPix projection uses an equal-area cylindrical projection for the
central 2/3 (by area) zone nearest the equator, and puts each polar region
into four right triangles in an interrupted Collignon projection. See
https://en.wikipedia.org/wiki/HEALPix for a picture of this worth many thousands
of comments. Here, y=0 at the equator; y=2 at the north pole, y=-2 at the
south pole, and 0 to 360 in RA is projected to 0 to 8 in x. RA and dec are
assumed to be in radians. */
void ra_dec_to_xy( const double ra, const double dec, double *x, double *y)
{
if( dec < -ASIN_TWO_THIRDS)
{
ra_dec_to_xy( ra, -dec, x, y);
*y = -*y;
return;
}
*x = fmod( ra * 4 / pi, 8.);
if( *x < 0.)
*x += 8.;
if( dec < ASIN_TWO_THIRDS)
*y = 1.5 * sin( dec);
else
{
const double z = sqrt(3. * (1. - sin( dec)));
const double center = floor( *x / 2.) * 2. + 1.;
*y = 2. - z;
*x = center + z * (*x - center);
}
}
#ifdef NOT_NEEDED_YET
/* So far, we aren't doing reverse conversions. Note that this is, */
/* as yet, untested code. Though it's a straightforward inverse of */
/* the above ra_dec_to_xy function. */
static void xy_to_ra_dec( const double x, const double y, double *ra, double *dec)
{
if( y < -1.)
{
xy_to_ra_dec( x, -y, ra, dec);
*dec = -*dec;
return;
}
if( y < 1.)
{
*ra = x * pi / 4.;
*dec = asin( y * 2. / 3.);
}
else
{
const double z = 2. - y;
const double center = floor( (x + 1.) / 2.) * 2.;
*ra = (center + (x - center) / z) * pi / 4.;
*dec = asin( 1. - z * z / 3.);
}
}
#endif
/* HEALPix pixels are squares (when projected) tilted at a 45-degree angle.
It helps to rotate them as shown below, with the origin shifted to the top
of the leftmost polar triangle. (In the following, to show the 192 cells of
N=2 with two digits, I mix hexadecimal and decimal; the cell numbers run
from ...98, 99, to a0, a1... and similarly up to the maximum j1 = cell 191.)
03 11 23 39
10 22 38 55 71
21 37 54 70 87 a3
36 53 69 86 a2 b9 d5
02 09 20 35 52 68 85 a1 b8 d4 f1 g7
i=0 1 2 3 4 5 6 08 19 34 51 67 84 a0 b7 d3 f0 g6 h9
18 33 50 66 83 99 b6 d2 e9 g5 h8 i7
32 49 65 82 98 b5 d1 e8 g4 h7 i6 j1
01 07 17 31 48 64 81 97 b4 d0 e7 g3
06 16 30 47 63 80 96 b3 c9 e6 g2 h6
15 29 46 62 79 95 b2 c8 e5 g1 h5 i5 -2 (i, j) rotated to keep cell
28 45 61 78 94 b1 c7 e4 g0 h4 i4 j0 -1 00 at the origin
00 05 14 27 44 60 77 93 b0 c6 e3 f9 j=0
04 13 26 43 59 76 92 a9 c5 e2 f8 h3 1
12 25 42 58 75 91 a8 c4 e1 f7 h2 i3 2 N = 4 = power of 2
24 41 57 74 90 a7 c3 e0 f6 h1 i2 i9 3
40 56 73 89 a6 c2 d9 f5 4
72 88 a5 c1 d8 f4 h0
a4 c0 d7 f3 g9 i1
d6 f2 g8 i0 i8
In this scheme, lines of constant latitude would run from lower left
to upper right; meridians would run (mostly) from lower right to upper
left, with some distortion in the polar zones.
NOTE: in unrotated form, tile 00 = north pole = (x=1, y=2); tile 24 =
(x=0, y=1); tile 96 = center = (x=2,y=0). Note north-south symmetry; if
(x, y) lands in tile n, (8-x, y) lands in tile n_tiles - 1 - n, where
n_tiles = 12N^2. Which is why the code handles only the northern polar
triangles and the equatorial zone; the southern triangles are handled
just by rotating the coordinates into the northern triangles. */
unsigned xy_to_healpix( const double x, const double y, const unsigned N)
{
int i, j, line;
unsigned rval = 0;
if( y <= -1.) /* use rotational symmetry for south triangles */
return( 12 * N * N - 1 - xy_to_healpix( 8. - x, -y, N));
// x -= 1.; /* translate to top of first polar triangle, tile 00 */
// y -= 2.;
i = (int)floor( (double)N * (x - y + 1.) / 2.);
j = (int)floor( -(double)N * (x + y - 3.) / 2.);
line = i + j;
if( (unsigned)line < N) /* north polar triangles */
rval = line * (line + 1) * 2 + i % N + (i / N) * (line + 1);
else
{ /* equatorial region */
const unsigned offset_along_line =
(unsigned)( i - (line + 1 - N) / 2) % (4 * N);
rval = N * (1 - N) * 2 + 4 * line * N + offset_along_line;
}
return( rval);
}
#include <stdio.h>
#include <stdlib.h>
#ifdef OLD_TEST_MAIN
int main( const int argc, const char **argv)
{
const double ra = atof( argv[1]) * pi / 180;
const double dec = atof( argv[2]) * pi / 180;
double x, y;
unsigned hp;
ra_dec_to_xy( ra, dec, &x, &y);
printf( "x = %f y = %f\n", x, y);
hp = xy_to_healpix( x, y, 4);
printf( "hp = %u\n", hp);
return( 0);
}
#endif /* #ifdef OLD_TEST_MAIN */
#ifdef OTHER_TEST_MAIN
/* Uses a list of HEALPix tile centers provided with the astrometric */
/* bias data, and checks that the above code gives the same tile */
/* numbers as the file does. "Noise" can be added to jostle the */
/* positions around within the tile. */
int main( const int argc, const char **argv)
{
FILE *ifile = fopen( "tiles.dat", "rb");
char buff[100];
const double noise = (argc > 1 ? atof( argv[1]) : 0) * pi / 180.;
while( fgets( buff, sizeof( buff), ifile))
if( *buff != '!')
{
unsigned tile_no, computed;
double ra, dec, x, y;
sscanf( buff, "%u %lf %lf", &tile_no, &ra, &dec);
if( noise)
{
ra += ((2. * (double)rand( ) / (double)RAND_MAX) - 1.) * noise
/ cos( dec);
dec += ((2. * (double)rand( ) / (double)RAND_MAX) - 1.) * noise;
}
ra_dec_to_xy( ra, dec, &x, &y);
computed = xy_to_healpix( x, y, 64);
if( computed != tile_no)
printf( "Computed %u\n%s", computed, buff);
}
fclose( ifile);
return( 0);
}
#endif /* #ifdef OTHER_TEST_MAIN */
|
