File indexing completed on 2024-06-09 15:20:31

0001 /*
0002     SPDX-FileCopyrightText: 2021 Hy Murveit <hy@murveit.com>
0003 
0004     SPDX-License-Identifier: GPL-2.0-or-later
0005 */
0006 
0007 /* Project Includes */
0008 #include "test_polaralign.h"
0009 #include "kstarsdata.h"
0010 #include "ksnumbers.h"
0011 #include "time/kstarsdatetime.h"
0012 #include "auxiliary/dms.h"
0013 #include "Options.h"
0014 #include <libnova/libnova.h>
0015 #include "fitsviewer/fitscommon.h"
0016 #include "fitsviewer/fitsdata.h"
0017 #include "rotations.h"
0018 
0019 using Rotations::V3;
0020 
0021 // Solver's solution. RA, DEC, & orientation in degrees, pixScale in arc-seconds/pixel,
0022 // and the sample's time in UTC.
0023 struct Solution
0024 {
0025     double ra, dec, orientation, pixScale;
0026     int year, month, day, hour, minute, second;
0027 };
0028 
0029 // A set of 3 solutions for the polar alignment.
0030 struct PaaData
0031 {
0032     Solution s1, s2, s3;
0033     int x, y, correctedX, correctedY;
0034 };
0035 
0036 constexpr int IMAGE_WIDTH = 4656;
0037 constexpr int IMAGE_HEIGHT = 3520;
0038 
0039 void loadDummyFits(QSharedPointer<FITSData> &image, const KStarsDateTime &time,
0040                    double ra, double dec, double orientation, double pixScale, bool eastToTheRight)
0041 {
0042     image.reset(new FITSData(FITS_NORMAL));
0043 
0044     // Borrow one of the other fits files in the test suite.
0045     // We don't use the contents of the image, but the image's width and height.
0046     // We do set the wcs data.
0047     QFuture<bool> worker = image->loadFromFile("ngc4535-autofocus1.fits");
0048     QTRY_VERIFY_WITH_TIMEOUT(worker.isFinished(), 60000);
0049     QVERIFY(worker.result());
0050     auto stats = image->getStatistics();
0051     stats.width = IMAGE_WIDTH;
0052     stats.height = IMAGE_HEIGHT;
0053     image->restoreStatistics(stats);
0054     image->setDateTime(time);
0055     image->injectWCS(orientation, ra, dec, pixScale, eastToTheRight);
0056     QVERIFY(image->loadWCS());
0057 }
0058 
0059 void setupData(Solution s, QSharedPointer<FITSData> &image, bool eastToTheRight)
0060 {
0061     KStarsDateTime time;
0062     time.setDate(QDate(s.year, s.month, s.day));
0063     time.setTime(QTime(s.hour, s.minute, s.second));
0064     time.setTimeSpec(Qt::UTC);
0065     loadDummyFits(image, time, s.ra, s.dec, s.orientation, s.pixScale, eastToTheRight);
0066 }
0067 
0068 void setupData(const PaaData &data, int sampleNum, QSharedPointer<FITSData> &image, bool eastToTheRight)
0069 {
0070     const Solution &s = sampleNum == 0 ? data.s1 : sampleNum == 1 ? data.s2 : data.s3;
0071     setupData(s, image, eastToTheRight);
0072 }
0073 
0074 TestPolarAlign::TestPolarAlign() : QObject()
0075 {
0076     Options::setUseRelativistic(false);
0077 }
0078 
0079 TestPolarAlign::~TestPolarAlign()
0080 {
0081 }
0082 
0083 void TestPolarAlign::compare(double a, double e, QString msg, int line, double tolerance)
0084 {
0085     QVERIFY2(std::fabs(a - e) < tolerance,
0086              qPrintable(QString("Line %1: %2: actual %3, expected %4 error %5 arc-seconds")
0087                         .arg(line).arg(msg).arg(a).arg(e)
0088                         .arg(((a - e) * 3600.0), 3, 'f', 1)));
0089 }
0090 
0091 bool TestPolarAlign::compare(const QPointF &point, double x, double y, double tolerance)
0092 {
0093     return ((std::fabs(point.x() - x) < tolerance) &&
0094             (std::fabs(point.y() - y) < tolerance));
0095 }
0096 
0097 namespace
0098 {
0099 void runPAA(const GeoLocation &geo, const PaaData &data, bool eastToTheRight = true)
0100 {
0101     constexpr bool canSkipPixelError = false;
0102 
0103     PolarAlign polarAlign(&geo);
0104 
0105     QSharedPointer<FITSData> image;
0106 
0107     setupData(data, 0, image, eastToTheRight);
0108     QVERIFY(polarAlign.addPoint(image));
0109 
0110     setupData(data, 1, image, eastToTheRight);
0111     QVERIFY(polarAlign.addPoint(image));
0112 
0113     setupData(data, 2, image, eastToTheRight);
0114     QVERIFY(polarAlign.addPoint(image));
0115 
0116     QVERIFY(polarAlign.findAxis());
0117 
0118     double axisAz, axisAlt;
0119     polarAlign.getAxis(&axisAz, &axisAlt);
0120 
0121     double azimuthError, altitudeError;
0122     polarAlign.calculateAzAltError(&azimuthError, &altitudeError);
0123     polarAlign.setMaxPixelSearchRange(std::hypot(azimuthError, altitudeError) + 1);
0124 
0125     QPointF corrected;
0126     if (data.x < 0 && data.y < 0)
0127     {
0128         QVERIFY(polarAlign.findCorrectedPixel(
0129                     image, QPointF(image->width() / 2, image->height() / 2), &corrected));
0130     }
0131     else
0132     {
0133         QVERIFY(polarAlign.findCorrectedPixel(
0134                     image, QPointF(data.x, data.y), &corrected));
0135         double azE, altE;
0136 
0137         // Just fix the altitude, not the azimuth, for the pixelError tests below.
0138         QPointF altCorrected;
0139         QVERIFY(polarAlign.findCorrectedPixel(
0140                     image, QPointF(data.x, data.y), &altCorrected, true));
0141 
0142         // Below pixelError tests, test our ability to estimate how much error is left to be corrected
0143         // as the user progresses along the polar alignment procedure.
0144         // We use pretty generous error margins below since pixelError is an approximation.
0145 
0146         // Test the entire path, start to end. We should see the full azimuth and altitude errors.
0147         QVERIFY(polarAlign.pixelError(image, QPointF(data.x, data.y), corrected, &azE, &altE));
0148         QVERIFY(fabs(azE - azimuthError) < .02);
0149         QVERIFY(fabs(altE - altitudeError) < .01);
0150 
0151         // Simulate that the user has started correcting, and is halfway on the alt-only segment,
0152         // there should be 1/2 the alt error left, and all of the az error.
0153         if (polarAlign.pixelError(image, QPointF((altCorrected.x() + data.x) / 2, (altCorrected.y() + data.y) / 2), corrected,
0154                                   &azE, &altE))
0155         {
0156             QVERIFY(fabs(azE - azimuthError) < .01);
0157             QVERIFY(fabs(altE - altitudeError / 2.0) < .01);
0158         }
0159         else QVERIFY(canSkipPixelError);
0160 
0161         // Simulate that the full alt correction path has been completed.
0162         // The azimuth error should not be fixed, but alt should be fully corrected.
0163         if (polarAlign.pixelError(image, altCorrected, corrected, &azE, &altE))
0164         {
0165             QVERIFY(fabs(azE - azimuthError) < .01);
0166             QVERIFY(fabs(altE) < .01);
0167         }
0168         else QVERIFY(canSkipPixelError);
0169 
0170         // Now simulate that the user has gone further, halfway through on the az path.
0171         // There should be 1/2 the az error left but none of the alt error should remain.
0172         if (polarAlign.pixelError(image, QPointF((corrected.x() + altCorrected.x()) / 2,
0173                                   (corrected.y() + altCorrected.y()) / 2), corrected, &azE, &altE))
0174         {
0175             QVERIFY(fabs(azE - azimuthError / 2) < .01);
0176             QVERIFY(fabs(altE) < .01);
0177         }
0178         else QVERIFY(canSkipPixelError);
0179 
0180         // At the end there should be no error left.
0181         if (polarAlign.pixelError(image, corrected, corrected, &azE, &altE))
0182         {
0183             QVERIFY(fabs(azE) < .01);
0184             QVERIFY(fabs(altE) < .01);
0185         }
0186         else QVERIFY(canSkipPixelError);
0187 
0188         // Test the alt-only path.
0189         // Using that correction path, there should be no Azimuth error, but alt error remains.
0190         if (polarAlign.pixelError(image, QPointF(data.x, data.y), altCorrected, &azE, &altE))
0191         {
0192             QVERIFY(fabs(azE) < .01);
0193             QVERIFY(fabs(altE - altitudeError) < .01);
0194         }
0195         else QVERIFY(canSkipPixelError);
0196 
0197     }
0198     // Some approximations in the time (no sub-second decimals), and the
0199     // findAltAz search resolution introduce a pixel or two of error.
0200     // This is way below the noise in polar alignment.
0201     QVERIFY(data.correctedX < 0 || abs(data.correctedX - corrected.x()) < 3);
0202     QVERIFY(data.correctedY < 0 || abs(data.correctedY - corrected.y()) < 3);
0203 }
0204 }  // namespace
0205 
0206 void TestPolarAlign::testRunPAA()
0207 {
0208     const GeoLocation siliconValley(dms(-122, 10), dms(37, 26, 30));
0209     const GeoLocation dallas(dms(-96, 49), dms(32, 46, 45));
0210     const GeoLocation kuwait(dms(47, 59), dms(29, 22, 43));
0211 
0212     // from log on 12/29 log_19-14-48.txt
0213     // First 14m az error, low alt error
0214     runPAA(siliconValley,
0215     {
0216         // ra0        dec0   orientation pixscale  date-time in UTC
0217         { 21.74338, 89.83996,  78.36826, 1.32548, 2020, 12, 30, 03, 16, 8},
0218         { 10.17858, 89.82383,  95.80819, 1.32543, 2020, 12, 30, 03, 16, 46},
0219         {358.41548, 89.82049, 113.44157, 1.32468, 2020, 12, 30, 03, 17, 18},
0220         // mount axis x,y  pole center in x,y
0221         1880, 1565, 1858, 1559
0222     });
0223 
0224     // Corrected the above, now very low error (0 0' 11")
0225     runPAA(siliconValley,
0226     {
0227         { 74.23824, 89.75702, 127.82348, 1.32612, 2020, 12, 30, 3, 28, 50},
0228         { 66.19830, 89.77333, 148.11177, 1.32484, 2020, 12, 30, 3, 29, 22},
0229         { 59.87800, 89.80076, 170.45317, 1.32384, 2020, 12, 30, 3, 29, 59},
0230         2238, 1225, 1842, 1576
0231     });
0232 
0233     // Manually increased the alt error, now 0 23' 50'
0234     runPAA(siliconValley,
0235     {
0236         { 46.61670, 89.41158,  98.86135, 1.32586,  2020, 12, 30, 3, 34, 02},
0237         { 43.13352, 89.41121, 123.87848, 1.32566, 2020, 12, 30, 3, 34, 39},
0238         { 40.04437, 89.42768, 149.79166, 1.32392, 2020, 12, 30, 3, 35, 11},
0239         1544, 415, 1842, 1585
0240     });
0241 
0242     // Fixed it again.
0243     runPAA(siliconValley,
0244     {
0245         { 75.31905, 89.76408, 126.18540, 1.32508, 2020, 12, 30, 3, 39, 33},
0246         { 67.12354, 89.77885, 145.52730, 1.32492, 2020, 12, 30, 3, 40, 04},
0247         { 60.22412, 89.80675, 168.40166, 1.32473, 2020, 12, 30, 3, 40, 37},
0248         2222, 1245, 1842, 1597
0249     });
0250 
0251     // From 12/31
0252     // Comments: the PA error with double precision, float precision, and original scheme.
0253     // Commented solution is the float solution.
0254 
0255     // double: 10'28" float: 11'45" orig: 9'42"
0256     runPAA(siliconValley,
0257     {
0258         { 63.69466, 89.77876, 124.27465, 1.32519, 2020, 12, 31, 18, 52, 42},
0259         { 54.75964, 89.79547, 143.06160, 1.32367, 2020, 12, 31, 18, 53, 15},
0260         { 47.33568, 89.82357, 165.12844, 1.32267, 2020, 12, 31, 18, 53, 52},
0261         2204, 1295, 1851, 1550
0262     });
0263 
0264     // double: 0'40", float: 0'00", orig: 0'41"
0265     runPAA(siliconValley,
0266     {
0267         { 21.22696, 89.82033, 79.74928, 1.32493, 2020, 12, 31, 19, 00, 52},
0268         { 10.14439, 89.80891, 97.13185, 1.32523, 2020, 12, 31, 19, 01, 25},
0269         { 359.05594, 89.81023, 114.39229, 1.32437, 2020, 12, 31, 19, 01, 58},
0270         1858, 1546, 1847, 1574
0271     });
0272 
0273     // double: 11'59", float: 11'01", orig: 11'08"
0274     runPAA(siliconValley,
0275     {
0276         { 266.87080, 89.98377, -35.40868, 1.33053, 2020, 12, 31, 19, 06, 46},
0277         { 293.43530, 89.94730, 19.32200, 1.32335, 2020, 12, 31, 19, 07, 24},
0278         { 285.55085, 89.91166, 40.10587, 1.32260, 2020, 12, 31, 19, 07, 56},
0279         2173, 1943, 1842, 1574
0280     });
0281 
0282     // double: 1'14", float: 0'00", orig: 0'57"
0283     runPAA(siliconValley,
0284     {
0285         { 22.74402, 89.82040, 78.32075, 1.32462, 2020, 12, 31, 19, 12, 34},
0286         { 11.69887, 89.80704, 95.88488, 1.32375, 2020, 12, 31, 19, 13, 11},
0287         { 0.95646, 89.80782, 113.01201, 1.32432, 2020, 12, 31, 19, 13, 43},
0288         1847, 1555, 1848, 1598
0289     });
0290 
0291     // double: 32'28", float: 38'37", orig: 30'19"
0292     runPAA(siliconValley,
0293     {
0294         { 317.02018, 89.44380, 10.85411, 1.32380, 2020, 12, 31, 19, 18, 00},
0295         { 316.55330, 89.40583, 38.75090, 1.32350, 2020, 12, 31, 19, 18, 33},
0296         { 314.63678, 89.37679, 64.65160, 1.32335, 2020, 12, 31, 19, 19, 06},
0297         795, 2485, 1826, 1595
0298     });
0299 
0300     // double: 1'34", float: 1'30", orig: 1'26"
0301     runPAA(siliconValley,
0302     {
0303         { 27.96856, 89.82500, 80.98310, 1.32435, 2020, 12, 31, 19, 22, 48},
0304         { 16.47574, 89.81367, 97.74878, 1.32394, 2020, 12, 31, 19, 23, 25},
0305         { 5.50816, 89.81455, 114.24102, 1.32390, 2020, 12, 31, 19, 23, 58},
0306         1868, 1552, 1831, 1604
0307     });
0308 
0309     // double: 0'25", float: 0'00", orig: 0'25"
0310     runPAA(siliconValley,
0311     {
0312         { 22.29262, 89.82396, 73.68740, 1.32484, 2020, 12, 31, 19, 29, 17},
0313         { 11.20186, 89.80799, 91.89482, 1.32477, 2020, 12, 31, 19, 29, 55},
0314         { 0.26659, 89.80523, 109.34529, 1.32489, 2020, 12, 31, 19, 30, 27},
0315         1828, 1584, 1831, 1602
0316     });
0317 
0318     // double: 17'35", float: 16'46", orig: 16'21"
0319     runPAA(siliconValley,
0320     {
0321         { 45.20803, 89.58706, 95.79552, 1.32496, 2020, 12, 31, 19, 32, 40},
0322         { 39.89580, 89.58748, 119.13221, 1.32328, 2020, 12, 31, 19, 33, 13},
0323         { 35.18193, 89.60542, 144.28407, 1.32255, 2020, 12, 31, 19, 33, 51},
0324         1700, 887, 1833, 1606
0325     });
0326 
0327     // double: 0'15", float: 0'00", orig: 0'04"
0328     runPAA(siliconValley,
0329     {
0330         { 22.79496, 89.83001, 71.18473, 1.32378, 2020, 12, 31, 19, 41, 16},
0331         { 11.65388, 89.81182, 89.34180, 1.32420, 2020, 12, 31, 19, 41, 54},
0332         { 1.02046, 89.80613, 106.01159, 1.32431, 2020, 12, 31, 19, 42, 26},
0333         1821, 1614, 1818, 1615
0334     });
0335 
0336     // Log from 1/5/20201
0337     runPAA(siliconValley,
0338     {
0339         { 33.77699, 89.85172, 55.15250, 1.32489, 2021, 01, 06, 05,  9, 23 },
0340         { 23.30335, 89.82440, 74.05880, 1.32592, 2021, 01, 06, 05,  9, 50 },
0341         { 12.72288, 89.81051, 91.35879, 1.32540, 2021, 01, 06, 05, 10, 15 },
0342         2328, 1760, 2331, 1783
0343     });
0344 
0345     // Note, far off pole
0346     runPAA(siliconValley,
0347     {
0348         { 50.76609, 49.91458, -8.37603, 1.32461, 2021, 01, 06, 05, 13, 32 },
0349         { 23.24014, 49.88528, -8.28620, 1.32496, 2021, 01, 06, 05, 13, 57 },
0350         { 354.47536, 49.87901, -8.20577, 1.32499, 2021, 01, 06, 05, 14, 24 },
0351         2328, 1760, 2320, 1795
0352     });
0353 
0354     // Note, far off pole
0355     runPAA(siliconValley,
0356     {
0357         { 50.73028, 49.49321, -8.27081, 1.32448, 2021, 01, 06, 05, 18, 31 },
0358         { 21.45768, 49.52983, -8.00221, 1.32467, 2021, 01, 06, 05, 18, 58 },
0359         { 353.23623, 49.64241, -7.79150, 1.32521, 2021, 01, 06, 05, 19, 24 },
0360         3007, 1256, 3516, 951
0361     });
0362 
0363     runPAA(siliconValley,
0364     {
0365         { 33.93431, 89.84705, 51.18593, 1.32465, 2021, 01, 06, 05, 25, 48 },
0366         { 25.03596, 89.82113, 69.85483, 1.32561, 2021, 01, 06, 05, 26, 13 },
0367         { 14.68785, 89.80521, 88.25234, 1.32571, 2021, 01, 06, 05, 26, 40 },
0368         2328, 1760, 2350, 1785
0369     });
0370 
0371     // Note, far off pole
0372     runPAA(siliconValley,
0373     {
0374         { 50.38595, 49.53670, -8.15094, 1.32484, 2021, 01, 06, 05, 30, 55 },
0375         { 22.40326, 49.57133, -8.12608, 1.32515, 2021, 01, 06, 05, 31, 21 },
0376         { 354.02099, 49.60593, -8.13631, 1.32476, 2021, 01, 06, 05, 31, 47 },
0377         2446, 2088, 2392, 2522
0378     });
0379 
0380 
0381     /*
0382     // pixelError has issues with this one. Letting it slide for now.
0383     runPAA(siliconValley,
0384     {
0385         { 80.96890, 11.74259, 171.42773, 1.32514, 2021, 01, 06, 05, 45, 05 },
0386         { 110.26384, 11.81600, 171.56395, 1.32316, 2021, 01, 06, 05, 45, 31 },
0387         { 138.33895, 11.78703, 171.85951, 1.32656, 2021, 01, 06, 05, 45, 57 },
0388         1087, 861, 454, 130
0389     });
0390     */
0391 
0392     runPAA(siliconValley,
0393     {
0394         { 38.40316, 89.32378, 49.26512, 1.32596, 2021, 01, 06, 05, 50, 17 },
0395         { 36.29482, 89.29488, 75.01016, 1.32645, 2021, 01, 06, 05, 50, 42 },
0396         { 33.52382, 89.28142, 100.14838, 1.32604, 2021, 01, 06, 05, 51, 8 },
0397         2328, 1760, 3746, 2172
0398     });
0399 
0400     // Note, far off pole
0401     runPAA(siliconValley,
0402     {
0403         { 51.11166, 48.88918, -8.54258, 1.32508, 2021, 01, 06, 05, 52, 43 },
0404         { 23.14732, 48.86777, -8.06648, 1.32514, 2021, 01, 06, 05, 53, 9 },
0405         { 355.71059, 48.98942, -7.64124, 1.32527, 2021, 01, 06, 05, 53, 35 },
0406         2205, 2214, 3203, 1043
0407     });
0408 
0409     // Note, far off pole
0410     runPAA(siliconValley,
0411     {
0412         { 53.79635, 49.62699, -8.35495, 1.32678, 2021, 01, 06, 06, 03, 39 },
0413         { 25.42441, 49.58901, -8.32128, 1.32477, 2021, 01, 06, 06, 04, 05 },
0414         { 357.52396, 49.57574, -8.24039, 1.32470, 2021, 01, 06, 06, 04, 30 },
0415         895, 1336, 882, 1347
0416     });
0417 
0418     // Jasem3
0419     runPAA(kuwait,
0420     {
0421         { 102.65020, 39.69621, 138.49586, 1.09119, 2021, 01, 11, 17, 25, 10 },
0422         { 86.97538, 39.64818, 138.52258, 1.09179, 2021, 01, 11, 17, 25, 29 },
0423         { 74.74086, 39.61638, 138.54977, 1.09189, 2021, 01, 11, 17, 25, 47 },
0424         2328, 1760, 2613, 2003
0425     });
0426     // Jasem2
0427     runPAA(kuwait,
0428     {
0429         { 75.36752, 67.31292, 138.39209, 1.09149, 2021, 01, 11, 17, 21, 12 },
0430         { 88.67953, 67.34817, 138.41916, 1.09133, 2021, 01, 11, 17, 21, 30 },
0431         { 101.85324, 67.38325, 138.31124, 1.09222, 2021, 01, 11, 17, 21, 48 },
0432         2328, 1760, 2369, 1689
0433     });
0434     // Jasem1
0435     runPAA(kuwait,
0436     {
0437         { 103.47913, 67.38804, 138.36922, 1.09212, 2021, 01, 11, 17, 19, 57 },
0438         { 89.15202, 67.34897, 138.37273, 1.09211, 2021, 01, 11, 17, 20, 18 },
0439         { 75.36751, 67.31304, 138.37937, 1.09212, 2021, 01, 11, 17, 20, 36 },
0440         2328, 1760, 2453, 1796
0441     });
0442 
0443     // Jo-nightly-1
0444     runPAA(dallas,
0445     {
0446         { 12.81781, 89.07239, -8.65391, 9.49554, 2021, 01, 12, 0, 27, 23 },
0447         { 348.85538, 89.03890, -3.73479, 9.49479, 2021, 01, 12, 0, 27, 43 },
0448         { 326.22608, 89.04224, 1.19850, 9.49679, 2021, 01, 12, 0, 28, 03 },
0449         2328, 1760, 2319, 1730
0450     });
0451 
0452     // Jo-nightly-2
0453     runPAA(dallas,
0454     {
0455         { 18.65836, 89.05113, -4.20730, 9.49647, 2021, 01, 12, 0, 33, 19 },
0456         { 40.85876, 89.08147, -10.33458, 9.49571, 2021, 01, 12, 0, 33, 40 },
0457         { 65.37414, 89.14252, -13.67848, 9.49793, 2021, 01, 12, 0, 33, 59 },
0458         2328, 1760, 2310, 1748
0459     });
0460 
0461     // Jo-nightly-3
0462     runPAA(dallas,
0463     {
0464         { 19.27912, 89.05092, -4.02575, 9.49570, 2021, 01, 12, 0, 34, 57 },
0465         { 356.21018, 89.04905, 0.67433, 9.49619, 2021, 01, 12, 0, 35, 16 },
0466         { 331.16747, 89.08275, 5.35544, 9.49534, 2021, 01, 12, 0, 35, 37 },
0467         2328, 1760, 2343, 1749
0468     });
0469 
0470     // Jo-nightly-4
0471     runPAA(dallas,
0472     {
0473         { 21.31811, 89.09443, -2.94373, 9.49555, 2021, 01, 12, 0, 38, 59 },
0474         { 45.33237, 89.11518, -8.11396, 9.49733, 2021, 01, 12, 0, 39, 25 },
0475         { 71.91959, 89.16055, -10.72869, 9.49602, 2021, 01, 12, 0, 39, 45 },
0476         2328, 1760, 2329, 1751
0477     });
0478 
0479     // Jo-nightly-5
0480     runPAA(dallas,
0481     {
0482         { 21.89789, 89.09423, -2.73671, 9.49587, 2021, 01, 12, 0, 40, 24 },
0483         { 357.18691, 89.09645, 0.78761, 9.49571, 2021, 01, 12, 0, 40, 46 },
0484         { 330.78136, 89.12673, 4.15432, 9.49656, 2021, 01, 12, 0, 41, 06 },
0485         2328, 1760, 2336, 1764
0486     });
0487 
0488     // Jo-nightly-6
0489     runPAA(dallas,
0490     {
0491         { 21.90875, 89.09512, -3.01211, 9.49640, 2021, 01, 12, 0, 41, 38 },
0492         { 46.14103, 89.11651, -8.15365, 9.49657, 2021, 01, 12, 0, 41, 58 },
0493         { 72.50313, 89.16190, -10.72389, 9.49792, 2021, 01, 12, 0, 42, 17 },
0494         2328, 1760, 2329, 1750
0495     });
0496 
0497     // Jo orig-1
0498     runPAA(dallas,
0499     {
0500         { 22.67449, 89.09498, -2.83722, 9.49689, 2021, 01, 12, 0, 43, 54 },
0501         { 358.94560, 89.09623, 0.54390, 9.49583, 2021, 01, 12, 0, 44, 13 },
0502         { 331.80520, 89.12571, 4.01328, 9.49679, 2021, 01, 12, 0, 44, 34 },
0503         2328, 1760, 2336, 1765
0504     });
0505 
0506     // reversed
0507     runPAA(dallas,
0508     {
0509         { 331.80520, 89.12571, 4.01328, 9.49679, 2021, 01, 12, 0, 44, 34 },
0510         { 358.94560, 89.09623, 0.54390, 9.49583, 2021, 01, 12, 0, 44, 13 },
0511         { 22.67449, 89.09498, -2.83722, 9.49689, 2021, 01, 12, 0, 43, 54 },
0512         2328, 1760, 2336, 1756
0513     });
0514 
0515     // Jo orig-2
0516     runPAA(dallas,
0517     {
0518         { 23.51969, 89.09648, -3.24191, 9.49659, 2021, 01, 12, 0, 49, 15 },
0519         { 48.00993, 89.11943, -8.37887, 9.49623, 2021, 01, 12, 0, 49, 37 },
0520         { 75.23851, 89.16742, -10.84337, 9.49660, 2021, 01, 12, 0, 49, 59 },
0521         2328, 1760, 2329, 1750
0522     });
0523 
0524     // reversed
0525     runPAA(dallas,
0526     {
0527         { 75.23851, 89.16742, -10.84337, 9.49660, 2021, 01, 12, 0, 49, 59 },
0528         { 48.00993, 89.11943, -8.37887, 9.49623, 2021, 01, 12, 0, 49, 37 },
0529         { 23.51969, 89.09648, -3.24191, 9.49659, 2021, 01, 12, 0, 49, 15 },
0530         2328, 1760, 2337, 1756
0531     });
0532 
0533     runPAA(siliconValley,
0534     {
0535         // ra0       dec0   orientation pixscale  date-time in UTC
0536         { 47.31826, 1.94192, -84.57643, 1.32459, 2021, 01, 26, 3, 15, 28},
0537         { 18.35969, 2.08571, -84.37967, 1.32451, 2021, 01, 26, 3, 15, 54},
0538         {350.43339, 2.27862, -84.29692, 1.32590, 2021, 01, 26, 3, 16, 19},
0539         // mount axis x,y  pole center in x,y
0540         3855, 1117, 4022, 1623
0541     });
0542 
0543     // Brett's RASA data. Flipped parity images.
0544     const GeoLocation irvine(dms(-117, 50), dms(33, 33));
0545     runPAA(irvine,
0546     {
0547         { 231.85829, 89.49201, 92.16176, 2.51940, 2021, 03, 27, 4, 51, 14},
0548         { 195.35970, 89.53665, 85.77720, 2.51975, 2021, 03, 27, 4, 51, 30},
0549         { 162.20876, 89.54505, 77.99854, 2.51962, 2021, 03, 27, 4, 51, 44},
0550 
0551         // low error--don't know true coords
0552         // Polar Alignment Error:  00° 01' 46\". Azimuth:  00° 01' 18\"  Altitude: -00° 01' 12\""
0553 
0554         // Made up the starting point. Using the calculated target.
0555         500, 500, 509, 537
0556     },
0557 
0558     false);
0559 
0560     runPAA(irvine,
0561     {
0562         { 247.48287, 89.46797, 107.38986, 2.51988, 2021, 03, 27, 4, 52, 49},
0563         { 211.13546, 89.60233, 102.12462, 2.51959, 2021, 03, 27, 4, 53, 06},
0564         { 163.42722, 89.67487, 83.48466, 2.51931, 2021, 03, 27, 4, 53, 21},
0565 
0566         // Polar Alignment Error:  00° 09' 38\". Azimuth:  00° 01' 15\"  Altitude: -00° 09' 33\"
0567         // 2911, 824, 2725, 959  // without flip
0568         // 2911, 824, 3084, 916  // where Brett wound up
0569         2911, 824, 3097, 958     // what this algorithm calculates
0570     },
0571     // False means the parity is flipped.
0572     false);
0573 
0574     runPAA(irvine,
0575     {
0576         { 233.33273, 89.35997, 91.72542, 2.51947, 2021, 03, 27, 4, 59, 10},
0577         { 205.15803, 89.42125, 93.16138, 2.51971, 2021, 03, 27, 4, 59, 27},
0578         { 176.87893, 89.48116, 91.40070, 2.51964, 2021, 03, 27, 4, 59, 42},
0579 
0580         // Polar Alignment Error:  00° 10' 43\". Azimuth:  00° 10' 40\"  Altitude: -00° 00' 57\""
0581         // 2757, 318, 2854, 504, flip  // without flip
0582         // 2757, 318, 2657, 450        // where Brett wound up
0583         2757, 318, 2661, 506           // what this algorithm calculates.
0584     },
0585     // False means the parity is flipped.
0586     false);
0587 
0588     // From Jasem's friend's log. Error was larger than hardcoded search limits.
0589     // Fixed to adjust search limits.
0590     runPAA(kuwait,
0591     {
0592         {180.65170, -48.38820, 39.56825, 0.75773, 2022, 04, 07, 20, 02, 49},
0593         {198.36159, -49.33399, 39.61439, 0.75786, 2022, 04, 07, 20, 03, 13},
0594         {215.61283, -50.20872, 39.26446, 0.75801, 2022, 04, 07, 20, 03, 37},
0595         1924, 1452, -1, -1
0596     });
0597 }
0598 
0599 struct RefreshSolution
0600 {
0601     double azErr, altErr;
0602     double azAdj, altAdj;
0603 };
0604 
0605 void TestPolarAlign::testRefreshCoords()
0606 {
0607     const GeoLocation siliconValley(dms(-122, 10), dms(37, 26, 30));
0608 
0609     // Taken from the 5/29 log. Times are UTC, so Pacific + 7.
0610     Solution p1 = {211.174, 60.8994, 69.66035, 0.39157, 2022, 5, 30, 5, 11, 11 };
0611     Solution p2 = {233.324, 60.632,  69.66035, 0.39157, 2022, 5, 30, 5, 11, 34 };
0612     Solution p3 = {254.451, 60.3434, 69.66035, 0.39157, 2022, 5, 30, 5, 11, 57 };
0613 
0614     QVector<Solution> ps = {p1, p2, p3};
0615 
0616     // right at start Estimated current adjustment: Az 0.0' Alt 0.0' residual 4a-s"
0617     Solution r1 =  {254.454,  60.346, 0, 0.391548, 2022, 5, 30, 5, 13, 3 };
0618     // refresh 25, Estimated current adjustment: Az 0.0' Alt -28.0' residual 23a-s"
0619     Solution r2 =  {253.841,  60.054, 0, 0.391548, 2022, 5, 30, 5, 14, 31 };
0620     // refresh 26, Estimated current adjustment: Az 0.0' Alt -28.0' residual 26a-s"
0621     Solution r3 =  {253.842,  60.054, 0, 0.391548, 2022, 5, 30, 5, 14, 34 };
0622 
0623     // refresh 27, Estimated current adjustment: Az 11.0' Alt -23.0' residual 220a-s"
0624     Solution r4 =  {253.769,  60.207, 0, 0.391548, 2022, 5, 30, 5, 14, 48 };
0625     // refresh 28, Estimated current adjustment: Az 10.0' Alt -22.0' residual 265a-s"
0626     Solution r5 =  {253.769,  60.206, 0, 0.391548, 2022, 5, 30, 5, 14, 52 };
0627     // refresh 29, Estimated current adjustment: Az 17.0' Alt -19.0' residual 409a-s"
0628     Solution r6 =  {253.724,  60.297, 0, 0.391548, 2022, 5, 30, 5, 15, 2 };
0629     // refresh 36, Estimated current adjustment: Az 27.0' Alt -15.0' residual 607a-s"
0630     Solution r7 =  {253.656,  60.429, 0, 0.391548, 2022, 5, 30, 5, 15, 28 };
0631 
0632     QVector<Solution> rs = {r1, r2, r3, r4, r5, r6, r7};
0633 
0634 
0635     QVector<RefreshSolution> refreshSolutions =
0636     {
0637         // az error| alt error| az adjust | alt adjust
0638         { 0.629487,  -0.468870,  0.000000,  -0.001389 },
0639         { 0.638877,  -0.010564, -0.005556,  -0.459722 },
0640         { 0.638877,  -0.010564, -0.005556,  -0.459722 },
0641         { 0.387476,  -0.011953,  0.245833,  -0.458333 },
0642         { 0.387500,  -0.009175,  0.245833,  -0.461111 },
0643         { 0.234734,  -0.007787,  0.398611,  -0.462500 },
0644         { 0.016678,  -0.007787,  0.616667,  -0.462500 },
0645     };
0646 
0647 
0648     PolarAlign polarAlign(&siliconValley);
0649     foreach (const Solution &p, ps)
0650     {
0651         QSharedPointer<FITSData> image;
0652         setupData(p, image, true);
0653         QVERIFY(polarAlign.addPoint(image));
0654     }
0655     polarAlign.findAxis();
0656 
0657     double initialAzimuthError, initialAltitudeError;
0658     polarAlign.calculateAzAltError(&initialAzimuthError, &initialAltitudeError);
0659     dms polarError(hypot(initialAzimuthError, initialAltitudeError));
0660     dms azError(initialAzimuthError), altError(initialAltitudeError);
0661 
0662     int i = 0;
0663     foreach (const Solution &r, rs)
0664     {
0665         KStarsDateTime newTime;
0666         newTime.setDate(QDate(r.year, r.month, r.day));
0667         newTime.setTime(QTime(r.hour, r.minute, r.second));
0668         newTime.setTimeSpec(Qt::UTC);
0669 
0670         SkyPoint solution;
0671         SkyPoint refreshPoint(r.ra / 15.0, r.dec);
0672         double azErr, altErr, azAdjustment, altAdjustment;
0673         polarAlign.processRefreshCoords(refreshPoint, newTime, &azErr, &altErr, &azAdjustment, &altAdjustment);
0674 
0675         constexpr double smallError = .0001;
0676         QVERIFY(fabs(azErr - refreshSolutions[i].azErr) < smallError);
0677         QVERIFY(fabs(altErr - refreshSolutions[i].altErr) < smallError);
0678         QVERIFY(fabs(azAdjustment - refreshSolutions[i].azAdj) < smallError);
0679         QVERIFY(fabs(altAdjustment - refreshSolutions[i].altAdj) < smallError);
0680         i++;
0681     }
0682 }
0683 
0684 void TestPolarAlign::getAzAlt(const KStarsDateTime &time, const GeoLocation &geo,
0685                               const QPointF &pixel, double ra, double dec, double orientation,
0686                               double pixScale, double *az, double *alt)
0687 {
0688     QSharedPointer<FITSData> image;
0689     loadDummyFits(image, time, ra, dec, orientation, pixScale, true);
0690 
0691     SkyPoint pt;
0692     PolarAlign polarAlign(&geo);
0693     QVERIFY(polarAlign.prepareAzAlt(image, pixel, &pt));
0694     *az = pt.az().Degrees();
0695     *alt = pt.alt().Degrees();
0696 }
0697 
0698 void TestPolarAlign::testAlt()
0699 {
0700     // Silicon Valley
0701     const GeoLocation geo(dms(-122, 10), dms(37, 26, 30));
0702     constexpr double pixScale = 1.32577;
0703 
0704     KStarsDateTime time1;
0705     time1.setDate(QDate(2020, 12, 27));
0706     time1.setTime(QTime(3, 32, 18));
0707     time1.setTimeSpec(Qt::UTC);
0708 
0709 
0710     double az1 = 0, alt1 = 0;
0711     getAzAlt(time1, geo, QPointF(IMAGE_WIDTH / 2, IMAGE_HEIGHT / 2), 20.26276, 89.84129, 75.75402, pixScale, &az1, &alt1);
0712 
0713     KStarsDateTime time2;
0714     time2.setDate(QDate(2020, 12, 27));
0715     time2.setTime(QTime(3, 42, 10));
0716     time2.setTimeSpec(Qt::UTC);
0717 
0718     double az2 = 0, alt2 = 0;
0719     getAzAlt(time2, geo, QPointF(IMAGE_WIDTH / 2, IMAGE_HEIGHT / 2), 26.66804, 89.61473, 79.69920, pixScale, &az2, &alt2);
0720 
0721     // The first set of coordinates and times were taken were sampled, and then the mount's
0722     // altitude knob was adjusted so that the telescope was pointing higher.
0723     // The second set of coordinates should indicate a significant change in altitude
0724     // with only a small change in azimuth.
0725     compare(az1,  0.0452539, "Azimuth1", __LINE__);
0726     compare(az2,  0.0523052, "Azimuth2", __LINE__);
0727     compare(alt1, 37.491241, "Altitude1", __LINE__);
0728     compare(alt2, 37.720853, "Altitude2", __LINE__);
0729 }
0730 
0731 void TestPolarAlign::testRotate_data()
0732 {
0733     QTest::addColumn<double>("az");
0734     QTest::addColumn<double>("alt");
0735     QTest::addColumn<double>("deltaAz");
0736     QTest::addColumn<double>("deltaAlt");
0737     QTest::addColumn<double>("azRotated");
0738     QTest::addColumn<double>("altRotated");
0739 
0740     //                       Az      Alt     dAz    dAlt      Az result     Alt result
0741 
0742     QTest::newRow("") <<    0.0 <<   0.0  << 0.0  << 0.0 <<     0.0      <<  0.0;
0743     QTest::newRow("") <<    0.0 <<   0.0  << 0.3  << 0.4 <<     0.3      <<  0.4;
0744     QTest::newRow("") <<    1.0 <<   1.0  << 0.0  << 0.0 <<     1.0      <<  1.0;
0745     QTest::newRow("") <<    1.0 <<   1.0  << 0.3  << 0.4 <<     1.300146 <<  1.399939;
0746 
0747     // Point in the west, north of the equator, curves up and south with increasing alt.
0748     QTest::newRow("") <<  -30.0 <<  60.0  << 0.0  << 1.0 <<   -30.893174 << 60.862134;
0749     QTest::newRow("") <<  -30.0 <<  60.0  << 0.0  << 2.0 <<   -31.843561 << 61.716007;
0750     QTest::newRow("") <<  -30.0 <<  60.0  << 0.0  << 3.0 <<   -32.855841 << 62.560844;
0751     QTest::newRow("") <<  -30.0 <<  60.0  << 0.0  << 4.0 <<   -33.935125 << 63.395779;
0752 
0753     // Point in the east, north of the equator, curves up and south with increasing alt.
0754     QTest::newRow("") <<   30.0 <<  60.0  << 0.0  << 1.0 <<    30.893174 << 60.862134;
0755     QTest::newRow("") <<   30.0 <<  60.0  << 0.0  << 2.0 <<    31.843561 << 61.716007;
0756     QTest::newRow("") <<   30.0 <<  60.0  << 0.0  << 3.0 <<    32.855841 << 62.560844;
0757     QTest::newRow("") <<   30.0 <<  60.0  << 0.0  << 4.0 <<    33.935125 << 63.395779;
0758 
0759     // Point in the west, south of the equator, curves down and south with increasing alt.
0760     QTest::newRow("") << -110.0 <<  60.0  << 0.0  << 1.0 <<  -111.607689 << 59.644787;
0761     QTest::newRow("") << -110.0 <<  60.0  << 0.0  << 2.0 <<  -113.174586 << 59.263816;
0762     QTest::newRow("") << -110.0 <<  60.0  << 0.0  << 3.0 <<  -114.699478 << 58.858039;
0763     QTest::newRow("") << -110.0 <<  60.0  << 0.0  << 4.0 <<  -116.181474 << 58.428421;
0764 
0765     // Point in the east, south of the equator, curves down and south with increasing alt.
0766     QTest::newRow("") <<  110.0 <<  60.0  << 0.0  << 1.0 <<   111.607689 << 59.644787;
0767     QTest::newRow("") <<  110.0 <<  60.0  << 0.0  << 2.0 <<   113.174586 << 59.263816;
0768     QTest::newRow("") <<  110.0 <<  60.0  << 0.0  << 3.0 <<   114.699478 << 58.858039;
0769     QTest::newRow("") <<  110.0 <<  60.0  << 0.0  << 4.0 <<   116.181474 << 58.428421;
0770 
0771     // First points of last 4 series. This time dAlt is negative so goes the opposite way.
0772     QTest::newRow("") <<  -30.0 <<  60.0  << 0.0  << -1.0 <<  -29.159759 << 59.130303;
0773     QTest::newRow("") <<   30.0 <<  60.0  << 0.0  << -1.0 <<   29.159759 << 59.130303;
0774     QTest::newRow("") << -110.0 <<  60.0  << 0.0  << -1.0 << -108.353075 << 60.328521;
0775     QTest::newRow("") <<  110.0 <<  60.0  << 0.0  << -1.0 <<  108.353075 << 60.328521;
0776 
0777     // Same as last 4, but with dAz = 1
0778     QTest::newRow("") <<  -30.0 <<  60.0  << 1.0  << -1.0 <<  -28.159759 << 59.130303;
0779     QTest::newRow("") <<   30.0 <<  60.0  << 1.0  << -1.0 <<   30.159759 << 59.130303;
0780     QTest::newRow("") << -110.0 <<  60.0  << 1.0  << -1.0 << -107.353075 << 60.328521;
0781     QTest::newRow("") <<  110.0 <<  60.0  << 1.0  << -1.0 <<  109.353075 << 60.328521;
0782 
0783     // ditto, but with dAz = -1
0784     QTest::newRow("") <<  -30.0 <<  60.0  << -1.0 << -1.0 <<  -30.159759 << 59.130303;
0785     QTest::newRow("") <<   30.0 <<  60.0  << -1.0 << -1.0 <<   28.159759 << 59.130303;
0786     QTest::newRow("") << -110.0 <<  60.0  << -1.0 << -1.0 << -109.353075 << 60.328521;
0787     QTest::newRow("") <<  110.0 <<  60.0  << -1.0 << -1.0 <<  107.353075 << 60.328521;
0788 
0789     // Some much larger rotations
0790     QTest::newRow("") <<    1.0 <<  1.0  <<  20.0 << 20.0 <<   21.070967 << 20.996804;
0791     QTest::newRow("") <<    1.0 <<  1.0  <<   0.0 << 20.0 <<    1.070967 << 20.996804;
0792     QTest::newRow("") <<    1.0 <<  1.0  << -20.0 << 20.0 <<  -18.929033 << 20.996804;
0793     QTest::newRow("") <<  -30.0 << 60.0  <<  20.0 << 20.0 <<  -46.116102 << 74.132541;
0794     QTest::newRow("") <<   30.0 << 60.0  <<   0.0 << 20.0 <<   66.116102 << 74.132541;
0795     QTest::newRow("") << -110.0 << 60.0  << -20.0 << 20.0 << -154.199340 << 49.052359;
0796     QTest::newRow("") <<  110.0 << 60.0  <<  20.0 << -20.0 <<  93.912714 << 60.725451;
0797 }
0798 
0799 void TestPolarAlign::testRotate()
0800 {
0801     QFETCH(double, az);
0802     QFETCH(double, alt);
0803     QFETCH(double, deltaAz);
0804     QFETCH(double, deltaAlt);
0805     QFETCH(double, azRotated);
0806     QFETCH(double, altRotated);
0807 
0808     // PolarAlign is needed to specifiy the hemisphere, which specifies how the
0809     // altitude rotation works (positive alt error causes opposite rotations in
0810     // the northern vs southern hemisphere.
0811     GeoLocation geo(dms(-122, 10), dms(37, 26, 30));
0812     PolarAlign polarAlign(&geo);
0813 
0814     QPointF point, rot;
0815     point = QPointF(az, alt);
0816     rot = QPointF(deltaAz, deltaAlt);
0817     QPointF result = Rotations::rotateRaAxis(point, rot);
0818     QVERIFY(compare(result, azRotated, altRotated, .001));
0819 }
0820 
0821 QTEST_GUILESS_MAIN(TestPolarAlign)