29#include "hipsmanager.h"
31#include "kstarsdata.h"
32#include "geolocation.h"
34static const double twothird = 2.0 / 3.0;
35static const double twopi = 6.283185307179586476925286766559005768394;
36static const double inv_halfpi = 0.6366197723675813430755350534900574;
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 };
41static 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)
52static const short utab[] =
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)
63static short xoffset[] = { -1, -1, 0, 1, 1, 1, 0, -1 };
64static short yoffset[] = { 0, 1, 1, 1, 0, -1, -1, -1 };
66static short facearray[9][12] =
68 { 8, 9, 10, 11, -1, -1, -1, -1, 10, 11, 8, 9 },
69 { 5, 6, 7, 4, 8, 9, 10, 11, 9, 10, 11, 8 },
70 { -1, -1, -1, -1, 5, 6, 7, 4, -1, -1, -1, -1 },
71 { 4, 5, 6, 7, 11, 8, 9, 10, 11, 8, 9, 10 },
72 { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 },
73 { 1, 2, 3, 0, 0, 1, 2, 3, 5, 6, 7, 4 },
74 { -1, -1, -1, -1, 7, 4, 5, 6, -1, -1, -1, -1 },
75 { 3, 0, 1, 2, 3, 0, 1, 2, 4, 5, 6, 7 },
76 { 2, 3, 0, 1, -1, -1, -1, -1, 0, 1, 2, 3 }
79static uchar swaparray[9][12] =
92static double fmodulo(
double v1,
double v2)
96 return (v1 < v2) ? v1 : fmod(v1, v2);
99 double tmp = fmod(v1, v2) + v2;
100 return (tmp == v2) ? 0. : tmp;
104void HEALPix::getCornerPoints(
int level,
int pix,
SkyPoint *skyCoords)
110 int nside = 1 <<
level;
111 boundaries(nside, pix, 1, v);
114 for (
int i = 0; i < 4; i++)
116 transformed[i].
setX(v[i].x());
117 transformed[i].
setY(v[i].y());
118 transformed[i].
setZ(v[i].z());
120 double ra = 0, de = 0;
121 xyz2sph(transformed[i], ra, de);
125 if (HIPSManager::Instance()->getCurrentFrame() == HIPSManager::HIPS_GALACTIC_FRAME)
127 dms galacticLong(ra);
131 skyCoords[i].
setRA0(skyCoords[i].ra());
136 skyCoords[i].
setRA0(ra / 15.0);
140 skyCoords[i].
updateCoords(KStarsData::Instance()->updateNum(),
false);
146void HEALPix::boundaries(qint32 nside, qint32 pix,
int step,
QVector3D *out)
150 nest2xyf(nside, pix, &ix, &iy, &fn);
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);
157 for (
int i = 0; i < step; ++i)
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);
166QVector3D HEALPix::toVec3(
double fx,
double fy,
int face)
168 double jr = jrll[face] - fx - fy;
173 bool lochave_sth =
false;
179 double tmp = nr * nr / 3.;
183 locsth = sqrt(tmp * (2. - tmp));
190 double tmp = nr * nr / 3.;
194 locsth = sqrt(tmp * (2. - tmp));
201 locz = (2 - jr) * 2. / 3.;
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;
209 double st = lochave_sth ? locsth : sqrt((1.0 - locz) * (1.0 + locz));
213 out.
setX(st * cos(locphi));
214 out.
setY(st * sin(locphi));
220void HEALPix::nest2xyf(
int nside,
int pix,
int *ix,
int *iy,
int *face_num)
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);
228 raw = (pix & 0x5555) | ((pix & 0x55550000) >> 15);
229 *iy = ctab[raw & 0xff] | (ctab[raw >> 8] << 4);
232static int64_t spread_bits64 (
int v)
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);
241static int xyf2nest (
int nside,
int ix,
int iy,
int face_num)
243 return (face_num * nside * nside) +
244 (utab[ix & 0xff] | (utab[ix >> 8] << 16)
245 | (utab[iy & 0xff] << 1) | (utab[iy >> 8] << 17));
270static int nside2order(qint32 nside)
274 while((nside >> (++i)) > 0);
291void HEALPix::neighbours(
int nside, qint32 ipix,
int *result)
293 int ix, iy, face_num;
295 int order = nside2order(nside);
297 nest2xyf(nside, ipix, &ix, &iy, &face_num);
299 qint32 nsm1 = nside - 1;
300 if ((ix > 0) && (ix < nsm1) && (iy > 0) && (iy < nsm1))
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;
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;
318 for (
int i = 0; i < 8; ++i)
320 int x = ix + xoffset[i];
321 int y = iy + yoffset[i];
344 int f = facearray[nbnum][face_num];
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);
357 result[i] = xyf2nest(nside, x, y, f);
365int HEALPix::ang2pix_nest_z_phi (qint32 nside_,
double z,
double phi)
368 double tt = fmodulo(phi, twopi) * inv_halfpi;
369 int face_num, ix, iy;
373 double temp1 = nside_ * (0.5 + tt);
374 double temp2 = nside_ * (z * 0.75);
375 int jp = (int)(temp1 - temp2);
376 int jm = (int)(temp1 + temp2);
377 int ifp = jp / nside_;
378 int ifm = jm / nside_;
379 face_num = (ifp == ifm) ? (ifp | 4) : ((ifp < ifm) ? ifp : (ifm + 8));
381 ix = jm & (nside_ - 1);
382 iy = nside_ - (jp & (nside_ - 1)) - 1;
386 int ntt = (int)tt, jp, jm;
388 if (ntt >= 4) ntt = 3;
390 tmp = nside_ * sqrt(3 * (1 - za));
392 jp = (int)(tp * tmp);
393 jm = (int)((1.0 - tp) * tmp);
394 if (jp >= nside_) jp = nside_ - 1;
395 if (jm >= nside_) jm = nside_ - 1;
399 ix = nside_ - jm - 1;
400 iy = nside_ - jp - 1;
410 return xyf2nest(nside_, ix, iy, face_num);
413int HEALPix::getPix(
int level,
double ra,
double dec)
415 int nside = 1 <<
level;
416 double polar[2] = { 0, 0 };
418 if (HIPSManager::Instance()->getCurrentFrame() == HIPSManager::HIPS_GALACTIC_FRAME)
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,
425 double rcb = cos(dec);
428 xyz = gl.mapVector(xyz);
429 xyz2sph(xyz, polar[1], polar[0]);
437 return ang2pix_nest_z_phi(nside, sin(polar[0]), polar[1]);
440void HEALPix::getPixChilds(
int pix,
int *childs)
442 childs[0] = pix * 4 + 0;
443 childs[1] = pix * 4 + 1;
444 childs[2] = pix * 4 + 2;
445 childs[3] = pix * 4 + 3;
448void HEALPix::xyz2sph(
const QVector3D &vec,
double &l,
double &b)
450 double rho = vec.
x() * vec.
x() + vec.
y() * vec.
y();
454 l = atan2(vec.
y(), vec.
x());
455 l -= floor(l / (M_PI * 2)) * (M_PI * 2);
456 b = atan2(vec.
z(), sqrt(rho));
467 b = (vec.
z() > 0.0) ? M_PI / 2. : -
dms::PI / 2.;
The sky coordinates of a point in the sky.
void B1950ToJ2000(void)
Exact precession from Besselian epoch 1950 to epoch J2000.
void EquatorialToHorizontal(const CachingDms *LST, const CachingDms *lat)
Determine the (Altitude, Azimuth) coordinates of the SkyPoint from its (RA, Dec) coordinates,...
void setRA0(dms r)
Sets RA0, the catalog Right Ascension.
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),...
void GalacticToEquatorial1950(const dms *galLong, const dms *galLat)
Computes equatorial coordinates referred to 1950 from galactic ones referred to epoch B1950.
void setDec0(dms d)
Sets Dec0, the catalog Declination.
An angle, stored as degrees, but expressible in many ways.
static constexpr double PI
PI is a const static member; it's public so that it can be used anywhere, as long as dms....
static constexpr double DegToRad
DegToRad is a const static member equal to the number of radians in one degree (dms::PI/180....
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)