File indexing completed on 2024-03-24 15:17:56
0001 /* 0002 SPDX-FileCopyrightText: 2001-2005 Jason Harris <jharris@30doradus.org> 0003 SPDX-FileCopyrightText: 2004-2005 Pablo de Vicente <p.devicente@wanadoo.es> 0004 0005 SPDX-License-Identifier: GPL-2.0-or-later 0006 */ 0007 0008 #include "skypoint.h" 0009 0010 #include "dms.h" 0011 #include "ksnumbers.h" 0012 #include "kstarsdatetime.h" 0013 #include "kssun.h" 0014 #include "kstarsdata.h" 0015 #include "Options.h" 0016 #include "skyobject.h" 0017 #include "skycomponents/skymapcomposite.h" 0018 0019 #include "config-kstars.h" 0020 0021 #include <KLocalizedString> 0022 0023 #include <QDebug> 0024 0025 #include <cmath> 0026 0027 0028 0029 // N.B. To use libnova depending on the availability of the package, 0030 // uncomment the following: 0031 /* 0032 #ifdef HAVE_LIBNOVA 0033 #define SKYPOINT_USE_LIBNOVA 1 0034 #endif 0035 */ 0036 0037 #ifdef SKYPOINT_USE_LIBNOVA 0038 #include <libnova/libnova.h> 0039 bool SkyPoint::implementationIsLibnova = true; 0040 #else 0041 bool SkyPoint::implementationIsLibnova = false; 0042 #endif 0043 0044 #ifdef PROFILE_COORDINATE_CONVERSION 0045 #include <ctime> // For profiling, remove if not profiling. 0046 long unsigned SkyPoint::eqToHzCalls = 0; 0047 double SkyPoint::cpuTime_EqToHz = 0.; 0048 #endif 0049 0050 KSSun *SkyPoint::m_Sun = nullptr; 0051 const double SkyPoint::altCrit = -1.0; 0052 0053 SkyPoint::SkyPoint() 0054 { 0055 // Default constructor. Set nonsense values 0056 RA0.setD(-1); // RA >= 0 always :-) 0057 Dec0.setD(180); // Dec is between -90 and 90 Degrees :-) 0058 RA = RA0; 0059 Dec = Dec0; 0060 lastPrecessJD = J2000; // By convention, we use J2000 coordinates 0061 } 0062 0063 void SkyPoint::set(const dms &r, const dms &d) 0064 { 0065 RA0 = RA = r; 0066 Dec0 = Dec = d; 0067 lastPrecessJD = J2000; // By convention, we use J2000 coordinates 0068 } 0069 0070 void SkyPoint::EquatorialToHorizontal(const dms *LST, const dms *lat) 0071 { 0072 // qDebug() << Q_FUNC_INFO << "NOTE: This EquatorialToHorizontal overload (using dms pointers instead of CachingDms pointers) is deprecated and should be replaced with CachingDms prototype wherever speed is desirable!"; 0073 CachingDms _LST(*LST), _lat(*lat); 0074 EquatorialToHorizontal(&_LST, &_lat); 0075 } 0076 0077 void SkyPoint::EquatorialToHorizontal(const CachingDms *LST, const CachingDms *lat) 0078 { 0079 #ifdef PROFILE_COORDINATE_CONVERSION 0080 std::clock_t start = std::clock(); 0081 #endif 0082 //Uncomment for spherical trig version 0083 double AltRad, AzRad; 0084 double sindec, cosdec, sinlat, coslat, sinHA, cosHA; 0085 double sinAlt, cosAlt; 0086 0087 CachingDms HourAngle = 0088 (*LST) - ra(); // Using CachingDms subtraction operator to find cos/sin of HourAngle without calling sincos() 0089 0090 lat->SinCos(sinlat, coslat); 0091 dec().SinCos(sindec, cosdec); 0092 HourAngle.SinCos(sinHA, cosHA); 0093 0094 sinAlt = sindec * sinlat + cosdec * coslat * cosHA; 0095 AltRad = asin(sinAlt); 0096 0097 cosAlt = sqrt( 0098 1 - 0099 sinAlt * 0100 sinAlt); // Avoid trigonometric function. Return value of asin is always in [-pi/2, pi/2] and in this domain cosine is always non-negative, so we can use this. 0101 if (cosAlt == 0.) 0102 cosAlt = cos(AltRad); 0103 0104 double arg = (sindec - sinlat * sinAlt) / (coslat * cosAlt); 0105 if (arg <= -1.0) 0106 AzRad = dms::PI; 0107 else if (arg >= 1.0) 0108 AzRad = 0.0; 0109 else 0110 AzRad = acos(arg); 0111 0112 if (sinHA > 0.0 && AzRad != 0.0) 0113 AzRad = 2.0 * dms::PI - AzRad; // resolve acos() ambiguity 0114 0115 Alt.setRadians(AltRad); 0116 Az.setRadians(AzRad); 0117 #ifdef PROFILE_COORDINATE_CONVERSION 0118 std::clock_t stop = std::clock(); 0119 cpuTime_EqToHz += double(stop - start) / double(CLOCKS_PER_SEC); // Accumulate time in seconds 0120 ++eqToHzCalls; 0121 #endif 0122 0123 // //Uncomment for XYZ version 0124 // double xr, yr, zr, xr1, zr1, sa, ca; 0125 // //Z-axis rotation by -LST 0126 // dms a = dms( -1.0*LST->Degrees() ); 0127 // a.SinCos( sa, ca ); 0128 // xr1 = m_X*ca + m_Y*sa; 0129 // yr = -1.0*m_X*sa + m_Y*ca; 0130 // zr1 = m_Z; 0131 // 0132 // //Y-axis rotation by lat - 90. 0133 // a = dms( lat->Degrees() - 90.0 ); 0134 // a.SinCos( sa, ca ); 0135 // xr = xr1*ca - zr1*sa; 0136 // zr = xr1*sa + zr1*ca; 0137 // 0138 // //FIXME: eventually, we will work with XYZ directly 0139 // Alt.setRadians( asin( zr ) ); 0140 // Az.setRadians( atan2( yr, xr ) ); 0141 } 0142 0143 void SkyPoint::HorizontalToEquatorial(const dms *LST, const dms *lat) 0144 { 0145 double HARad, DecRad; 0146 double sinlat, coslat, sinAlt, cosAlt, sinAz, cosAz; 0147 double sinDec, cosDec; 0148 0149 lat->SinCos(sinlat, coslat); 0150 alt().SinCos(sinAlt, cosAlt); 0151 Az.SinCos(sinAz, cosAz); 0152 0153 sinDec = sinAlt * sinlat + cosAlt * coslat * cosAz; 0154 DecRad = asin(sinDec); 0155 cosDec = cos(DecRad); 0156 Dec.setRadians(DecRad); 0157 0158 double x = (sinAlt - sinlat * sinDec) / (coslat * cosDec); 0159 0160 //Under certain circumstances, x can be very slightly less than -1.0000, or slightly 0161 //greater than 1.0000, leading to a crash on acos(x). However, the value isn't 0162 //*really* out of range; it's a kind of roundoff error. 0163 if (x < -1.0 && x > -1.000001) 0164 HARad = dms::PI; 0165 else if (x > 1.0 && x < 1.000001) 0166 HARad = 0.0; 0167 else if (x < -1.0) 0168 { 0169 //qWarning() << "Coordinate out of range."; 0170 HARad = dms::PI; 0171 } 0172 else if (x > 1.0) 0173 { 0174 //qWarning() << "Coordinate out of range."; 0175 HARad = 0.0; 0176 } 0177 else 0178 HARad = acos(x); 0179 0180 if (sinAz > 0.0) 0181 HARad = 2.0 * dms::PI - HARad; // resolve acos() ambiguity 0182 0183 RA.setRadians(LST->radians() - HARad); 0184 RA.reduceToRange(dms::ZERO_TO_2PI); 0185 } 0186 0187 void SkyPoint::findEcliptic(const CachingDms *Obliquity, dms &EcLong, dms &EcLat) 0188 { 0189 double sinRA, cosRA, sinOb, cosOb, sinDec, cosDec; 0190 ra().SinCos(sinRA, cosRA); 0191 dec().SinCos(sinDec, cosDec); 0192 Obliquity->SinCos(sinOb, cosOb); 0193 0194 double ycosDec = sinRA * cosOb * cosDec + sinDec * sinOb; 0195 double ELongRad = atan2(ycosDec, cosDec * cosRA); 0196 EcLong.setRadians(ELongRad); 0197 EcLong.reduceToRange(dms::ZERO_TO_2PI); 0198 EcLat.setRadians(asin(sinDec * cosOb - cosDec * sinOb * sinRA)); // FIXME: Haversine 0199 } 0200 0201 void SkyPoint::setFromEcliptic(const CachingDms *Obliquity, const dms &EcLong, const dms &EcLat) 0202 { 0203 double sinLong, cosLong, sinLat, cosLat, sinObliq, cosObliq; 0204 EcLong.SinCos(sinLong, cosLong); 0205 EcLat.SinCos(sinLat, cosLat); 0206 Obliquity->SinCos(sinObliq, cosObliq); 0207 0208 // double sinDec = sinLat * cosObliq + cosLat * sinObliq * sinLong; 0209 0210 double ycosLat = sinLong * cosObliq * cosLat - sinLat * sinObliq; 0211 // double RARad = atan2( y, cosLong ); 0212 RA.setUsing_atan2(ycosLat, cosLong * cosLat); 0213 RA.reduceToRange(dms::ZERO_TO_2PI); 0214 // Dec.setUsing_asin(sinDec); 0215 0216 // Use Haversine to set declination 0217 Dec.setRadians(dms::PI / 2.0 - 2.0 * asin(sqrt(0.5 * ( 0218 1.0 - sinLat * cosObliq 0219 - cosLat * sinObliq * sinLong 0220 )))); 0221 } 0222 0223 void SkyPoint::precess(const KSNumbers *num) 0224 { 0225 double cosRA0, sinRA0, cosDec0, sinDec0; 0226 const Eigen::Matrix3d &precessionMatrix = num->p2(); 0227 Eigen::Vector3d v, s; 0228 0229 RA0.SinCos(sinRA0, cosRA0); 0230 Dec0.SinCos(sinDec0, cosDec0); 0231 0232 s[0] = cosRA0 * cosDec0; 0233 s[1] = sinRA0 * cosDec0; 0234 s[2] = sinDec0; 0235 0236 // NOTE: Rotation matrices are the fastest way to do rotations on 0237 // a vector. Quaternions need more multiplications. The rotation 0238 // matrix compensates in some sense by having more 'precomputed' 0239 // multiplications. The matrix elements seem to cache nicely, so 0240 // there isn't much overhead in accessing them. 0241 0242 //Multiply P2 and s to get v, the vector representing the new coords. 0243 // for ( unsigned int i=0; i<3; ++i ) { 0244 // v[i] = 0.0; 0245 // for (uint j=0; j< 3; ++j) { 0246 // v[i] += num->p2( j, i )*s[j]; 0247 // } 0248 // } 0249 v.noalias() = precessionMatrix * s; 0250 0251 //Extract RA, Dec from the vector: 0252 RA.setUsing_atan2(v[1], v[0]); 0253 RA.reduceToRange(dms::ZERO_TO_2PI); 0254 Dec.setUsing_asin(v[2]); 0255 } 0256 0257 SkyPoint SkyPoint::deprecess(const KSNumbers *num, long double epoch) 0258 { 0259 SkyPoint p1(RA, Dec); 0260 long double now = num->julianDay(); 0261 p1.precessFromAnyEpoch(now, epoch); 0262 if ((std::isnan(RA0.Degrees()) || std::isnan(Dec0.Degrees())) || 0263 (!std::isnan(Dec0.Degrees()) && fabs(Dec0.Degrees()) > 90.0)) 0264 { 0265 // We have invalid RA0 and Dec0, so set them if epoch = J2000. Otherwise, do not touch. 0266 if (epoch == J2000L) 0267 { 0268 RA0 = p1.ra(); 0269 Dec0 = p1.dec(); 0270 } 0271 } 0272 return p1; 0273 } 0274 0275 void SkyPoint::nutate(const KSNumbers *num, const bool reverse) 0276 { 0277 //Step 2: Nutation 0278 if (fabs(Dec.Degrees()) < 80.0) //approximate method 0279 { 0280 double cosRA, sinRA, cosDec, sinDec, tanDec; 0281 double cosOb, sinOb; 0282 double dRA, dDec; 0283 0284 RA.SinCos(sinRA, cosRA); 0285 Dec.SinCos(sinDec, cosDec); 0286 0287 tanDec = sinDec / cosDec; 0288 0289 // Equ 23.1 in Jean Meeus' book 0290 // nut_ecliptic / num->obliquity() is called epsilon in Meeus 0291 // nut.longitude / num->dEcLong() is called delta psi in Meeus 0292 // nut.obliquity / num->dObliq() is called delta epsilon in Meeus 0293 // Meeus notes that these expressions are invalid if the star is 0294 // close to one of the celestial poles 0295 0296 // N.B. These expressions are valid for FK5 coordinates 0297 // (presumably also valid for ICRS by extension), but not for 0298 // FK4. See the "Important Remark" on Page 152 of Jean Meeus' 0299 // book. 0300 0301 #ifdef SKYPOINT_USE_LIBNOVA 0302 // code lifted from libnova ln_get_equ_nut, tailored to our needs 0303 // with the option to add or remove nutation 0304 struct ln_nutation nut; 0305 ln_get_nutation (num->julianDay(), &nut); // FIXME: Is this cached, or is it a slow call? If it is slow, we should move it to KSNumbers 0306 0307 double nut_ecliptic = ln_deg_to_rad(nut.ecliptic + nut.obliquity); 0308 sinOb = sin(nut_ecliptic); 0309 cosOb = cos(nut_ecliptic); 0310 0311 dRA = nut.longitude * (cosOb + sinOb * sinRA * tanDec) - nut.obliquity * cosRA * tanDec; 0312 dDec = nut.longitude * (sinOb * cosRA) + nut.obliquity * sinRA; 0313 #else 0314 num->obliquity()->SinCos(sinOb, cosOb); 0315 0316 // N.B. num->dEcLong() and num->dObliq() methods return in 0317 // degrees, whereby the corrections dRA and dDec are also in 0318 // degrees. 0319 dRA = num->dEcLong() * (cosOb + sinOb * sinRA * tanDec) - num->dObliq() * cosRA * tanDec; 0320 dDec = num->dEcLong() * (sinOb * cosRA) + num->dObliq() * sinRA; 0321 #endif 0322 // the sign changed to remove nutation 0323 if (reverse) 0324 { 0325 dRA = -dRA; 0326 dDec = -dDec; 0327 } 0328 RA.setD(RA.Degrees() + dRA); 0329 Dec.setD(Dec.Degrees() + dDec); 0330 } 0331 else //exact method 0332 { 0333 // NOTE: Meeus declares that you must add Δψ to the ecliptic 0334 // longitude of the body to get a more accurate precession 0335 // result, but fails to emphasize that the NCP of the two 0336 // coordinates systems is different. The (RA, Dec) without 0337 // nutation computed, i.e. the mean place of the sky point is 0338 // referenced to the un-nutated geocentric frame (without the 0339 // obliquity correction), whereas the (RA, Dec) after nutation 0340 // is applied is referenced to the nutated geocentric 0341 // frame. This is more clearly explained in the "Explanatory 0342 // Supplement to the Astronomical Almanac" by 0343 // K. P. Seidelmann, which can be borrowed on the internet 0344 // archive as of this writing, see page 114: 0345 // https://archive.org/details/explanatorysuppl00pken 0346 // 0347 // The rotation matrix formulation in (3.222-3) and the figure 0348 // 3.222.1 make this clear 0349 0350 // TODO apply reverse test above 80 degrees 0351 dms EcLong, EcLat; 0352 CachingDms obliquityWithoutNutation(*num->obliquity() - dms(num->dObliq())); 0353 0354 if (reverse) 0355 { 0356 //Subtract dEcLong from the Ecliptic Longitude 0357 findEcliptic(num->obliquity(), EcLong, EcLat); 0358 dms newLong(EcLong.Degrees() - num->dEcLong()); 0359 setFromEcliptic(&obliquityWithoutNutation, newLong, EcLat); // FIXME: Check 0360 } 0361 else 0362 { 0363 //Add dEcLong to the Ecliptic Longitude 0364 findEcliptic(&obliquityWithoutNutation, EcLong, EcLat); 0365 dms newLong(EcLong.Degrees() + num->dEcLong()); 0366 setFromEcliptic(num->obliquity(), newLong, EcLat); 0367 } 0368 } 0369 } 0370 0371 SkyPoint SkyPoint::moveAway(const SkyPoint &from, double dist) const 0372 { 0373 CachingDms lat1, dtheta; 0374 0375 if (dist == 0.0) 0376 { 0377 qDebug() << Q_FUNC_INFO << "moveAway called with zero distance!"; 0378 return *this; 0379 } 0380 0381 double dst = fabs(dist * dms::DegToRad / 3600.0); // In radian 0382 0383 // Compute the bearing angle w.r.t. the RA axis ("latitude") 0384 CachingDms dRA(ra() - from.ra()); 0385 CachingDms dDec(dec() - from.dec()); 0386 double bearing = atan2(dRA.sin() / dRA.cos(), dDec.sin()); // Do not use dRA = PI / 2!! 0387 //double bearing = atan2( dDec.radians() , dRA.radians() ); 0388 0389 // double dir0 = (dist >= 0 ) ? bearing : bearing + dms::PI; // in radian 0390 double dir0 = bearing + std::signbit(dist) * dms::PI; // might be faster? 0391 double sinDst = sin(dst), cosDst = cos(dst); 0392 0393 lat1.setUsing_asin(dec().sin() * cosDst + dec().cos() * sinDst * cos(dir0)); 0394 dtheta.setUsing_atan2(sin(dir0) * sinDst * dec().cos(), cosDst - dec().sin() * lat1.sin()); 0395 0396 return SkyPoint(ra() + dtheta, lat1); 0397 } 0398 0399 bool SkyPoint::checkBendLight() 0400 { 0401 // First see if we are close enough to the sun to bother about the 0402 // gravitational lensing effect. We correct for the effect at 0403 // least till b = 10 solar radii, where the effect is only about 0404 // 0.06". Assuming min. sun-earth distance is 200 solar radii. 0405 static const dms maxAngle(1.75 * (30.0 / 200.0) / dms::DegToRad); 0406 0407 if (!m_Sun) 0408 { 0409 SkyComposite *skycomopsite = KStarsData::Instance()->skyComposite(); 0410 0411 if (skycomopsite == nullptr) 0412 return false; 0413 0414 m_Sun = dynamic_cast<KSSun *>(skycomopsite->findByName(i18n("Sun"))); 0415 0416 if (m_Sun == nullptr) 0417 return false; 0418 } 0419 0420 // TODO: This can be optimized further. We only need a ballpark estimate of the distance to the sun to start with. 0421 return (fabs(angularDistanceTo(static_cast<const SkyPoint *>(m_Sun)).Degrees()) <= 0422 maxAngle.Degrees()); // NOTE: dynamic_cast is slow and not important here. 0423 } 0424 0425 bool SkyPoint::bendlight() 0426 { 0427 // NOTE: This should be applied before aberration 0428 // NOTE: One must call checkBendLight() before unnecessarily calling this. 0429 // We correct for GR effects 0430 0431 // NOTE: This code is buggy. The sun needs to be initialized to 0432 // the current epoch -- but we are not certain that this is the 0433 // case. We have, as of now, no way of telling if the sun is 0434 // initialized or not. If we initialize the sun here, we will be 0435 // slowing down the program rather substantially and potentially 0436 // introducing bugs. Therefore, we just ignore this problem, and 0437 // hope that whenever the user is interested in seeing the effects 0438 // of GR, we have the sun initialized correctly. This is usually 0439 // the case. When the sun is not correctly initialized, rearth() 0440 // is not computed, so we just assume it is nominally equal to 1 0441 // AU to get a reasonable estimate. 0442 Q_ASSERT(m_Sun); 0443 double corr_sec = 1.75 * m_Sun->physicalSize() / 0444 ((std::isfinite(m_Sun->rearth()) ? m_Sun->rearth() : 1) * AU_KM * 0445 angularDistanceTo(static_cast<const SkyPoint *>(m_Sun)).sin()); 0446 Q_ASSERT(corr_sec > 0); 0447 0448 SkyPoint sp = moveAway(*m_Sun, corr_sec); 0449 setRA(sp.ra()); 0450 setDec(sp.dec()); 0451 return true; 0452 } 0453 0454 void SkyPoint::aberrate(const KSNumbers *num, bool reverse) 0455 { 0456 #ifdef SKYPOINT_USE_LIBNOVA 0457 ln_equ_posn pos { RA.Degrees(), Dec.Degrees() }; 0458 ln_equ_posn abPos { 0, 0 }; 0459 ln_get_equ_aber(&pos, num->julianDay(), &abPos); 0460 if (reverse) 0461 { 0462 RA.setD(RA.Degrees() * 2 - abPos.ra); 0463 Dec.setD(Dec.Degrees() * 2 - abPos.dec); 0464 } 0465 else 0466 { 0467 RA.setD(abPos.ra); 0468 Dec.setD(abPos.dec); 0469 } 0470 0471 #else 0472 // N.B. These expressions are valid for FK5 coordinates 0473 // (presumably also valid for ICRS by extension), but not for 0474 // FK4. See the "Important Remark" on Page 152 of Jean Meeus' 0475 // book. 0476 0477 // N.B. Even though Meeus does not say this explicitly, these 0478 // equations must not apply near the pole. Therefore, we fall-back 0479 // to the expressions provided by Meeus in ecliptic coordinates 0480 // (Equ 23.2) when we are near the pole. 0481 0482 double K = num->constAberr().Degrees(); //constant of aberration 0483 double e = num->earthEccentricity(); // eccentricity of Earth's orbit 0484 0485 if (fabs(Dec.Degrees()) < 80.0) 0486 { 0487 0488 double cosRA, sinRA, cosDec, sinDec; 0489 double cosL, sinL, cosP, sinP; 0490 double cosOb, sinOb; 0491 0492 0493 RA.SinCos(sinRA, cosRA); 0494 Dec.SinCos(sinDec, cosDec); 0495 0496 num->obliquity()->SinCos(sinOb, cosOb); 0497 // double tanOb = sinOb/cosOb; 0498 0499 num->sunTrueLongitude().SinCos(sinL, cosL); 0500 num->earthPerihelionLongitude().SinCos(sinP, cosP); 0501 0502 //Step 3: Aberration 0503 // These are the expressions given in Jean Meeus, Equation (23.3) 0504 0505 // double dRA = -1.0 * K * ( cosRA * cosL * cosOb + sinRA * sinL )/cosDec 0506 // + e * K * ( cosRA * cosP * cosOb + sinRA * sinP )/cosDec; 0507 0508 // double dDec = -1.0 * K * ( cosL * cosOb * ( tanOb * cosDec - sinRA * sinDec ) + cosRA * sinDec * sinL ) 0509 // + e * K * ( cosP * cosOb * ( tanOb * cosDec - sinRA * sinDec ) + cosRA * sinDec * sinP ); 0510 0511 // However, we have factorized the expressions below by pulling 0512 // out common factors to make it more optimized, in case the 0513 // compiler fails to spot these optimizations. 0514 0515 // N.B. I had failed to factor out the expressions correctly, 0516 // making mistakes, in c5e709bd91, which should now be 0517 // fixed. --asimha 0518 0519 // FIXME: Check if the unit tests have sufficient coverage to 0520 // check this expression 0521 0522 double dRA = (K / cosDec) * ( 0523 cosRA * cosOb * (e * cosP - cosL) 0524 + sinRA * (e * sinP - sinL) 0525 ); 0526 double dDec = K * ( 0527 (sinOb * cosDec - cosOb * sinDec * sinRA) * (e * cosP - cosL) 0528 + cosRA * sinDec * (e * sinP - sinL) 0529 ); 0530 0531 // N.B. Meeus points out that the result has the same units as 0532 // K, so the corrections are in degrees. 0533 0534 if (reverse) 0535 { 0536 dRA = -dRA; 0537 dDec = -dDec; 0538 } 0539 RA.setD(RA.Degrees() + dRA); 0540 Dec.setD(Dec.Degrees() + dDec); 0541 } 0542 else 0543 { 0544 dms EcLong, EcLat; 0545 double sinEcLat, cosEcLat; 0546 const auto &L = num->sunTrueLongitude(); 0547 const auto &P = num->earthPerihelionLongitude(); 0548 0549 findEcliptic(num->obliquity(), EcLong, EcLat); 0550 EcLat.SinCos(sinEcLat, cosEcLat); 0551 0552 double sin_L_minus_EcLong, cos_L_minus_EcLong; 0553 double sin_P_minus_EcLong, cos_P_minus_EcLong; 0554 (L - EcLong).SinCos(sin_L_minus_EcLong, cos_L_minus_EcLong); 0555 (P - EcLong).SinCos(sin_P_minus_EcLong, cos_P_minus_EcLong); 0556 0557 // Equation (23.2) in Meeus 0558 // N.B. dEcLong, dEcLat are in degrees, because K is expressed in degrees. 0559 double dEcLong = (K / cosEcLat) * (e * cos_P_minus_EcLong - cos_L_minus_EcLong); 0560 double dEcLat = K * sinEcLat * (e * sin_P_minus_EcLong - sin_L_minus_EcLong); 0561 0562 // Note: These are approximate corrections, so it is 0563 // appropriate to change their sign to reverse them to first 0564 // order in the corrections. As a result, the forward and 0565 // reverse calculations will not be exact inverses, but only 0566 // approximate inverses. 0567 if (reverse) 0568 { 0569 dEcLong = -dEcLong; 0570 dEcLat = -dEcLat; 0571 } 0572 0573 // Update the ecliptic coordinates to their new values 0574 EcLong.setD(EcLong.Degrees() + dEcLong); 0575 EcLat.setD(EcLat.Degrees() + dEcLat); 0576 setFromEcliptic(num->obliquity(), EcLong, EcLat); 0577 } 0578 #endif 0579 } 0580 0581 // Note: This method is one of the major rate determining factors in how fast the map pans / zooms in or out 0582 void SkyPoint::updateCoords(const KSNumbers *num, bool /*includePlanets*/, const CachingDms *lat, const CachingDms *LST, 0583 bool forceRecompute) 0584 { 0585 //Correct the catalog coordinates for the time-dependent effects 0586 //of precession, nutation and aberration 0587 bool recompute, lens; 0588 0589 // NOTE: The same short-circuiting checks are also implemented in 0590 // StarObject::JITUpdate(), even before calling 0591 // updateCoords(). While this is code-duplication, these bits of 0592 // code need to be really optimized, at least for stars. For 0593 // optimization purposes, the code is left duplicated in two 0594 // places. Please be wary of changing one without changing the 0595 // other. 0596 0597 Q_ASSERT(std::isfinite(lastPrecessJD)); 0598 0599 if (Options::useRelativistic() && checkBendLight()) 0600 { 0601 recompute = true; 0602 lens = true; 0603 } 0604 else 0605 { 0606 recompute = (Options::alwaysRecomputeCoordinates() || forceRecompute || 0607 std::abs(lastPrecessJD - num->getJD()) >= 0.00069444); // Update once per solar minute 0608 lens = false; 0609 } 0610 if (recompute) 0611 { 0612 precess(num); 0613 nutate(num); 0614 if (lens) 0615 bendlight(); // FIXME: Shouldn't we apply this on the horizontal coordinates? 0616 aberrate(num); 0617 lastPrecessJD = num->getJD(); 0618 Q_ASSERT(std::isfinite(RA.Degrees()) && std::isfinite(Dec.Degrees())); 0619 } 0620 0621 if (lat || LST) 0622 qWarning() << i18n("lat and LST parameters should only be used in KSPlanetBase objects."); 0623 } 0624 0625 void SkyPoint::precessFromAnyEpoch(long double jd0, long double jdf) 0626 { 0627 double cosRA, sinRA, cosDec, sinDec; 0628 double v[3], s[3]; 0629 0630 RA = RA0; 0631 Dec = Dec0; // Is this necessary? 0632 0633 if (jd0 == jdf) 0634 return; 0635 0636 RA.SinCos(sinRA, cosRA); 0637 Dec.SinCos(sinDec, cosDec); 0638 0639 if (jd0 == B1950L) 0640 { 0641 B1950ToJ2000(); 0642 jd0 = J2000L; 0643 RA.SinCos(sinRA, cosRA); 0644 Dec.SinCos(sinDec, cosDec); 0645 } 0646 0647 if (jd0 != jdf) 0648 { 0649 // The original coordinate is referred to the FK5 system and 0650 // is NOT J2000. 0651 if (jd0 != J2000L) 0652 { 0653 //v is a column vector representing input coordinates. 0654 v[0] = cosRA * cosDec; 0655 v[1] = sinRA * cosDec; 0656 v[2] = sinDec; 0657 0658 //Need to first precess to J2000.0 coords 0659 //s is the product of P1 and v; s represents the 0660 //coordinates precessed to J2000 0661 KSNumbers num(jd0); 0662 for (unsigned int i = 0; i < 3; ++i) 0663 { 0664 s[i] = num.p1(0, i) * v[0] + num.p1(1, i) * v[1] + num.p1(2, i) * v[2]; 0665 } 0666 0667 //Input coords already in J2000, set s accordingly. 0668 } 0669 else 0670 { 0671 s[0] = cosRA * cosDec; 0672 s[1] = sinRA * cosDec; 0673 s[2] = sinDec; 0674 } 0675 0676 if (jdf == B1950L) 0677 { 0678 RA.setRadians(atan2(s[1], s[0])); 0679 Dec.setRadians(asin(s[2])); 0680 J2000ToB1950(); 0681 0682 return; 0683 } 0684 0685 KSNumbers num(jdf); 0686 for (unsigned int i = 0; i < 3; ++i) 0687 { 0688 v[i] = num.p2(0, i) * s[0] + num.p2(1, i) * s[1] + num.p2(2, i) * s[2]; 0689 } 0690 0691 RA.setUsing_atan2(v[1], v[0]); 0692 Dec.setUsing_asin(v[2]); 0693 0694 RA.reduceToRange(dms::ZERO_TO_2PI); 0695 0696 return; 0697 } 0698 } 0699 0700 void SkyPoint::apparentCoord(long double jd0, long double jdf) 0701 { 0702 precessFromAnyEpoch(jd0, jdf); 0703 KSNumbers num(jdf); 0704 nutate(&num); 0705 if (Options::useRelativistic() && checkBendLight()) 0706 bendlight(); 0707 aberrate(&num); 0708 } 0709 0710 SkyPoint SkyPoint::catalogueCoord(long double jdf) 0711 { 0712 KSNumbers num(jdf); 0713 0714 // remove abberation 0715 aberrate(&num, true); 0716 0717 // remove nutation 0718 nutate(&num, true); 0719 0720 // remove precession 0721 // the start position needs to be in RA0,Dec0 0722 RA0 = RA; 0723 Dec0 = Dec; 0724 // from now to J2000 0725 precessFromAnyEpoch(jdf, static_cast<long double>(J2000)); 0726 // the J2000 position is in RA,Dec, move to RA0, Dec0 0727 RA0 = RA; 0728 Dec0 = Dec; 0729 lastPrecessJD = J2000; 0730 0731 SkyPoint sp(RA0, Dec0); 0732 return sp; 0733 } 0734 0735 void SkyPoint::Equatorial1950ToGalactic(dms &galLong, dms &galLat) 0736 { 0737 double a = 192.25; 0738 double sinb, cosb, sina_RA, cosa_RA, sinDEC, cosDEC, tanDEC; 0739 0740 dms c(303.0); 0741 dms b(27.4); 0742 tanDEC = tan(Dec.radians()); 0743 0744 b.SinCos(sinb, cosb); 0745 dms(a - RA.Degrees()).SinCos(sina_RA, cosa_RA); 0746 Dec.SinCos(sinDEC, cosDEC); 0747 0748 galLong.setRadians(c.radians() - atan2(sina_RA, cosa_RA * sinb - tanDEC * cosb)); 0749 galLong.reduceToRange(dms::ZERO_TO_2PI); 0750 0751 galLat.setRadians(asin(sinDEC * sinb + cosDEC * cosb * cosa_RA)); 0752 } 0753 0754 void SkyPoint::GalacticToEquatorial1950(const dms *galLong, const dms *galLat) 0755 { 0756 double a = 123.0; 0757 double sinb, cosb, singLat, cosgLat, tangLat, singLong_a, cosgLong_a; 0758 0759 dms c(12.25); 0760 dms b(27.4); 0761 tangLat = tan(galLat->radians()); 0762 galLat->SinCos(singLat, cosgLat); 0763 0764 dms(galLong->Degrees() - a).SinCos(singLong_a, cosgLong_a); 0765 b.SinCos(sinb, cosb); 0766 0767 RA.setRadians(c.radians() + atan2(singLong_a, cosgLong_a * sinb - tangLat * cosb)); 0768 RA.reduceToRange(dms::ZERO_TO_2PI); 0769 0770 Dec.setRadians(asin(singLat * sinb + cosgLat * cosb * cosgLong_a)); 0771 } 0772 0773 void SkyPoint::B1950ToJ2000(void) 0774 { 0775 double cosRA, sinRA, cosDec, sinDec; 0776 // double cosRA0, sinRA0, cosDec0, sinDec0; 0777 double v[3], s[3]; 0778 0779 // 1984 January 1 0h 0780 KSNumbers num(2445700.5L); 0781 0782 // Eterms due to aberration 0783 addEterms(); 0784 RA.SinCos(sinRA, cosRA); 0785 Dec.SinCos(sinDec, cosDec); 0786 0787 // Precession from B1950 to J1984 0788 s[0] = cosRA * cosDec; 0789 s[1] = sinRA * cosDec; 0790 s[2] = sinDec; 0791 for (unsigned int i = 0; i < 3; ++i) 0792 { 0793 v[i] = num.p2b(0, i) * s[0] + num.p2b(1, i) * s[1] + num.p2b(2, i) * s[2]; 0794 } 0795 0796 // RA zero-point correction at 1984 day 1, 0h. 0797 RA.setRadians(atan2(v[1], v[0])); 0798 Dec.setRadians(asin(v[2])); 0799 0800 RA.setH(RA.Hours() + 0.06390 / 3600.); 0801 RA.SinCos(sinRA, cosRA); 0802 Dec.SinCos(sinDec, cosDec); 0803 0804 s[0] = cosRA * cosDec; 0805 s[1] = sinRA * cosDec; 0806 s[2] = sinDec; 0807 0808 // Precession from 1984 to J2000. 0809 0810 for (unsigned int i = 0; i < 3; ++i) 0811 { 0812 v[i] = num.p1(0, i) * s[0] + num.p1(1, i) * s[1] + num.p1(2, i) * s[2]; 0813 } 0814 0815 RA.setRadians(atan2(v[1], v[0])); 0816 Dec.setRadians(asin(v[2])); 0817 } 0818 0819 void SkyPoint::J2000ToB1950(void) 0820 { 0821 double cosRA, sinRA, cosDec, sinDec; 0822 // double cosRA0, sinRA0, cosDec0, sinDec0; 0823 double v[3], s[3]; 0824 0825 // 1984 January 1 0h 0826 KSNumbers num(2445700.5L); 0827 0828 RA.SinCos(sinRA, cosRA); 0829 Dec.SinCos(sinDec, cosDec); 0830 0831 s[0] = cosRA * cosDec; 0832 s[1] = sinRA * cosDec; 0833 s[2] = sinDec; 0834 0835 // Precession from J2000 to 1984 day, 0h. 0836 0837 for (unsigned int i = 0; i < 3; ++i) 0838 { 0839 v[i] = num.p2(0, i) * s[0] + num.p2(1, i) * s[1] + num.p2(2, i) * s[2]; 0840 } 0841 0842 RA.setRadians(atan2(v[1], v[0])); 0843 Dec.setRadians(asin(v[2])); 0844 0845 // RA zero-point correction at 1984 day 1, 0h. 0846 0847 RA.setH(RA.Hours() - 0.06390 / 3600.); 0848 RA.SinCos(sinRA, cosRA); 0849 Dec.SinCos(sinDec, cosDec); 0850 0851 // Precession from B1950 to J1984 0852 0853 s[0] = cosRA * cosDec; 0854 s[1] = sinRA * cosDec; 0855 s[2] = sinDec; 0856 for (unsigned int i = 0; i < 3; ++i) 0857 { 0858 v[i] = num.p1b(0, i) * s[0] + num.p1b(1, i) * s[1] + num.p1b(2, i) * s[2]; 0859 } 0860 0861 RA.setRadians(atan2(v[1], v[0])); 0862 Dec.setRadians(asin(v[2])); 0863 0864 // Eterms due to aberration 0865 subtractEterms(); 0866 } 0867 0868 SkyPoint SkyPoint::Eterms(void) 0869 { 0870 double sd, cd, sinEterm, cosEterm; 0871 dms raTemp, raDelta, decDelta; 0872 0873 Dec.SinCos(sd, cd); 0874 raTemp.setH(RA.Hours() + 11.25); 0875 raTemp.SinCos(sinEterm, cosEterm); 0876 0877 raDelta.setH(0.0227 * sinEterm / (3600. * cd)); 0878 decDelta.setD(0.341 * cosEterm * sd / 3600. + 0.029 * cd / 3600.); 0879 0880 return SkyPoint(raDelta, decDelta); 0881 } 0882 0883 void SkyPoint::addEterms(void) 0884 { 0885 SkyPoint spd = Eterms(); 0886 0887 RA = RA + spd.ra(); 0888 Dec = Dec + spd.dec(); 0889 } 0890 0891 void SkyPoint::subtractEterms(void) 0892 { 0893 SkyPoint spd = Eterms(); 0894 0895 RA = RA - spd.ra(); 0896 Dec = Dec - spd.dec(); 0897 } 0898 0899 dms SkyPoint::angularDistanceTo(const SkyPoint *sp, double *const positionAngle) const 0900 { 0901 // double dalpha = sp->ra().radians() - ra().radians() ; 0902 // double ddelta = sp->dec().radians() - dec().radians(); 0903 CachingDms dalpha = sp->ra() - ra(); 0904 CachingDms ddelta = sp->dec() - dec(); 0905 0906 // double sa = sin(dalpha/2.); 0907 // double sd = sin(ddelta/2.); 0908 0909 // double hava = sa*sa; 0910 // double havd = sd*sd; 0911 0912 // Compute the haversin directly: 0913 double hava = (1 - dalpha.cos()) / 2.; 0914 double havd = (1 - ddelta.cos()) / 2.; 0915 0916 // Haversine law 0917 double aux = havd + (sp->dec().cos()) * dec().cos() * hava; 0918 0919 dms angDist; 0920 angDist.setRadians(2. * fabs(asin(sqrt(aux)))); 0921 0922 if (positionAngle) 0923 { 0924 // Also compute the position angle of the line from this SkyPoint to sp 0925 //*positionAngle = acos( tan(-ddelta)/tan( angDist.radians() ) ); // FIXME: Might fail for large ddelta / zero angDist 0926 //if( -dalpha < 0 ) 0927 // *positionAngle = 2*M_PI - *positionAngle; 0928 *positionAngle = 0929 atan2f(dalpha.sin(), (dec().cos()) * tan(sp->dec().radians()) - (dec().sin()) * dalpha.cos()) * 180 / M_PI; 0930 } 0931 return angDist; 0932 } 0933 0934 double SkyPoint::vRSun(long double jd0) 0935 { 0936 double ca, sa, cd, sd, vsun; 0937 double cosRA, sinRA, cosDec, sinDec; 0938 0939 /* Sun apex (where the sun goes) coordinates */ 0940 0941 dms asun(270.9592); // Right ascention: 18h 3m 50.2s [J2000] 0942 dms dsun(30.00467); // Declination: 30^o 0' 16.8'' [J2000] 0943 vsun = 20.; // [km/s] 0944 0945 asun.SinCos(sa, ca); 0946 dsun.SinCos(sd, cd); 0947 0948 /* We need an auxiliary SkyPoint since we need the 0949 * source referred to the J2000 equinox and we do not want to overwrite 0950 * the current values 0951 */ 0952 0953 SkyPoint aux; 0954 aux.set(RA0, Dec0); 0955 0956 aux.precessFromAnyEpoch(jd0, J2000L); 0957 0958 aux.ra().SinCos(sinRA, cosRA); 0959 aux.dec().SinCos(sinDec, cosDec); 0960 0961 /* Computation is done performing the scalar product of a unitary vector 0962 in the direction of the source with the vector velocity of Sun, both being in the 0963 LSR reference system: 0964 Vlsr = Vhel + Vsun.u_radial => 0965 Vlsr = Vhel + vsun(cos D cos A,cos D sen A,sen D).(cos d cos a,cos d sen a,sen d) 0966 Vhel = Vlsr - Vsun.u_radial 0967 */ 0968 0969 return vsun * (cd * cosDec * (cosRA * ca + sa * sinRA) + sd * sinDec); 0970 } 0971 0972 double SkyPoint::vHeliocentric(double vlsr, long double jd0) 0973 { 0974 return vlsr - vRSun(jd0); 0975 } 0976 0977 double SkyPoint::vHelioToVlsr(double vhelio, long double jd0) 0978 { 0979 return vhelio + vRSun(jd0); 0980 } 0981 0982 double SkyPoint::vREarth(long double jd0) 0983 { 0984 double sinRA, sinDec, cosRA, cosDec; 0985 0986 /* u_radial = unitary vector in the direction of the source 0987 Vlsr = Vhel + Vsun.u_radial 0988 = Vgeo + VEarth.u_radial + Vsun.u_radial => 0989 0990 Vgeo = (Vlsr -Vsun.u_radial) - VEarth.u_radial 0991 = Vhel - VEarth.u_radial 0992 = Vhel - (vx, vy, vz).(cos d cos a,cos d sen a,sen d) 0993 */ 0994 0995 /* We need an auxiliary SkyPoint since we need the 0996 * source referred to the J2000 equinox and we do not want to overwrite 0997 * the current values 0998 */ 0999 1000 SkyPoint aux(RA0, Dec0); 1001 1002 aux.precessFromAnyEpoch(jd0, J2000L); 1003 1004 aux.ra().SinCos(sinRA, cosRA); 1005 aux.dec().SinCos(sinDec, cosDec); 1006 1007 /* vEarth is referred to the J2000 equinox, hence we need that 1008 the source coordinates are also in the same reference system. 1009 */ 1010 1011 KSNumbers num(jd0); 1012 return num.vEarth(0) * cosDec * cosRA + num.vEarth(1) * cosDec * sinRA + num.vEarth(2) * sinDec; 1013 } 1014 1015 double SkyPoint::vGeocentric(double vhelio, long double jd0) 1016 { 1017 return vhelio - vREarth(jd0); 1018 } 1019 1020 double SkyPoint::vGeoToVHelio(double vgeo, long double jd0) 1021 { 1022 return vgeo + vREarth(jd0); 1023 } 1024 1025 double SkyPoint::vRSite(double vsite[3]) 1026 { 1027 double sinRA, sinDec, cosRA, cosDec; 1028 1029 RA.SinCos(sinRA, cosRA); 1030 Dec.SinCos(sinDec, cosDec); 1031 1032 return vsite[0] * cosDec * cosRA + vsite[1] * cosDec * sinRA + vsite[2] * sinDec; 1033 } 1034 1035 double SkyPoint::vTopoToVGeo(double vtopo, double vsite[3]) 1036 { 1037 return vtopo + vRSite(vsite); 1038 } 1039 1040 double SkyPoint::vTopocentric(double vgeo, double vsite[3]) 1041 { 1042 return vgeo - vRSite(vsite); 1043 } 1044 1045 bool SkyPoint::checkCircumpolar(const dms *gLat) const 1046 { 1047 return fabs(dec().Degrees()) > (90 - fabs(gLat->Degrees())); 1048 } 1049 1050 dms SkyPoint::altRefracted() const 1051 { 1052 return refract(Alt, Options::useRefraction()); 1053 } 1054 1055 void SkyPoint::setAltRefracted(dms alt_apparent) 1056 { 1057 setAlt(unrefract(alt_apparent, Options::useRefraction())); 1058 } 1059 1060 void SkyPoint::setAltRefracted(double alt_apparent) 1061 { 1062 setAlt(unrefract(alt_apparent, Options::useRefraction())); 1063 } 1064 1065 double SkyPoint::refractionCorr(double alt) 1066 { 1067 return 1.02 / tan(dms::DegToRad * (alt + 10.3 / (alt + 5.11))) / 60; 1068 } 1069 1070 double SkyPoint::refract(const double alt, bool conditional) 1071 { 1072 if (!conditional) 1073 { 1074 return alt; 1075 } 1076 static double corrCrit = SkyPoint::refractionCorr(SkyPoint::altCrit); 1077 1078 if (alt > SkyPoint::altCrit) 1079 return (alt + SkyPoint::refractionCorr(alt)); 1080 else 1081 return (alt + 1082 corrCrit * (alt + 90) / 1083 (SkyPoint::altCrit + 90)); // Linear extrapolation from corrCrit at altCrit to 0 at -90 degrees 1084 } 1085 1086 // Found uncorrected value by solving equation. This is OK since 1087 // unrefract is never called in loops with the potential exception of 1088 // slewing. 1089 // 1090 // Convergence is quite fast just a few iterations. 1091 double SkyPoint::unrefract(const double alt, bool conditional) 1092 { 1093 if (!conditional) 1094 { 1095 return alt; 1096 } 1097 double h0 = alt; 1098 double h1 = 1099 alt - 1100 (refract(h0) - 1101 h0); // It's probably okay to add h0 in refract() and subtract it here, since refract() is called way more frequently. 1102 1103 while (fabs(h1 - h0) > 1e-4) 1104 { 1105 h0 = h1; 1106 h1 = alt - (refract(h0) - h0); 1107 } 1108 return h1; 1109 } 1110 1111 dms SkyPoint::findAltitude(const SkyPoint *p, const KStarsDateTime &dt, const GeoLocation *geo, const double hour) 1112 { 1113 Q_ASSERT(p); 1114 if (!p) 1115 return dms(NaN::d); 1116 1117 // Jasem 2015-08-24 Using correct procedure to find altitude 1118 return SkyPoint::timeTransformed(p, dt, geo, hour).alt(); 1119 } 1120 1121 SkyPoint SkyPoint::timeTransformed(const SkyPoint *p, const KStarsDateTime &dt, const GeoLocation *geo, 1122 const double hour) 1123 { 1124 Q_ASSERT(p); 1125 if (!p) 1126 return SkyPoint(NaN::d, NaN::d); 1127 1128 // Jasem 2015-08-24 Using correct procedure to find altitude 1129 SkyPoint sp = *p; // make a copy 1130 KStarsDateTime targetDateTime = dt.addSecs(hour * 3600.0); 1131 dms LST = geo->GSTtoLST(targetDateTime.gst()); 1132 sp.EquatorialToHorizontal(&LST, geo->lat()); 1133 return sp; 1134 } 1135 1136 double SkyPoint::maxAlt(const dms &lat) const 1137 { 1138 double retval = (lat.Degrees() + 90. - dec().Degrees()); 1139 if (retval > 90.) 1140 retval = 180. - retval; 1141 return retval; 1142 } 1143 1144 double SkyPoint::minAlt(const dms &lat) const 1145 { 1146 double retval = (lat.Degrees() - 90. + dec().Degrees()); 1147 if (retval < -90.) 1148 retval = 180. + retval; 1149 return retval; 1150 } 1151 1152 dms SkyPoint::parallacticAngle(const CachingDms &LST, const CachingDms &lat) 1153 { 1154 // N.B. Technically, we could use the form 1155 // cos(angle) = (sin(φ) - sin(h) sin(δ))/(cos(h) cos(δ)) 1156 // where h = altitude, δ = declination, φ = latitude, 1157 // and then disambiguate the sign as 1158 // if (az().reduce() < 180°) angle = -angle; 1159 // However, acos(...) is inaccurate when cosine is nearly flat, i.e. near 0° and 180°. 1160 // It is therefore better to go through some extra pain to use atan2() 1161 1162 // Therefore we use the form shown in Jean Meeus' book (14.1) 1163 dms HA = LST - ra(); 1164 double tan_lat = lat.sin() / lat.cos(); 1165 double angle = atan2( // Measured CW on sky map (See Meeus' Fig on Pg 99) 1166 HA.sin(), 1167 tan_lat * dec().cos() - HA.cos() * dec().sin() 1168 ); 1169 return dms(angle / dms::DegToRad); 1170 } 1171 1172 1173 #ifndef KSTARS_LITE 1174 QDBusArgument &operator<<(QDBusArgument &argument, const SkyPoint &source) 1175 { 1176 argument.beginStructure(); 1177 argument << source.ra().Hours() << source.dec().Degrees(); 1178 argument.endStructure(); 1179 return argument; 1180 } 1181 1182 const QDBusArgument &operator>>(const QDBusArgument &argument, SkyPoint &dest) 1183 { 1184 double ra, dec; 1185 argument.beginStructure(); 1186 argument >> ra >> dec; 1187 argument.endStructure(); 1188 dest = SkyPoint(ra, dec); 1189 return argument; 1190 } 1191 #endif 1192