1 /*
2 * Stellarium
3 * Copyright (C) 2017 Guillaume Chereau
4 *
5 * This program is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU General Public License
7 * as published by the Free Software Foundation; either version 2
8 * of the License, or (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
18 */
19
20 #include <math.h>
21
22 #ifndef M_PI
23 #define M_PI 3.14159265358979323846264338327950288
24 #endif
25
26 /* utab[m] = (short)(
27 (m&0x1 ) | ((m&0x2 ) << 1) | ((m&0x4 ) << 2) | ((m&0x8 ) << 3)
28 | ((m&0x10) << 4) | ((m&0x20) << 5) | ((m&0x40) << 6) | ((m&0x80) << 7)); */
29 static const short utab[]={
30 #define Z(a) 0x##a##0, 0x##a##1, 0x##a##4, 0x##a##5
31 #define Y(a) Z(a##0), Z(a##1), Z(a##4), Z(a##5)
32 #define X(a) Y(a##0), Y(a##1), Y(a##4), Y(a##5)
33 X(0),X(1),X(4),X(5)
34 #undef X
35 #undef Y
36 #undef Z
37 };
38
39 /* ctab[m] = (short)(
40 (m&0x1 ) | ((m&0x2 ) << 7) | ((m&0x4 ) >> 1) | ((m&0x8 ) << 6)
41 | ((m&0x10) >> 2) | ((m&0x20) << 5) | ((m&0x40) >> 3) | ((m&0x80) << 4)); */
42 static const short ctab[]={
43 #define Z(a) a,a+1,a+256,a+257
44 #define Y(a) Z(a),Z(a+2),Z(a+512),Z(a+514)
45 #define X(a) Y(a),Y(a+4),Y(a+1024),Y(a+1028)
46 X(0),X(8),X(2048),X(2056)
47 #undef X
48 #undef Y
49 #undef Z
50 };
51
52 // Position of the healpix faces.
53 static const int FACES[12][2] = {{1, 0}, {3, 0}, {5, 0}, {7, 0},
54 {0, -1}, {2, -1}, {4, -1}, {6, -1},
55 {1, -2}, {3, -2}, {5, -2}, {7, -2}};
56
healpix_xyf2nest(int nside,int ix,int iy,int face_num)57 int healpix_xyf2nest(int nside, int ix, int iy, int face_num)
58 {
59 return (face_num*nside*nside) +
60 (utab[ix&0xff] | (utab[ix>>8]<<16)
61 | (utab[iy&0xff]<<1) | (utab[iy>>8]<<17));
62 }
63
nest2xyf(int nside,int pix,int * ix,int * iy,int * face_num)64 static void nest2xyf(int nside, int pix, int *ix, int *iy, int *face_num)
65 {
66 int npface_ = nside * nside, raw;
67 *face_num = pix / npface_;
68 pix &= (npface_ - 1);
69 raw = (pix & 0x5555) | ((pix & 0x55550000) >> 15);
70 *ix = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
71 pix >>= 1;
72 raw = (pix & 0x5555) | ((pix & 0x55550000) >> 15);
73 *iy = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
74 }
75
76 // Create a 3x3 mat that transforms uv texture position to healpix xy
77 // coordinates.
healpix_get_mat3(int nside,int pix,double out[3][3])78 void healpix_get_mat3(int nside, int pix, double out[3][3])
79 {
80 int ix, iy, face;
81 nest2xyf(nside, pix, &ix, &iy, &face);
82 out[0][0] = +M_PI / 4 / nside;
83 out[0][1] = +M_PI / 4 / nside;
84 out[0][2] = 0;
85 out[1][0] = -M_PI / 4 / nside;
86 out[1][1] = +M_PI / 4 / nside;
87 out[1][2] = 0;
88 out[2][0] = (FACES[face][0] + (ix - iy + 0.0) / nside) * M_PI / 4;
89 out[2][1] = (FACES[face][1] + (ix + iy + 0.0) / nside) * M_PI / 4;
90 out[2][2] = 1;
91 }
92
healpix_xy2_z_phi(const double xy[2],double * z,double * phi)93 static void healpix_xy2_z_phi(const double xy[2], double *z, double *phi)
94 {
95 double x = xy[0], y = xy[1];
96 double sigma, xc;
97
98 if (fabs(y) > M_PI / 4) {
99 // Polar
100 sigma = 2 - fabs(y * 4) / M_PI;
101 *z = (y > 0 ? 1 : -1) * (1 - sigma * sigma / 3);
102 xc = -M_PI + (2 * floor((x + M_PI) * 4 / (2 * M_PI)) + 1) * M_PI / 4;
103 *phi = sigma ? (xc + (x - xc) / sigma) : x;
104 } else {
105 // Equatorial
106 *phi = x;
107 *z = y * 8 / (M_PI * 3);
108 }
109 }
110
healpix_xy2ang(const double xy[2],double * theta,double * phi)111 void healpix_xy2ang(const double xy[2], double *theta, double *phi)
112 {
113 double z;
114 healpix_xy2_z_phi(xy, &z, phi);
115 *theta = acos(z);
116 }
117
healpix_xy2vec(const double xy[2],double out[3])118 void healpix_xy2vec(const double xy[2], double out[3])
119 {
120 double z, phi, stheta;
121 healpix_xy2_z_phi(xy, &z, &phi);
122 stheta = sqrt((1 - z) * (1 + z));
123 out[0] = stheta * cos(phi);
124 out[1] = stheta * sin(phi);
125 out[2] = z;
126 }
127
healpix_pix2vec(int nside,int pix,double out[3])128 void healpix_pix2vec(int nside, int pix, double out[3])
129 {
130 int ix, iy, face;
131 double xy[2];
132 nest2xyf(nside, pix, &ix, &iy, &face);
133 xy[0] = (FACES[face][0] + (ix - iy + 0.0) / nside) * M_PI / 4;
134 xy[1] = (FACES[face][1] + (ix + iy + 1.0) / nside) * M_PI / 4;
135 healpix_xy2vec(xy, out);
136 }
137
healpix_pix2ang(int nside,int pix,double * theta,double * phi)138 void healpix_pix2ang(int nside, int pix, double *theta, double *phi)
139 {
140 int ix, iy, face;
141 double xy[2];
142 nest2xyf(nside, pix, &ix, &iy, &face);
143 xy[0] = (FACES[face][0] + (ix - iy + 0.0) / nside) * M_PI / 4;
144 xy[1] = (FACES[face][1] + (ix + iy + 1.0) / nside) * M_PI / 4;
145 healpix_xy2ang(xy, theta, phi);
146 }
147