File indexing completed on 2024-03-24 03:46:43
0001 /* ----------------------------------------------------------------------------- 0002 0003 SPDX-FileCopyrightText: 1997-2016 Martin Reinecke, 0004 SPDX-FileCopyrightText: 1997-2016 Krzysztof M. Gorski 0005 SPDX-FileCopyrightText: 1997-2016 Eric Hivon 0006 SPDX-FileCopyrightText: Benjamin D. Wandelt 0007 SPDX-FileCopyrightText: Anthony J. Banday, 0008 SPDX-FileCopyrightText: Matthias Bartelmann, 0009 SPDX-FileCopyrightText: Reza Ansari 0010 SPDX-FileCopyrightText: Kenneth M. Ganga 0011 0012 This file is part of HEALPix. 0013 0014 Based on work by Pavel Mraz from SkyTechX. 0015 0016 Modified by Jasem Mutlaq for KStars: 0017 SPDX-FileCopyrightText: Jasem Mutlaq <mutlaqja@ikarustech.com> 0018 0019 SPDX-License-Identifier: GPL-2.0-or-later 0020 0021 For more information about HEALPix see https://healpix.sourceforge.io/ 0022 0023 ---------------------------------------------------------------------------*/ 0024 0025 #include "healpix.h" 0026 0027 #include <QMatrix4x4> 0028 0029 #include "hipsmanager.h" 0030 #include "skypoint.h" 0031 #include "kstarsdata.h" 0032 #include "geolocation.h" 0033 0034 static const double twothird = 2.0 / 3.0; 0035 static const double twopi = 6.283185307179586476925286766559005768394; 0036 static const double inv_halfpi = 0.6366197723675813430755350534900574; 0037 0038 static const int jrll[] = { 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4 }; 0039 static const int jpll[] = { 1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7 }; 0040 0041 static const short ctab[] = 0042 { 0043 #define Z(a) a,a+1,a+256,a+257 0044 #define Y(a) Z(a),Z(a+2),Z(a+512),Z(a+514) 0045 #define X(a) Y(a),Y(a+4),Y(a+1024),Y(a+1028) 0046 X(0), X(8), X(2048), X(2056) 0047 #undef X 0048 #undef Y 0049 #undef Z 0050 }; 0051 0052 static const short utab[] = 0053 { 0054 #define Z(a) 0x##a##0, 0x##a##1, 0x##a##4, 0x##a##5 0055 #define Y(a) Z(a##0), Z(a##1), Z(a##4), Z(a##5) 0056 #define X(a) Y(a##0), Y(a##1), Y(a##4), Y(a##5) 0057 X(0), X(1), X(4), X(5) 0058 #undef X 0059 #undef Y 0060 #undef Z 0061 }; 0062 0063 static short xoffset[] = { -1, -1, 0, 1, 1, 1, 0, -1 }; 0064 static short yoffset[] = { 0, 1, 1, 1, 0, -1, -1, -1 }; 0065 0066 static short facearray[9][12] = 0067 { 0068 { 8, 9, 10, 11, -1, -1, -1, -1, 10, 11, 8, 9 }, // S 0069 { 5, 6, 7, 4, 8, 9, 10, 11, 9, 10, 11, 8 }, // SE 0070 { -1, -1, -1, -1, 5, 6, 7, 4, -1, -1, -1, -1 }, // E 0071 { 4, 5, 6, 7, 11, 8, 9, 10, 11, 8, 9, 10 }, // SW 0072 { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 }, // center 0073 { 1, 2, 3, 0, 0, 1, 2, 3, 5, 6, 7, 4 }, // NE 0074 { -1, -1, -1, -1, 7, 4, 5, 6, -1, -1, -1, -1 }, // W 0075 { 3, 0, 1, 2, 3, 0, 1, 2, 4, 5, 6, 7 }, // NW 0076 { 2, 3, 0, 1, -1, -1, -1, -1, 0, 1, 2, 3 } 0077 }; // N 0078 0079 static uchar swaparray[9][12] = 0080 { 0081 { 0, 0, 3 }, // S 0082 { 0, 0, 6 }, // SE 0083 { 0, 0, 0 }, // E 0084 { 0, 0, 5 }, // SW 0085 { 0, 0, 0 }, // center 0086 { 5, 0, 0 }, // NE 0087 { 0, 0, 0 }, // W 0088 { 6, 0, 0 }, // NW 0089 { 3, 0, 0 } 0090 }; // N 0091 0092 static double fmodulo(double v1, double v2) 0093 { 0094 if (v1 >= 0) 0095 { 0096 return (v1 < v2) ? v1 : fmod(v1, v2); 0097 } 0098 0099 double tmp = fmod(v1, v2) + v2; 0100 return (tmp == v2) ? 0. : tmp; 0101 } 0102 0103 0104 void HEALPix::getCornerPoints(int level, int pix, SkyPoint *skyCoords) 0105 { 0106 QVector3D v[4]; 0107 // Transform from HealPIX convention to KStars 0108 QVector3D transformed[4]; 0109 0110 int nside = 1 << level; 0111 boundaries(nside, pix, 1, v); 0112 0113 // From rectangular coordinates to Sky coordinates 0114 for (int i = 0; i < 4; i++) 0115 { 0116 transformed[i].setX(v[i].x()); 0117 transformed[i].setY(v[i].y()); 0118 transformed[i].setZ(v[i].z()); 0119 0120 double ra = 0, de = 0; 0121 xyz2sph(transformed[i], ra, de); 0122 de /= dms::DegToRad; 0123 ra /= dms::DegToRad; 0124 0125 if (HIPSManager::Instance()->getCurrentFrame() == HIPSManager::HIPS_GALACTIC_FRAME) 0126 { 0127 dms galacticLong(ra); 0128 dms galacticLat(de); 0129 skyCoords[i].GalacticToEquatorial1950(&galacticLong, &galacticLat); 0130 skyCoords[i].B1950ToJ2000(); 0131 skyCoords[i].setRA0(skyCoords[i].ra()); 0132 skyCoords[i].setDec0(skyCoords[i].dec()); 0133 } 0134 else 0135 { 0136 skyCoords[i].setRA0(ra / 15.0); 0137 skyCoords[i].setDec0(de); 0138 } 0139 0140 skyCoords[i].updateCoords(KStarsData::Instance()->updateNum(), false); 0141 skyCoords[i].EquatorialToHorizontal(KStarsData::Instance()->lst(), KStarsData::Instance()->geo()->lat()); 0142 } 0143 0144 } 0145 0146 void HEALPix::boundaries(qint32 nside, qint32 pix, int step, QVector3D *out) 0147 { 0148 int ix, iy, fn; 0149 0150 nest2xyf(nside, pix, &ix, &iy, &fn); 0151 0152 double dc = 0.5 / nside; 0153 double xc = (ix + 0.5) / nside; 0154 double yc = (iy + 0.5) / nside; 0155 double d = 1. / (step * nside); 0156 0157 for (int i = 0; i < step; ++i) 0158 { 0159 out[0] = toVec3(xc + dc - i * d, yc + dc, fn); 0160 out[1] = toVec3(xc - dc, yc + dc - i * d, fn); 0161 out[2] = toVec3(xc - dc + i * d, yc - dc, fn); 0162 out[3] = toVec3(xc + dc, yc - dc + i * d, fn); 0163 } 0164 } 0165 0166 QVector3D HEALPix::toVec3(double fx, double fy, int face) 0167 { 0168 double jr = jrll[face] - fx - fy; 0169 0170 double locz; 0171 double locsth; 0172 double locphi; 0173 bool lochave_sth = false; 0174 0175 double nr; 0176 if (jr < 1) 0177 { 0178 nr = jr; 0179 double tmp = nr * nr / 3.; 0180 locz = 1 - tmp; 0181 if (locz > 0.99) 0182 { 0183 locsth = sqrt(tmp * (2. - tmp)); 0184 lochave_sth = true; 0185 } 0186 } 0187 else if (jr > 3) 0188 { 0189 nr = 4 - jr; 0190 double tmp = nr * nr / 3.; 0191 locz = tmp - 1; 0192 if (locz < -0.99) 0193 { 0194 locsth = sqrt(tmp * (2. - tmp)); 0195 lochave_sth = true; 0196 } 0197 } 0198 else 0199 { 0200 nr = 1; 0201 locz = (2 - jr) * 2. / 3.; 0202 } 0203 0204 double tmp = jpll[face] * nr + fx - fy; 0205 if (tmp < 0) tmp += 8; 0206 if (tmp >= 8) tmp -= 8; 0207 locphi = (nr < 1e-15) ? 0 : (0.5 * dms::PI / 2.0 * tmp) / nr; 0208 0209 double st = lochave_sth ? locsth : sqrt((1.0 - locz) * (1.0 + locz)); 0210 0211 QVector3D out; 0212 0213 out.setX(st * cos(locphi)); 0214 out.setY(st * sin(locphi)); 0215 out.setZ(locz); 0216 0217 return out; 0218 } 0219 0220 void HEALPix::nest2xyf(int nside, int pix, int *ix, int *iy, int *face_num) 0221 { 0222 int npface_ = nside * nside, raw; 0223 *face_num = pix / npface_; 0224 pix &= (npface_ - 1); 0225 raw = (pix & 0x5555) | ((pix & 0x55550000) >> 15); 0226 *ix = ctab[raw & 0xff] | (ctab[raw >> 8] << 4); 0227 pix >>= 1; 0228 raw = (pix & 0x5555) | ((pix & 0x55550000) >> 15); 0229 *iy = ctab[raw & 0xff] | (ctab[raw >> 8] << 4); 0230 } 0231 0232 static int64_t spread_bits64 (int v) 0233 { 0234 return (int64_t)(utab[ v & 0xff]) 0235 | ((int64_t)(utab[(v >> 8) & 0xff]) << 16) 0236 | ((int64_t)(utab[(v >> 16) & 0xff]) << 32) 0237 | ((int64_t)(utab[(v >> 24) & 0xff]) << 48); 0238 } 0239 0240 0241 static int xyf2nest (int nside, int ix, int iy, int face_num) 0242 { 0243 return (face_num * nside * nside) + 0244 (utab[ix & 0xff] | (utab[ix >> 8] << 16) 0245 | (utab[iy & 0xff] << 1) | (utab[iy >> 8] << 17)); 0246 } 0247 0248 0249 /* 0250 int static leadingZeros(qint32 value) 0251 { 0252 int leadingZeros = 0; 0253 while(value != 0) 0254 { 0255 value = value >> 1; 0256 leadingZeros++; 0257 } 0258 0259 return (32 - leadingZeros); 0260 } 0261 */ 0262 0263 /* 0264 static int ilog2(qint32 arg) 0265 { 0266 return 32 - leadingZeros(qMax(arg, 1)); 0267 } 0268 */ 0269 0270 static int nside2order(qint32 nside) 0271 { 0272 { 0273 int i = 0; 0274 while((nside >> (++i)) > 0); 0275 return --i; 0276 } 0277 } 0278 0279 /** Returns the neighboring pixels of ipix. 0280 This method works in both RING and NEST schemes, but is 0281 considerably faster in the NEST scheme. 0282 @param nside defines the size of the neighboring area 0283 @param ipix the requested pixel number. 0284 @param result the result array 0285 @return array with indices of the neighboring pixels. 0286 The returned array contains (in this order) 0287 the pixel numbers of the SW, W, NW, N, NE, E, SE and S neighbor 0288 of ipix. If a neighbor does not exist (this can only happen 0289 for the W, N, E and S neighbors), its entry is set to -1. */ 0290 0291 void HEALPix::neighbours(int nside, qint32 ipix, int *result) 0292 { 0293 int ix, iy, face_num; 0294 0295 int order = nside2order(nside); 0296 0297 nest2xyf(nside, ipix, &ix, &iy, &face_num); 0298 0299 qint32 nsm1 = nside - 1; 0300 if ((ix > 0) && (ix < nsm1) && (iy > 0) && (iy < nsm1)) 0301 { 0302 qint32 fpix = (qint32)(face_num) << (2 * order), 0303 px0 = spread_bits64(ix ), py0 = spread_bits64(iy ) << 1, 0304 pxp = spread_bits64(ix + 1), pyp = spread_bits64(iy + 1) << 1, 0305 pxm = spread_bits64(ix - 1), pym = spread_bits64(iy - 1) << 1; 0306 0307 result[0] = fpix + pxm + py0; 0308 result[1] = fpix + pxm + pyp; 0309 result[2] = fpix + px0 + pyp; 0310 result[3] = fpix + pxp + pyp; 0311 result[4] = fpix + pxp + py0; 0312 result[5] = fpix + pxp + pym; 0313 result[6] = fpix + px0 + pym; 0314 result[7] = fpix + pxm + pym; 0315 } 0316 else 0317 { 0318 for (int i = 0; i < 8; ++i) 0319 { 0320 int x = ix + xoffset[i]; 0321 int y = iy + yoffset[i]; 0322 int nbnum = 4; 0323 if (x < 0) 0324 { 0325 x += nside; 0326 nbnum -= 1; 0327 } 0328 else if (x >= nside) 0329 { 0330 x -= nside; 0331 nbnum += 1; 0332 } 0333 if (y < 0) 0334 { 0335 y += nside; 0336 nbnum -= 3; 0337 } 0338 else if (y >= nside) 0339 { 0340 y -= nside; 0341 nbnum += 3; 0342 } 0343 0344 int f = facearray[nbnum][face_num]; 0345 0346 if (f >= 0) 0347 { 0348 int bits = swaparray[nbnum][face_num >> 2]; 0349 if ((bits & 1) > 0) x = (int)(nside - x - 1); 0350 if ((bits & 2) > 0) y = (int)(nside - y - 1); 0351 if ((bits & 4) > 0) 0352 { 0353 int tint = x; 0354 x = y; 0355 y = tint; 0356 } 0357 result[i] = xyf2nest(nside, x, y, f); 0358 } 0359 else 0360 result[i] = -1; 0361 } 0362 } 0363 } 0364 0365 int HEALPix::ang2pix_nest_z_phi (qint32 nside_, double z, double phi) 0366 { 0367 double za = fabs(z); 0368 double tt = fmodulo(phi, twopi) * inv_halfpi; /* in [0,4) */ 0369 int face_num, ix, iy; 0370 0371 if (za <= twothird) /* Equatorial region */ 0372 { 0373 double temp1 = nside_ * (0.5 + tt); 0374 double temp2 = nside_ * (z * 0.75); 0375 int jp = (int)(temp1 - temp2); /* index of ascending edge line */ 0376 int jm = (int)(temp1 + temp2); /* index of descending edge line */ 0377 int ifp = jp / nside_; /* in {0,4} */ 0378 int ifm = jm / nside_; 0379 face_num = (ifp == ifm) ? (ifp | 4) : ((ifp < ifm) ? ifp : (ifm + 8)); 0380 0381 ix = jm & (nside_ - 1); 0382 iy = nside_ - (jp & (nside_ - 1)) - 1; 0383 } 0384 else /* polar region, za > 2/3 */ 0385 { 0386 int ntt = (int)tt, jp, jm; 0387 double tp, tmp; 0388 if (ntt >= 4) ntt = 3; 0389 tp = tt - ntt; 0390 tmp = nside_ * sqrt(3 * (1 - za)); 0391 0392 jp = (int)(tp * tmp); /* increasing edge line index */ 0393 jm = (int)((1.0 - tp) * tmp); /* decreasing edge line index */ 0394 if (jp >= nside_) jp = nside_ - 1; /* for points too close to the boundary */ 0395 if (jm >= nside_) jm = nside_ - 1; 0396 if (z >= 0) 0397 { 0398 face_num = ntt; /* in {0,3} */ 0399 ix = nside_ - jm - 1; 0400 iy = nside_ - jp - 1; 0401 } 0402 else 0403 { 0404 face_num = ntt + 8; /* in {8,11} */ 0405 ix = jp; 0406 iy = jm; 0407 } 0408 } 0409 0410 return xyf2nest(nside_, ix, iy, face_num); 0411 } 0412 0413 int HEALPix::getPix(int level, double ra, double dec) 0414 { 0415 int nside = 1 << level; 0416 double polar[2] = { 0, 0 }; 0417 0418 if (HIPSManager::Instance()->getCurrentFrame() == HIPSManager::HIPS_GALACTIC_FRAME) 0419 { 0420 static QMatrix4x4 gl(-0.0548762f, -0.873437f, -0.483835f, 0, 0421 0.4941100f, -0.444830f, 0.746982f, 0, 0422 -0.8676660f, -0.198076f, 0.455984f, 0, 0423 0, 0, 0, 1); 0424 0425 double rcb = cos(dec); 0426 QVector3D xyz = QVector3D(rcb * cos(ra), rcb * sin(ra), sin(dec)); 0427 0428 xyz = gl.mapVector(xyz); 0429 xyz2sph(xyz, polar[1], polar[0]); 0430 } 0431 else 0432 { 0433 polar[0] = dec; 0434 polar[1] = ra; 0435 } 0436 0437 return ang2pix_nest_z_phi(nside, sin(polar[0]), polar[1]); 0438 } 0439 0440 void HEALPix::getPixChilds(int pix, int *childs) 0441 { 0442 childs[0] = pix * 4 + 0; 0443 childs[1] = pix * 4 + 1; 0444 childs[2] = pix * 4 + 2; 0445 childs[3] = pix * 4 + 3; 0446 } 0447 0448 void HEALPix::xyz2sph(const QVector3D &vec, double &l, double &b) 0449 { 0450 double rho = vec.x() * vec.x() + vec.y() * vec.y(); 0451 0452 if (rho > 0) 0453 { 0454 l = atan2(vec.y(), vec.x()); 0455 l -= floor(l / (M_PI * 2)) * (M_PI * 2); 0456 b = atan2(vec.z(), sqrt(rho)); 0457 } 0458 else 0459 { 0460 l = 0.0; 0461 if (vec.z() == 0.0) 0462 { 0463 b = 0.0; 0464 } 0465 else 0466 { 0467 b = (vec.z() > 0.0) ? M_PI / 2. : -dms::PI / 2.; 0468 } 0469 } 0470 } 0471 0472 0473