File indexing completed on 2024-04-28 11:26:07
0001 /* 0002 SPDX-FileCopyrightText: 2021 Akarsh Simha <akarsh@kde.org> 0003 0004 SPDX-License-Identifier: GPL-2.0-or-later 0005 */ 0006 0007 0008 #include "test_starobject.h" 0009 0010 #include "skyobjects/skypoint.h" 0011 #include "skyobjects/starobject.h" 0012 #include "ksnumbers.h" 0013 #include "time/kstarsdatetime.h" 0014 #include "auxiliary/dms.h" 0015 #include "nan.h" 0016 #include "Options.h" 0017 #include "config-kstars.h" 0018 0019 #ifdef HAVE_LIBERFA 0020 #include <erfa.h> 0021 #endif 0022 0023 #include <utility> 0024 0025 TestStarObject::TestStarObject() : QObject() 0026 { 0027 useRelativistic = Options::useRelativistic(); 0028 Options::setUseRelativistic(false); 0029 } 0030 0031 TestStarObject::~TestStarObject() 0032 { 0033 Options::setUseRelativistic(useRelativistic); 0034 } 0035 0036 // FIXME: Code duplication -- duplicated from TestSkyPoint::compare 0037 void TestStarObject::compare(QString msg, double ra1, double dec1, double ra2, double dec2, double err) 0038 { 0039 qDebug() << qPrintable(QString("%1 pos1 %2, %3 pos2 %4, %5 errors %6, %7 secs").arg(msg) 0040 .arg(ra1).arg(dec1).arg(ra2).arg(dec2) 0041 .arg((ra1 - ra2) * 3600.0, 6, 'f', 1).arg((dec1 - dec2) * 3600., 6, 'f', 1)); 0042 //QString str; 0043 //str.sprintf("%s pos1 %f, %f pos2 %f, %f errors %.1f, %.1f sec", msg.data(), ra1, dec1, ra2, dec2, (ra1 - ra2) *3600, (dec1 - dec2) * 3600 ); 0044 //qDebug() << str;// << msg << "SP " << ra1 << ", " << dec1 << " Sp0 " << ra2 << ", " << dec2 0045 //<< " errors " << ((ra1 - ra2) * 3600) << ", " << ((dec1 - dec2) * 3600) << " arcsec"; 0046 0047 double errRa = err / cos(dec1 * M_PI / 180.0); 0048 0049 QVERIFY2(fabs(ra1 - ra2) < errRa, qPrintable(QString("Ra %1, %2 error %3").arg(ra1).arg(ra2).arg(((ra1 - ra2) * 3600.0), 6, 0050 'f', 1))); 0051 QVERIFY2(fabs(dec1 - dec2) < err, qPrintable(QString("Dec %1, %2 error %3").arg(dec1).arg(dec2).arg((dec1 - dec2) * 3600., 0052 6, 'f', 1))); 0053 } 0054 0055 // FIXME: Code duplication -- duplicated from TestSkyPoint::compare 0056 void TestStarObject::compare(QString msg, double number1, double number2, double tolerance) 0057 { 0058 qDebug() << qPrintable(QString("%1 num1 %2 num2 %3 error %4").arg(msg) 0059 .arg(number1).arg(number2).arg(fabs(number1 - number2))); 0060 QVERIFY2(fabs(number1 - number2) < tolerance, 0061 qPrintable(QString("number1 %1, number2 %2; error %3").arg(number1).arg(number2).arg(fabs(number1 - number2)))); 0062 } 0063 0064 void TestStarObject::testUpdateCoordsStepByStep() 0065 { 0066 Options::setUseRelativistic(false); 0067 0068 /* 0069 * Run the worked example shown in Jean Meeus' "Astronomical 0070 * Algorithms" 2nd Edition, step-by-step 0071 * 0072 * The test-case is based on a combination of worked examples 21.b 0073 * and 23.a which concern the apparent place of the star theta 0074 * Persei on 2028 Nov 13.19 TD 0075 */ 0076 0077 // N.B. We can neglect the difference between TD and regular UTC 0078 // time for the purposes of this check. 0079 0080 // We check most results to arcsecond tolerance 0081 constexpr double subarcsecond_tolerance = 0.1 / 3600.0; 0082 0083 KStarsDateTime dt = KStarsDateTime::fromString("2028-11-13T04:33"); 0084 const KSNumbers num(dt.djd()); 0085 StarObject p; 0086 0087 /* Example 21.b */ 0088 0089 // ICRS coordinates 0090 CachingDms ra0 = dms::fromString("02:44:11.986", false); 0091 CachingDms dec0 = dms::fromString("+49:13:42.48", true); 0092 p.setRA0(ra0); 0093 p.setDec0(dec0); 0094 p.setProperMotion(0.03425 * 15.0 * 1000.0 * dec0.cos(), -0.0895 * 1000.0); // in mas/yr 0095 0096 // Mean equatorial coordinates 0097 p.getIndexCoords(&num, ra0, dec0); 0098 0099 // NOTE: We use a more accurate version of proper motion 0100 // correction than Meeus' example, whereby we will likely have 0101 // some difference 0102 compare( 0103 "Results of proper motion application", 0104 ra0.Degrees(), dec0.Degrees(), 0105 41.054063, 49.227750, 0106 subarcsecond_tolerance 0107 ); 0108 0109 // Set the proper motion corrected coordinates into the StarObject 0110 p.setRA0(ra0); 0111 p.setDec0(dec0); 0112 0113 qDebug() << "If the above test passed, our implementation of proper motion is likely correct!"; 0114 0115 // Compute precession 0116 p.precessFromAnyEpoch(J2000L, 0117 num.julianDay()); // Hmm... for some reason, SkyPoint::precess is protected, but this is public 0118 0119 // Compare the resulting mean equatorial coordinates with Meeus 0120 compare( 0121 "Mean equatorial coordinates", 0122 p.ra().Degrees(), p.dec().Degrees(), 0123 41.547214, 49.348483, 0124 subarcsecond_tolerance 0125 ); 0126 0127 /* Example 23.a */ 0128 dms ra = dms::fromString("02:46:11.331", false); // Meeus 23.a mean equatorial coordinates 0129 dms dec = dms::fromString("+49:20:54.54", true); // Meeus 23.a mean equatorial coordinates 0130 0131 // Proceed to compute the nutation and aberration, checking the inputs to the computation as well 0132 compare("Obliquity of the Ecliptic, ε in degrees", num.obliquity()->Degrees(), 23.436, 1e-3); // epsilon 0133 compare("Nutation in the longitude, Δψ in arcsec", num.dEcLong() * 3600.0, 14.861, 1.0); // Delta psi (arcsec) 0134 compare("Nutation in the obliquity, Δε in arcsec", num.dObliq() * 3600.0, 2.705, 1.0); // Delta epsilon (arcsec) 0135 compare("Eccentricity of earth orbit, e (unitless)", num.earthEccentricity(), 0.01669649, 1e-4); // e 0136 compare("True longitude of the Sun, L or ⊙ in degrees", num.sunTrueLongitude().reduce().Degrees(), 231.328, 0.1); // L 0137 // FIXME: Below has very large tolerance! 0138 compare("Longitude of earth perihelion, P or π in degrees", num.earthPerihelionLongitude().reduce().Degrees(), 103.434, 0139 0.5); // P / pi 0140 0141 p.nutate(&num); 0142 dms dra1(p.ra() - ra), ddec1(p.dec() - dec); // delta alpha 1, delta delta 1 0143 dms ra1(p.ra()), dec1(p.dec()); 0144 0145 p.aberrate(&num); 0146 dms dra2(p.ra() - ra1), ddec2(p.dec() - dec1); // delta alpha 2, delta delta 2 0147 0148 // N.B. We are checking Meeus' approximate method with a very 0149 // tight tolerance. If we "upgrade" to a better method such as 0150 // matrices for nutation or Ron-Vondrák method for aberration, 0151 // these checks will fail even though the result will be 0152 // correct! This is acknowledged in the following flag. 0153 0154 bool implementationMatchesMeeus = false; 0155 if (!SkyPoint::implementationIsLibnova) 0156 { 0157 implementationMatchesMeeus = true; 0158 } 0159 0160 const double meeusCheckTolerance = (implementationMatchesMeeus ? 0.01 : 0.5); 0161 0162 // N.B. We are doing this slightly differently from 0163 // Meeus since the corrections of aberration are applied on 0164 // the nutation-corrected coordinates and not on the mean 0165 // equatorial coordinates, but this apparently does not matter 0166 0167 0168 compare("RA Correction from Nutation in arcsec", dra1.Degrees() * 3600.0, 15.843, meeusCheckTolerance); // arcseconds 0169 compare("Dec Correction from Nutation in arcsec", ddec1.Degrees() * 3600.0, 6.218, meeusCheckTolerance); // arcseconds 0170 0171 compare("RA Correction from Aberration in arcsec", dra2.Degrees() * 3600.0, 30.045, meeusCheckTolerance); // arcseconds 0172 compare("Dec Correction from Aberration in arcsec", ddec2.Degrees() * 3600.0, 6.697, meeusCheckTolerance); // arcseconds 0173 0174 qDebug() << "If the above tests have passed, our implementation matches the approximate method given in Meeus to " << 0175 meeusCheckTolerance << " arcseconds"; 0176 0177 /* End-to-end result: Combination of Example 21.b and 23.a against StarObject::updateCoords */ 0178 ra0 = dms::fromString("02:44:11.986", false); 0179 dec0 = dms::fromString("+49:13:42.48", true); 0180 p.setRA0(ra0); 0181 p.setDec0(dec0); 0182 0183 // Note: Since Meeus directly adds the proper motion correction to 0184 // the RA, it must be without the cos(dec) factor, whereas KStars 0185 // expects it WITH the cos(dec) factor 0186 p.setProperMotion(0.03425 * 15.0 * 1000.0 * dec0.cos(), -0.0895 * 1000.0); // in mas/yr 0187 p.updateCoordsNow(&num); 0188 compare( 0189 "End-to-end computation of StarObject::updateCoords on Meeus Example 21.b + 23.a", 0190 p.ra().Degrees(), p.dec().Degrees(), 0191 dms::fromString("02:46:14.390", false).Degrees(), dms::fromString("+49:21:07.45", true).Degrees(), 0192 subarcsecond_tolerance 0193 ); 0194 } 0195 0196 0197 void TestStarObject::testUpdateCoords() 0198 { 0199 /* 0200 * End-to-end check on StarObject::updateCoords() against examples 0201 * pulled from various professional sources 0202 */ 0203 0204 struct TestCase 0205 { 0206 KStarsDateTime dt; 0207 dms RA0, Dec0; 0208 dms RA, Dec; 0209 double pmRa, pmDec; 0210 0211 TestCase(const QString &dtstr, const QString &icrs, const QString &apparent, double pmRa_ = 0., double pmDec_ = 0.) 0212 : dt(KStarsDateTime::fromString(dtstr)), pmRa(pmRa_), pmDec(pmDec_) 0213 { 0214 auto icrs_parts = icrs.split("|"); 0215 bool result = RA0.setFromString(icrs_parts.at(0), false); 0216 result = result && Dec0.setFromString(icrs_parts.at(1), true); 0217 0218 auto apparent_parts = apparent.split("|"); 0219 result = result && RA.setFromString(apparent_parts.at(0), false); 0220 result = result && Dec.setFromString(apparent_parts.at(1), true); 0221 if (!result) 0222 { 0223 qDebug() << "NOTE: TEST CASE IS BROKEN!"; 0224 QVERIFY(false); 0225 } 0226 } 0227 }; // FIXME: Move to TestStarObject::testUpdateCoords_data() 0228 0229 constexpr double subarcsecond_tolerance = 0.5 / 3600.0; 0230 0231 QList<TestCase> testCases; 0232 0233 // The following test cases are constructed by taking apparent 0234 // coordinates from the "Apparent Places of Fundamnetal Stars" 0235 // published by Astronomiches Rechen-Institut Heidelberg: 0236 // https://www.worldcat.org/title/apparent-places-of-fundamental-stars/oclc/1722620 0237 0238 // We query stars from the FK5 catalog (and assume FK5 ~ ICRS to 0239 // the accuracy promised by KStars) at the following interface to 0240 // obtain apparent positions: 0241 // https://wwwadd.zah.uni-heidelberg.de/datenbanken/ariapfs/query.php.en 0242 0243 // The ICRS coordinates and proper motions are obtained via SIMBAD 0244 0245 // N.B. For the test-cases after 2000, please note that we must 0246 // take the RA results from the "Equinox" method as that refers to 0247 // GCRS. The other ones seem to be for the CIRS which is a 0248 // different frame. 0249 0250 testCases 0251 // Before 2000 0252 << TestCase("1998-01-25T19:12", "03:43:14.90|-09:45:48.21", "03:43:09.58|-09:46:27.82", -93.16, 0253 743.64) // FK5 135 == HIP17378 0254 << TestCase("1998-01-25T00:00", "02:31:49.09|+89:15:50.79", "02:30:20.47|+89:15:33.43", 44.48, -11.85) // FK5 907 == Polaris 0255 0256 // After 2000 0257 << TestCase("2010-11-14T00:00", "03:24:19.37|+49:51:40.24", "03:25:09.64|+49:54:05.23", 23.75, -26.23) // FK5 120 == Mirfak 0258 << TestCase("2021-06-10T13:12", "06:23:57.11|-52:41:44.38", "06:24:23.09|-52:42:31.17", 19.9, 23.2) // FK5 245 == Canopus 0259 << TestCase("2021-01-13T19:26", "02:31:49.09|+89:15:50.79", "02:58:49.17|+89:21:23.23", 44.48, -11.85) // FK5 907 == Polaris 0260 0261 // High PM stars 0262 << TestCase("2021-12-16T14:38", "20 14 16.62|+15 11 51.37", "20 15 15.56|+15 15 54.17", 55.03, 0263 58.14) // FK5 1526 == rho Aql 0264 << TestCase("2021-12-16T13:41", "19 23 53.17|-40 36 57.37", "19 25 21.39|-40 34 32.46", 30.49, 0265 -119.21) // FK5 728 == alp Sgr 0266 << TestCase("2021-10-22T00:00", "02 07 10.40|+23 27 44.70", "02 08 24.64|+23 33 55.57", 188.55, 0267 -148.08) // FK5 74 == alp Ari 0268 0269 // South pole 0270 << TestCase("2021-11-30T16:48", "21 08 46.84|-88 57 23.41", "21 26 11.11|-88 52 16.42", 26.671, 5.612) // FK5 923 == sig Oct 0271 0272 // Near NEP (for high aberration) 0273 << TestCase("2021-03-30T06:43", "19 12 33.30|+67 39 41.54", "19 12 32.63|+67 41 30.58", 95.74, 91.92) // FK5 723 == del Dra 0274 << TestCase("2021-12-22T13:12", "19 12 33.30|+67 39 41.54", "19 12 29.58|+67 42 00.62", 95.74, 91.92) // FK5 723 == del Dra 0275 ; 0276 0277 // TODO: Add more test cases in the 1980s within FK5 reference frame 0278 0279 int i = 0; 0280 for (auto &testCase : testCases) 0281 { 0282 StarObject s {testCase.RA0, testCase.Dec0, 0.0, "", "", "K0", testCase.pmRa, testCase.pmDec, 0.0, false, false, 0}; 0283 KSNumbers num(testCase.dt.djd()); 0284 qDebug() << "Computing apparent position for (" << testCase.RA0.Degrees() << ", " << testCase.Dec0.Degrees() << ")"; 0285 s.updateCoordsNow(&num); 0286 compare( 0287 QString("Testcase %1").arg(i), s.ra().Degrees(), s.dec().Degrees(), testCase.RA.Degrees(), testCase.Dec.Degrees(), 0288 subarcsecond_tolerance 0289 ); 0290 i++; 0291 } 0292 0293 } 0294 0295 #ifdef HAVE_LIBERFA 0296 void TestStarObject::compareProperMotionAgainstErfa_data() 0297 { 0298 QTest::addColumn<QString>("ra"); 0299 QTest::addColumn<QString>("dec"); 0300 QTest::addColumn<double>("pmRA"); 0301 QTest::addColumn<double>("pmDec"); 0302 QTest::addColumn<double>("epoch"); 0303 0304 // A selection of real stars spread across the sky 0305 QTest::newRow("Star 1") << "22 10 11.9852752" << "+06 11 52.307750" << 282.18 << 30.46 << 2028.3; 0306 QTest::newRow("Star 2") << "02 44 11.9870420" << "+49 13 42.411120" << 334.66 << -89.99 << 2074.8; 0307 QTest::newRow("Star 3") << "00 01 35.7015845" << "-77 03 56.609236" << -57.30 << -177.06 << 1934.2; 0308 QTest::newRow("Star 4") << "06 23 57.10988" << "-52 41 44.3810" << 19.93 << 23.24 << 1901.5; 0309 QTest::newRow("Star 5") << "14 51 38.3028905" << "-43 34 31.296546" << -25.20 << -27.13 << 1995.2; 0310 QTest::newRow("Star 6") << "01 09 43.92388" << "+35 37 14.0075" << 175.90 << -112.20 << 2140.3; 0311 QTest::newRow("Star 7") << "04 16 01.588231" << "-51 29 11.919063" << 99.463 << 183.353 << 2093.8; 0312 QTest::newRow("Star 8") << "14 39 29.71993" << "-60 49 55.9990" << -3608. << 686.0 << 2023.0; 0313 } 0314 0315 void TestStarObject::compareProperMotionAgainstErfa() 0316 { 0317 QFETCH(QString, ra); 0318 QFETCH(QString, dec); 0319 QFETCH(double, pmRA); 0320 QFETCH(double, pmDec); 0321 QFETCH(double, epoch); 0322 0323 dms _ra, _dec; 0324 bool result = _ra.setFromString(ra, false); 0325 result = result && _dec.setFromString(dec, true); 0326 Q_ASSERT(result); 0327 0328 auto startPoint = std::make_pair(_ra, _dec); 0329 long double jdf = KStarsDateTime::epochToJd(epoch); 0330 0331 std::pair<dms, dms> erfaResult; 0332 { 0333 /* 0334 * startPoint : pair of dms: (ra, dec) 0335 * pmRA : mas/yr, inclusive of cos(dec) factor 0336 * pmDec : mas/yr 0337 * jdf : final (target) epoch 0338 */ 0339 0340 constexpr double masToRad = dms::DegToRad / (3600.0 * 1000.0); 0341 double ra1 = startPoint.first.reduce().radians(); 0342 double dec1 = startPoint.second.radians(); 0343 double pmr1 = pmRA * masToRad / cos(dec1); 0344 double pmd1 = pmDec * masToRad; 0345 0346 double resRa {NaN::d}, resDec {NaN::d}, resPmRa {NaN::d}, resPmDec {NaN::d}, resPx, resRv; 0347 /* 0348 * Note: eraStarpm does not like a zero parallax. It replaces 0349 * it with a 1e-7 parallax (which would be fine), but then 0350 * concludes that the spatial velocity corresponding to the 0351 * provided proper motions is relativistic as a result. So a 0352 * typical parallax needs to be provided that is reasonable 0353 * for most high-PM stars. I've chosen to set that to 0.1, 0354 * which should correspond to 10pc ~ 30ly. This distance `r` 0355 * will divide out as long as the radial velocity is zero: in 0356 * other words, it does not matter what size of sphere we pick 0357 * 0358 * Another note: The documentation reads "The RA proper motion 0359 * is in terms of coordinate angle, not true angle. If the 0360 * catalog uses arcseconds for both RA and Dec proper motions, 0361 * the RA proper motion will need to be divided by cos(Dec) 0362 * before use." 0363 */ 0364 int result = eraStarpm( 0365 ra1, dec1, pmr1, pmd1, 0366 0.1, 0., // No parallax, radial velocity info 0367 J2000, 0, // Starting epoch is J2000 0368 double(jdf), 0, // Ending epoch is jdf 0369 &resRa, &resDec, &resPmRa, &resPmDec, 0370 &resPx, &resRv 0371 ); // return value irrelevant to us 0372 0373 if (result > 0) 0374 { 0375 qDebug() << "ERFA reports status " << result << ", correction (degrees) dRA = " << (resRa - ra1) / dms::DegToRad << ", dDec = " << (resDec - dec1) / dms::DegToRad; 0376 } 0377 0378 erfaResult = std::make_pair(dms { resRa / dms::DegToRad }, dms { resDec / dms::DegToRad }); 0379 } 0380 0381 CachingDms ra_f, dec_f; 0382 KSNumbers num(jdf); 0383 StarObject s {startPoint.first, startPoint.second, 0.0, "", "", "K0", pmRA, pmDec, 0.0, false, false, 0}; 0384 s.getIndexCoords(&num, ra_f, dec_f); // dms version 0385 double ra_fd, dec_fd; 0386 s.getIndexCoords(&num, &ra_fd, &dec_fd); // Degrees version 0387 0388 compare("KStars internal consistency: getIndexCoords degree vs dms implementation", 0389 ra_f.reduce().Degrees(), dec_f.Degrees(), 0390 dms(ra_fd).reduce().Degrees(), dec_fd, 0391 1e-8); // Same algorithm, should match exactly! 0392 0393 constexpr double arcsecond_tolerance = 1.0/3600.0; 0394 compare( 0395 "Proper Motion KStars vs ERFA", 0396 ra_f.reduce().Degrees(), dec_f.Degrees(), 0397 erfaResult.first.reduce().Degrees(), erfaResult.second.Degrees(), 0398 arcsecond_tolerance 0399 ); 0400 0401 } 0402 #endif 0403 0404 QTEST_GUILESS_MAIN(TestStarObject)