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 <[email protected]>
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 
34 static const double twothird = 2.0 / 3.0;
35 static const double twopi = 6.283185307179586476925286766559005768394;
36 static const double inv_halfpi = 0.6366197723675813430755350534900574;
37 
38 static const int jrll[] = { 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4 };
39 static const int jpll[] = { 1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7 };
40 
41 static 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 
52 static 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 
63 static short xoffset[] = { -1, -1, 0, 1, 1, 1, 0, -1 };
64 static short yoffset[] = { 0, 1, 1, 1, 0, -1, -1, -1 };
65 
66 static 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 
79 static 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 
92 static 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 
104 void 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 
146 void 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 
166 QVector3D 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 
220 void 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 
232 static 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 
241 static 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 /*
250 int 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 /*
264 static int ilog2(qint32 arg)
265 {
266  return 32 - leadingZeros(qMax(arg, 1));
267 }
268 */
269 
270 static 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 
291 void 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 
365 int 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 
413 int 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 
440 void 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 
448 void 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 
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:380
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
Stores dms coordinates for a point in the sky. for converting between coordinate systems.
Definition: skypoint.h:44
void setDec0(dms d)
Sets Dec0, the catalog Declination.
Definition: skypoint.h:119
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:385
QStringView level(QStringView ifopt)
void setRA0(dms r)
Sets RA0, the catalog Right Ascension.
Definition: skypoint.h:94
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
KGUIADDONS_EXPORT QColor tint(const QColor &base, const QColor &color, qreal amount=0.3)
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 B1950ToJ2000(void)
Exact precession from Besselian epoch 1950 to epoch J2000.
Definition: skypoint.cpp:773
QTextStream & dec(QTextStream &stream)
GeoCoordinates geo(const QVariant &location)
An angle, stored as degrees, but expressible in many ways.
Definition: dms.h:37
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-2022 The KDE developers.
Generated on Sat Aug 13 2022 04:01:54 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.