← Back to branch summary

~guillaume-chereau/stellarium/hips

« back to all changes in this revision

Viewing changes to src/core/healpix.c

  • Committer: Guillaume Chereau
  • Date: 2017-08-30 08:46:08 UTC
  • Revision ID: guillaume@noctua-software.com-20170830084608-ldyvkulw4vzxplyh
First draft of HiPS surveys support

This is based on the Toast survey support.  In fact I replaced the dss
button to use HiPS instead of toast, since I think eventually we should
only use HiPS instead of toast.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/core/healpix.c
 
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
/* utab[m] = (short)(
 
23
      (m&0x1 )       | ((m&0x2 ) << 1) | ((m&0x4 ) << 2) | ((m&0x8 ) << 3)
 
24
    | ((m&0x10) << 4) | ((m&0x20) << 5) | ((m&0x40) << 6) | ((m&0x80) << 7)); */
 
25
static const short utab[]={
 
26
#define Z(a) 0x##a##0, 0x##a##1, 0x##a##4, 0x##a##5
 
27
#define Y(a) Z(a##0), Z(a##1), Z(a##4), Z(a##5)
 
28
#define X(a) Y(a##0), Y(a##1), Y(a##4), Y(a##5)
 
29
X(0),X(1),X(4),X(5)
 
30
#undef X
 
31
#undef Y
 
32
#undef Z
 
33
};
 
34
 
 
35
/* ctab[m] = (short)(
 
36
       (m&0x1 )       | ((m&0x2 ) << 7) | ((m&0x4 ) >> 1) | ((m&0x8 ) << 6)
 
37
    | ((m&0x10) >> 2) | ((m&0x20) << 5) | ((m&0x40) >> 3) | ((m&0x80) << 4)); */
 
38
static const short ctab[]={
 
39
#define Z(a) a,a+1,a+256,a+257
 
40
#define Y(a) Z(a),Z(a+2),Z(a+512),Z(a+514)
 
41
#define X(a) Y(a),Y(a+4),Y(a+1024),Y(a+1028)
 
42
X(0),X(8),X(2048),X(2056)
 
43
#undef X
 
44
#undef Y
 
45
#undef Z
 
46
};
 
47
 
 
48
// Position of the healpix faces.
 
49
static const int FACES[12][2] = {{1,  0}, {3,  0}, {5,  0}, {7,  0},
 
50
                                 {0, -1}, {2, -1}, {4, -1}, {6, -1},
 
51
                                 {1, -2}, {3, -2}, {5, -2}, {7, -2}};
 
52
 
 
53
int healpix_xyf2nest(int nside, int ix, int iy, int face_num)
 
54
{
 
55
  return (face_num*nside*nside) +
 
56
      (utab[ix&0xff] | (utab[ix>>8]<<16)
 
57
    | (utab[iy&0xff]<<1) | (utab[iy>>8]<<17));
 
58
}
 
59
 
 
60
static void nest2xyf(int nside, int pix, int *ix, int *iy, int *face_num)
 
61
{
 
62
    int npface_ = nside * nside, raw;
 
63
    *face_num = pix / npface_;
 
64
    pix &= (npface_ - 1);
 
65
    raw = (pix & 0x5555) | ((pix & 0x55550000) >> 15);
 
66
    *ix = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
 
67
    pix >>= 1;
 
68
    raw = (pix & 0x5555) | ((pix & 0x55550000) >> 15);
 
69
    *iy = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
 
70
}
 
71
 
 
72
// Create a 3x3 mat that transforms uv texture position to healpix xy
 
73
// coordinates.
 
74
void healpix_get_mat3(int nside, int pix, double out[3][3])
 
75
{
 
76
    int ix, iy, face;
 
77
    nest2xyf(nside, pix, &ix, &iy, &face);
 
78
    out[0][0] = +M_PI / 4 / nside;
 
79
    out[0][1] = +M_PI / 4 / nside;
 
80
    out[0][2] = 0;
 
81
    out[1][0] = -M_PI / 4 / nside;
 
82
    out[1][1] = +M_PI / 4 / nside;
 
83
    out[1][2] = 0;
 
84
    out[2][0] = (FACES[face][0] + (ix - iy + 0.0) / nside) * M_PI / 4;
 
85
    out[2][1] = (FACES[face][1] + (ix + iy + 0.0) / nside) * M_PI / 4;
 
86
    out[2][2] = 1;
 
87
}
 
88
 
 
89
static void healpix_xy2_z_phi(const double xy[2], double *z, double *phi)
 
90
{
 
91
    double x = xy[0], y = xy[1];
 
92
    double sigma, xc;
 
93
 
 
94
    if (fabs(y) > M_PI / 4) {
 
95
        // Polar
 
96
        sigma = 2 - fabs(y * 4) / M_PI;
 
97
        *z = (y > 0 ? 1 : -1) * (1 - sigma * sigma / 3);
 
98
        xc = -M_PI + (2 * floor((x + M_PI) * 4 / (2 * M_PI)) + 1) * M_PI / 4;
 
99
        *phi = sigma ? (xc + (x - xc) / sigma) : x;
 
100
    } else {
 
101
        // Equatorial
 
102
        *phi = x;
 
103
        *z = y * 8 / (M_PI * 3);
 
104
    }
 
105
}
 
106
 
 
107
void healpix_xy2ang(const double xy[2], double *theta, double *phi)
 
108
{
 
109
    double z;
 
110
    healpix_xy2_z_phi(xy, &z, phi);
 
111
    *theta = acos(z);
 
112
}
 
113
 
 
114
void healpix_xy2vec(const double xy[2], double out[3])
 
115
{
 
116
    double z, phi, stheta;
 
117
    healpix_xy2_z_phi(xy, &z, &phi);
 
118
    stheta = sqrt((1 - z) * (1 + z));
 
119
    out[0] = stheta * cos(phi);
 
120
    out[1] = stheta * sin(phi);
 
121
    out[2] = z;
 
122
}
 
123
 
 
124
void healpix_pix2vec(int nside, int pix, double out[3])
 
125
{
 
126
    int ix, iy, face;
 
127
    double xy[2];
 
128
    nest2xyf(nside, pix, &ix, &iy, &face);
 
129
    xy[0] = (FACES[face][0] + (ix - iy + 0.0) / nside) * M_PI / 4;
 
130
    xy[1] = (FACES[face][1] + (ix + iy + 1.0) / nside) * M_PI / 4;
 
131
    healpix_xy2vec(xy, out);
 
132
}
 
133
 
 
134
void healpix_pix2ang(int nside, int pix, double *theta, double *phi)
 
135
{
 
136
    int ix, iy, face;
 
137
    double xy[2];
 
138
    nest2xyf(nside, pix, &ix, &iy, &face);
 
139
    xy[0] = (FACES[face][0] + (ix - iy + 0.0) / nside) * M_PI / 4;
 
140
    xy[1] = (FACES[face][1] + (ix + iy + 1.0) / nside) * M_PI / 4;
 
141
    healpix_xy2ang(xy, theta, phi);
 
142
}