Kstars

healpix.cpp
1/* -----------------------------------------------------------------------------
2
3 SPDX-FileCopyrightText: 1997-2016 Martin Reinecke,
4 SPDX-FileCopyrightText: 1997-2016 Krzysztof M. Gorski
5 SPDX-FileCopyrightText: 1997-2016 Eric Hivon
6 SPDX-FileCopyrightText: Benjamin D. Wandelt
7 SPDX-FileCopyrightText: Anthony J. Banday,
8 SPDX-FileCopyrightText: Matthias Bartelmann,
9 SPDX-FileCopyrightText: Reza Ansari
10 SPDX-FileCopyrightText: Kenneth M. Ganga
11
12 This file is part of HEALPix.
13
14 Based on work by Pavel Mraz from SkyTechX.
15
16 Modified by Jasem Mutlaq for KStars:
17 SPDX-FileCopyrightText: Jasem Mutlaq <mutlaqja@ikarustech.com>
18
19 SPDX-License-Identifier: GPL-2.0-or-later
20
21 For more information about HEALPix see https://healpix.sourceforge.io/
22
23 ---------------------------------------------------------------------------*/
24
25#include "healpix.h"
26
27#include <QMatrix4x4>
28
29#include "hipsmanager.h"
30#include "skypoint.h"
31#include "kstarsdata.h"
32#include "geolocation.h"
33
34static const double twothird = 2.0 / 3.0;
35static const double twopi = 6.283185307179586476925286766559005768394;
36static const double inv_halfpi = 0.6366197723675813430755350534900574;
37
38static const int jrll[] = { 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4 };
39static const int jpll[] = { 1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7 };
40
41static const short ctab[] =
42{
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
52static const short utab[] =
53{
54#define Z(a) 0x##a##0, 0x##a##1, 0x##a##4, 0x##a##5
55#define Y(a) Z(a##0), Z(a##1), Z(a##4), Z(a##5)
56#define X(a) Y(a##0), Y(a##1), Y(a##4), Y(a##5)
57 X(0), X(1), X(4), X(5)
58#undef X
59#undef Y
60#undef Z
61};
62
63static short xoffset[] = { -1, -1, 0, 1, 1, 1, 0, -1 };
64static short yoffset[] = { 0, 1, 1, 1, 0, -1, -1, -1 };
65
66static short facearray[9][12] =
67{
68 { 8, 9, 10, 11, -1, -1, -1, -1, 10, 11, 8, 9 }, // S
69 { 5, 6, 7, 4, 8, 9, 10, 11, 9, 10, 11, 8 }, // SE
70 { -1, -1, -1, -1, 5, 6, 7, 4, -1, -1, -1, -1 }, // E
71 { 4, 5, 6, 7, 11, 8, 9, 10, 11, 8, 9, 10 }, // SW
72 { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 }, // center
73 { 1, 2, 3, 0, 0, 1, 2, 3, 5, 6, 7, 4 }, // NE
74 { -1, -1, -1, -1, 7, 4, 5, 6, -1, -1, -1, -1 }, // W
75 { 3, 0, 1, 2, 3, 0, 1, 2, 4, 5, 6, 7 }, // NW
76 { 2, 3, 0, 1, -1, -1, -1, -1, 0, 1, 2, 3 }
77}; // N
78
79static uchar swaparray[9][12] =
80{
81 { 0, 0, 3 }, // S
82 { 0, 0, 6 }, // SE
83 { 0, 0, 0 }, // E
84 { 0, 0, 5 }, // SW
85 { 0, 0, 0 }, // center
86 { 5, 0, 0 }, // NE
87 { 0, 0, 0 }, // W
88 { 6, 0, 0 }, // NW
89 { 3, 0, 0 }
90}; // N
91
92static double fmodulo(double v1, double v2)
93{
94 if (v1 >= 0)
95 {
96 return (v1 < v2) ? v1 : fmod(v1, v2);
97 }
98
99 double tmp = fmod(v1, v2) + v2;
100 return (tmp == v2) ? 0. : tmp;
101}
102
103
104void HEALPix::getCornerPoints(int level, int pix, SkyPoint *skyCoords)
105{
106 QVector3D v[4];
107 // Transform from HealPIX convention to KStars
108 QVector3D transformed[4];
109
110 int nside = 1 << level;
111 boundaries(nside, pix, 1, v);
112
113 // From rectangular coordinates to Sky coordinates
114 for (int i = 0; i < 4; i++)
115 {
116 transformed[i].setX(v[i].x());
117 transformed[i].setY(v[i].y());
118 transformed[i].setZ(v[i].z());
119
120 double ra = 0, de = 0;
121 xyz2sph(transformed[i], ra, de);
122 de /= dms::DegToRad;
123 ra /= dms::DegToRad;
124
125 if (HIPSManager::Instance()->getCurrentFrame() == HIPSManager::HIPS_GALACTIC_FRAME)
126 {
127 dms galacticLong(ra);
128 dms galacticLat(de);
129 skyCoords[i].GalacticToEquatorial1950(&galacticLong, &galacticLat);
130 skyCoords[i].B1950ToJ2000();
131 skyCoords[i].setRA0(skyCoords[i].ra());
132 skyCoords[i].setDec0(skyCoords[i].dec());
133 }
134 else
135 {
136 skyCoords[i].setRA0(ra / 15.0);
137 skyCoords[i].setDec0(de);
138 }
139
140 skyCoords[i].updateCoords(KStarsData::Instance()->updateNum(), false);
141 skyCoords[i].EquatorialToHorizontal(KStarsData::Instance()->lst(), KStarsData::Instance()->geo()->lat());
142 }
143
144}
145
146void HEALPix::boundaries(qint32 nside, qint32 pix, int step, QVector3D *out)
147{
148 int ix, iy, fn;
149
150 nest2xyf(nside, pix, &ix, &iy, &fn);
151
152 double dc = 0.5 / nside;
153 double xc = (ix + 0.5) / nside;
154 double yc = (iy + 0.5) / nside;
155 double d = 1. / (step * nside);
156
157 for (int i = 0; i < step; ++i)
158 {
159 out[0] = toVec3(xc + dc - i * d, yc + dc, fn);
160 out[1] = toVec3(xc - dc, yc + dc - i * d, fn);
161 out[2] = toVec3(xc - dc + i * d, yc - dc, fn);
162 out[3] = toVec3(xc + dc, yc - dc + i * d, fn);
163 }
164}
165
166QVector3D HEALPix::toVec3(double fx, double fy, int face)
167{
168 double jr = jrll[face] - fx - fy;
169
170 double locz;
171 double locsth;
172 double locphi;
173 bool lochave_sth = false;
174
175 double nr;
176 if (jr < 1)
177 {
178 nr = jr;
179 double tmp = nr * nr / 3.;
180 locz = 1 - tmp;
181 if (locz > 0.99)
182 {
183 locsth = sqrt(tmp * (2. - tmp));
184 lochave_sth = true;
185 }
186 }
187 else if (jr > 3)
188 {
189 nr = 4 - jr;
190 double tmp = nr * nr / 3.;
191 locz = tmp - 1;
192 if (locz < -0.99)
193 {
194 locsth = sqrt(tmp * (2. - tmp));
195 lochave_sth = true;
196 }
197 }
198 else
199 {
200 nr = 1;
201 locz = (2 - jr) * 2. / 3.;
202 }
203
204 double tmp = jpll[face] * nr + fx - fy;
205 if (tmp < 0) tmp += 8;
206 if (tmp >= 8) tmp -= 8;
207 locphi = (nr < 1e-15) ? 0 : (0.5 * dms::PI / 2.0 * tmp) / nr;
208
209 double st = lochave_sth ? locsth : sqrt((1.0 - locz) * (1.0 + locz));
210
211 QVector3D out;
212
213 out.setX(st * cos(locphi));
214 out.setY(st * sin(locphi));
215 out.setZ(locz);
216
217 return out;
218}
219
220void HEALPix::nest2xyf(int nside, int pix, int *ix, int *iy, int *face_num)
221{
222 int npface_ = nside * nside, raw;
223 *face_num = pix / npface_;
224 pix &= (npface_ - 1);
225 raw = (pix & 0x5555) | ((pix & 0x55550000) >> 15);
226 *ix = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
227 pix >>= 1;
228 raw = (pix & 0x5555) | ((pix & 0x55550000) >> 15);
229 *iy = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
230}
231
232static int64_t spread_bits64 (int v)
233{
234 return (int64_t)(utab[ v & 0xff])
235 | ((int64_t)(utab[(v >> 8) & 0xff]) << 16)
236 | ((int64_t)(utab[(v >> 16) & 0xff]) << 32)
237 | ((int64_t)(utab[(v >> 24) & 0xff]) << 48);
238}
239
240
241static int xyf2nest (int nside, int ix, int iy, int face_num)
242{
243 return (face_num * nside * nside) +
244 (utab[ix & 0xff] | (utab[ix >> 8] << 16)
245 | (utab[iy & 0xff] << 1) | (utab[iy >> 8] << 17));
246}
247
248
249/*
250int static leadingZeros(qint32 value)
251{
252 int leadingZeros = 0;
253 while(value != 0)
254 {
255 value = value >> 1;
256 leadingZeros++;
257 }
258
259 return (32 - leadingZeros);
260}
261*/
262
263/*
264static int ilog2(qint32 arg)
265{
266 return 32 - leadingZeros(qMax(arg, 1));
267}
268*/
269
270static int nside2order(qint32 nside)
271{
272 {
273 int i = 0;
274 while((nside >> (++i)) > 0);
275 return --i;
276 }
277}
278
279/** Returns the neighboring pixels of ipix.
280 This method works in both RING and NEST schemes, but is
281 considerably faster in the NEST scheme.
282 @param nside defines the size of the neighboring area
283 @param ipix the requested pixel number.
284 @param result the result array
285 @return array with indices of the neighboring pixels.
286 The returned array contains (in this order)
287 the pixel numbers of the SW, W, NW, N, NE, E, SE and S neighbor
288 of ipix. If a neighbor does not exist (this can only happen
289 for the W, N, E and S neighbors), its entry is set to -1. */
290
291void HEALPix::neighbours(int nside, qint32 ipix, int *result)
292{
293 int ix, iy, face_num;
294
295 int order = nside2order(nside);
296
297 nest2xyf(nside, ipix, &ix, &iy, &face_num);
298
299 qint32 nsm1 = nside - 1;
300 if ((ix > 0) && (ix < nsm1) && (iy > 0) && (iy < nsm1))
301 {
302 qint32 fpix = (qint32)(face_num) << (2 * order),
303 px0 = spread_bits64(ix ), py0 = spread_bits64(iy ) << 1,
304 pxp = spread_bits64(ix + 1), pyp = spread_bits64(iy + 1) << 1,
305 pxm = spread_bits64(ix - 1), pym = spread_bits64(iy - 1) << 1;
306
307 result[0] = fpix + pxm + py0;
308 result[1] = fpix + pxm + pyp;
309 result[2] = fpix + px0 + pyp;
310 result[3] = fpix + pxp + pyp;
311 result[4] = fpix + pxp + py0;
312 result[5] = fpix + pxp + pym;
313 result[6] = fpix + px0 + pym;
314 result[7] = fpix + pxm + pym;
315 }
316 else
317 {
318 for (int i = 0; i < 8; ++i)
319 {
320 int x = ix + xoffset[i];
321 int y = iy + yoffset[i];
322 int nbnum = 4;
323 if (x < 0)
324 {
325 x += nside;
326 nbnum -= 1;
327 }
328 else if (x >= nside)
329 {
330 x -= nside;
331 nbnum += 1;
332 }
333 if (y < 0)
334 {
335 y += nside;
336 nbnum -= 3;
337 }
338 else if (y >= nside)
339 {
340 y -= nside;
341 nbnum += 3;
342 }
343
344 int f = facearray[nbnum][face_num];
345
346 if (f >= 0)
347 {
348 int bits = swaparray[nbnum][face_num >> 2];
349 if ((bits & 1) > 0) x = (int)(nside - x - 1);
350 if ((bits & 2) > 0) y = (int)(nside - y - 1);
351 if ((bits & 4) > 0)
352 {
353 int tint = x;
354 x = y;
355 y = tint;
356 }
357 result[i] = xyf2nest(nside, x, y, f);
358 }
359 else
360 result[i] = -1;
361 }
362 }
363}
364
365int HEALPix::ang2pix_nest_z_phi (qint32 nside_, double z, double phi)
366{
367 double za = fabs(z);
368 double tt = fmodulo(phi, twopi) * inv_halfpi; /* in [0,4) */
369 int face_num, ix, iy;
370
371 if (za <= twothird) /* Equatorial region */
372 {
373 double temp1 = nside_ * (0.5 + tt);
374 double temp2 = nside_ * (z * 0.75);
375 int jp = (int)(temp1 - temp2); /* index of ascending edge line */
376 int jm = (int)(temp1 + temp2); /* index of descending edge line */
377 int ifp = jp / nside_; /* in {0,4} */
378 int ifm = jm / nside_;
379 face_num = (ifp == ifm) ? (ifp | 4) : ((ifp < ifm) ? ifp : (ifm + 8));
380
381 ix = jm & (nside_ - 1);
382 iy = nside_ - (jp & (nside_ - 1)) - 1;
383 }
384 else /* polar region, za > 2/3 */
385 {
386 int ntt = (int)tt, jp, jm;
387 double tp, tmp;
388 if (ntt >= 4) ntt = 3;
389 tp = tt - ntt;
390 tmp = nside_ * sqrt(3 * (1 - za));
391
392 jp = (int)(tp * tmp); /* increasing edge line index */
393 jm = (int)((1.0 - tp) * tmp); /* decreasing edge line index */
394 if (jp >= nside_) jp = nside_ - 1; /* for points too close to the boundary */
395 if (jm >= nside_) jm = nside_ - 1;
396 if (z >= 0)
397 {
398 face_num = ntt; /* in {0,3} */
399 ix = nside_ - jm - 1;
400 iy = nside_ - jp - 1;
401 }
402 else
403 {
404 face_num = ntt + 8; /* in {8,11} */
405 ix = jp;
406 iy = jm;
407 }
408 }
409
410 return xyf2nest(nside_, ix, iy, face_num);
411}
412
413int HEALPix::getPix(int level, double ra, double dec)
414{
415 int nside = 1 << level;
416 double polar[2] = { 0, 0 };
417
418 if (HIPSManager::Instance()->getCurrentFrame() == HIPSManager::HIPS_GALACTIC_FRAME)
419 {
420 static QMatrix4x4 gl(-0.0548762f, -0.873437f, -0.483835f, 0,
421 0.4941100f, -0.444830f, 0.746982f, 0,
422 -0.8676660f, -0.198076f, 0.455984f, 0,
423 0, 0, 0, 1);
424
425 double rcb = cos(dec);
426 QVector3D xyz = QVector3D(rcb * cos(ra), rcb * sin(ra), sin(dec));
427
428 xyz = gl.mapVector(xyz);
429 xyz2sph(xyz, polar[1], polar[0]);
430 }
431 else
432 {
433 polar[0] = dec;
434 polar[1] = ra;
435 }
436
437 return ang2pix_nest_z_phi(nside, sin(polar[0]), polar[1]);
438}
439
440void HEALPix::getPixChilds(int pix, int *childs)
441{
442 childs[0] = pix * 4 + 0;
443 childs[1] = pix * 4 + 1;
444 childs[2] = pix * 4 + 2;
445 childs[3] = pix * 4 + 3;
446}
447
448void HEALPix::xyz2sph(const QVector3D &vec, double &l, double &b)
449{
450 double rho = vec.x() * vec.x() + vec.y() * vec.y();
451
452 if (rho > 0)
453 {
454 l = atan2(vec.y(), vec.x());
455 l -= floor(l / (M_PI * 2)) * (M_PI * 2);
456 b = atan2(vec.z(), sqrt(rho));
457 }
458 else
459 {
460 l = 0.0;
461 if (vec.z() == 0.0)
462 {
463 b = 0.0;
464 }
465 else
466 {
467 b = (vec.z() > 0.0) ? M_PI / 2. : -dms::PI / 2.;
468 }
469 }
470}
471
472
473
The sky coordinates of a point in the sky.
Definition skypoint.h:45
void B1950ToJ2000(void)
Exact precession from Besselian epoch 1950 to epoch J2000.
Definition skypoint.cpp:773
void EquatorialToHorizontal(const CachingDms *LST, const CachingDms *lat)
Determine the (Altitude, Azimuth) coordinates of the SkyPoint from its (RA, Dec) coordinates,...
Definition skypoint.cpp:77
void setRA0(dms r)
Sets RA0, the catalog Right Ascension.
Definition skypoint.h:94
virtual void updateCoords(const KSNumbers *num, bool includePlanets=true, const CachingDms *lat=nullptr, const CachingDms *LST=nullptr, bool forceRecompute=false)
Determine the current coordinates (RA, Dec) from the catalog coordinates (RA0, Dec0),...
Definition skypoint.cpp:582
void GalacticToEquatorial1950(const dms *galLong, const dms *galLat)
Computes equatorial coordinates referred to 1950 from galactic ones referred to epoch B1950.
Definition skypoint.cpp:754
void setDec0(dms d)
Sets Dec0, the catalog Declination.
Definition skypoint.h:119
An angle, stored as degrees, but expressible in many ways.
Definition dms.h:38
static constexpr double PI
PI is a const static member; it's public so that it can be used anywhere, as long as dms....
Definition dms.h:385
static constexpr double DegToRad
DegToRad is a const static member equal to the number of radians in one degree (dms::PI/180....
Definition dms.h:390
KGUIADDONS_EXPORT QColor tint(const QColor &base, const QColor &color, qreal amount=0.3)
GeoCoordinates geo(const QVariant &location)
QStringView level(QStringView ifopt)
QTextStream & dec(QTextStream &stream)
void setX(float x)
void setY(float y)
void setZ(float z)
float x() const const
float y() const const
float z() const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Mon Nov 4 2024 16:38:43 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.