Add new polaralign class including test harness

Tests/CMakeLists.txt

@@ -4,6 +4,7 @@ include_directories(
     ${kstars_SOURCE_DIR}/kstars/skycomponents
     ${kstars_SOURCE_DIR}/kstars/auxiliary
     ${kstars_SOURCE_DIR}/kstars/time
+    ${kstars_SOURCE_DIR}/kstars/ekos/align
     )
 
 if(WCSLIB_FOUND)
@@ -32,6 +33,7 @@ SET( TEST_KSLITE_LIBRARIES
 
 add_subdirectory(auxiliary)
 add_subdirectory(skyobjects)
+add_subdirectory(polaralign)
 
 IF (CFITSIO_FOUND)
     add_subdirectory(fitsviewer)

Tests/polaralign/CMakeLists.txt

@@ -0,0 +1,3 @@
+ADD_EXECUTABLE( test_polaralign test_polaralign.cpp )
+TARGET_LINK_LIBRARIES( test_polaralign ${TEST_LIBRARIES})
+ADD_TEST( NAME TestPolarAlign COMMAND test_polaralign )

Tests/polaralign/test_polaralign.cpp

@@ -0,0 +1,161 @@
+/***************************************************************************
+                   test_polaralign.cpp  -  KStars Planetarium
+                             -------------------
+    begin                : Tue 27 Sep 2016 20:54:28 CDT
+    copyright            : (c) 2016 by Akarsh Simha
+    email                : akarsh.simha@kdemail.net
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+/* Project Includes */
+#include "test_polaralign.h"
+#include "ksnumbers.h"
+#include "time/kstarsdatetime.h"
+#include "auxiliary/dms.h"
+#include "Options.h"
+#include <libnova/libnova.h>
+
+TestPolarAlign::TestPolarAlign() : QObject()
+{
+}
+
+TestPolarAlign::~TestPolarAlign()
+{
+}
+
+void TestPolarAlign::compare(QVector3D v, double x, double y, double z)
+{
+    QVERIFY2(std::fabs(v.x() - x) < 0.00001,
+             qPrintable(QString("dc.x %1, x %2 error %3").arg(v.x()).arg(x).arg(((v.x() - x) * 3600.0), 6,'f', 1)));
+    QVERIFY2(std::fabs(v.y() - y) < 0.00001,
+             qPrintable(QString("dc.y %1, y %2 error %3").arg(v.y()).arg(y).arg(((v.y() - y) * 3600.0), 6,'f', 1)));
+    QVERIFY2(std::fabs(v.z() - z) < 0.00001,
+             qPrintable(QString("dc.z %1, z %2 error %3").arg(v.z()).arg(z).arg(((v.z() - z) * 3600.0), 6,'f', 1)));
+}
+
+void TestPolarAlign::compare(double a, double e, QString msg)
+{
+    QVERIFY2(std::fabs(a - e) < 0.0003,
+             qPrintable(QString("%1: actual %2, expected %3 error %4").arg(msg).arg(a).arg(e).arg(((a - e) * 3600.0), 6,'f', 1)));
+}
+
+
+void TestPolarAlign::testDirCos_data()
+{
+    QTest::addColumn<double>("Ha");
+    QTest::addColumn<double>("Dec");
+    QTest::addColumn<double>("X");
+    QTest::addColumn<double>("Y");
+    QTest::addColumn<double>("Z");
+
+    QTest::newRow("HaDec0") << 0.0 << 0.0  << 1.0 << 0.0 << 0.0;
+    QTest::newRow("Ha0Dec45") << 0.0 << 45.0  << 0.707107 << 0.0 << 0.707107;
+    QTest::newRow("Ha6Dec45") << 6.0 << 45.0 << 0.0 << 0.707107 << 0.707107;
+    QTest::newRow("Ha-6Dec45") << -6.0 << 45.0 << 0.0 << -0.707107 << 0.707107;
+    QTest::newRow("at Pole") << 0.0 << 90.0 << 0.0 << 0.0 << 1.0;
+    QTest::newRow("near S Pole") << -3.0 << -85.0 << 0.0616284 << -0.0616284 << -0.996195;
+}
+
+void TestPolarAlign::testDirCos()
+{
+    dms h;
+    dms d;
+    QFETCH(double, Ha);
+    QFETCH(double, Dec);
+    h.setH(Ha);
+    d.setD(Dec);
+
+
+    QVector3D dc;
+    dc = PolarAlign::dirCos(h, d);
+
+
+    QFETCH(double, X);
+    QFETCH(double, Y);
+    QFETCH(double, Z);
+
+    compare(dc, X, Y, Z);
+
+    SkyPoint sp(Ha, Dec);
+    dc = PolarAlign::dirCos(sp);
+    QCOMPARE(dc.length(), 1.0);
+    compare(dc, X, Y, Z);
+}
+
+void TestPolarAlign::testPriSec_data()
+{
+    testDirCos_data();
+}
+
+void TestPolarAlign::testPriSec()
+{
+    QFETCH(double, Ha);
+    QFETCH(double, Dec);
+    QFETCH(double, X);
+    QFETCH(double, Y);
+    QFETCH(double, Z);
+    QVector3D dc(X, Y, Z);
+    dms p = PolarAlign::primary(dc);
+    compare(p.HoursHa(), Ha, "Ha");
+    compare(PolarAlign::secondary(dc).Degrees(), Dec, "Dec");
+}
+
+
+void TestPolarAlign::testPoleAxis_data()
+{
+    QTest::addColumn<double>("Ha1");
+    QTest::addColumn<double>("Dec1");
+    QTest::addColumn<double>("Ha2");
+    QTest::addColumn<double>("Dec2");
+    QTest::addColumn<double>("Ha3");
+    QTest::addColumn<double>("Dec3");
+    QTest::addColumn<double>("X");
+    QTest::addColumn<double>("Y");
+    QTest::addColumn<double>("Z");
+
+    QTest::newRow("Ha-606Dec0") << -6.0 << 0.0 << 0.0 << 0.0  << 6.0 << 0.0  << 0.0 << 0.0 << 1.0;
+    QTest::newRow("Ha20-2Dec89") << 2.0 << 89.0 << 0.0 << 89.0  << -2.0 << 89.0  << 0.0 << 0.0 << -1.0;
+    QTest::newRow("Ha0-22Dec89v") << 0.0 << 89.0 << -2.0 << 89.1  << 2.0 << 88.9  << -0.0006 << -0.003386 << -0.99999;
+    QTest::newRow("Ha2-20Dec-89") << 2.0 << -89.0  << -2.0 << -89.0 << 0.0 << -89.0  << 0.0 << 0.0 << 1.0;
+    QTest::newRow("Ha20-2Dec89") << 2.0 << 89.0  << 0.0 << 89.0 << -2.0 << 88.0  << 0.05633683 << 0.0150954 << 0.998298;
+    // failure cases, 2 or more points the same should ruturn a null matrix
+    QTest::newRow("Ha000Dec0") << 0.0 << 0.0 << 0.0 << 0.0  << 0.0 << 0.0  << 0.0 << 0.0 << 0.0;
+    QTest::newRow("Ha100Dec0") << 1.0 << 0.0 << 0.0 << 0.0  << 0.0 << 0.0  << 0.0 << 0.0 << 0.0;
+    QTest::newRow("Ha110Dec0") << 1.0 << 0.0 << 1.0 << 0.0  << 0.0 << 0.0  << 0.0 << 0.0 << 0.0;
+    QTest::newRow("Ha011Dec0") << 0.0 << 0.0 << 1.0 << 0.0  << 1.0 << 0.0  << 0.0 << 0.0 << 0.0;
+}
+
+void TestPolarAlign::testPoleAxis()
+{
+    QFETCH(double, Ha1);
+    QFETCH(double, Dec1);
+    QFETCH(double, Ha2);
+    QFETCH(double, Dec2);
+    QFETCH(double, Ha3);
+    QFETCH(double, Dec3);
+
+    QFETCH(double, X);
+    QFETCH(double, Y);
+    QFETCH(double, Z);
+
+    SkyPoint p1(Ha1, Dec1);
+    SkyPoint p2(Ha2, Dec2);
+    SkyPoint p3(Ha3, Dec3);
+
+
+    QVector3D pa = PolarAlign::poleAxis(p1, p2, p3);
+
+    compare(pa.x(), X, "X");
+    compare(pa.y(), Y, "Y");
+    compare(pa.z(), Z, "Z");
+}
+
+QTEST_GUILESS_MAIN(TestPolarAlign)

Tests/polaralign/test_polaralign.h

@@ -0,0 +1,58 @@
+/***************************************************************************
+                     test_polaralign.h  -  KStars Planetarium
+                             -------------------
+    begin                : Tue 27 Sep 2016 20:51:21 CDT
+    copyright            : (c) 2020 by Chris Rowland
+    email                : 
+***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TEST_POLARALIGN_H
+#define TEST_POLARALIGN_H
+
+#include <QtTest/QtTest>
+#include <QDebug>
+#include <QString>
+
+#define UNIT_TEST
+
+#include "../../kstars/ekos/align/polaralign.h"
+
+/**
+ * @class TestPolarAlign
+ * @short Tests for some polar align operations
+ * @author Chris Rowland
+ */
+
+class TestPolarAlign : public QObject
+{
+        Q_OBJECT
+
+    public:
+        TestPolarAlign();
+        ~TestPolarAlign() override;
+
+private slots:
+        void testDirCos_data();
+        void testDirCos();
+        void testPriSec_data();
+        void testPriSec();
+        void testPoleAxis_data();
+        void testPoleAxis();
+
+private:
+        void compare(QVector3D v, double x, double y, double z);
+        void compare(double a, double e, QString msg = "");
+        void compare(float a, double e, QString msg = "") { compare (static_cast<double>(a), e, msg); }
+};
+
+
+#endif

kstars/CMakeLists.txt

@@ -193,6 +193,7 @@ if (INDI_FOUND)
             ekos/align/onlineastrometryparser.cpp
             ekos/align/remoteastrometryparser.cpp
             ekos/align/astapastrometryparser.cpp
+            ekos/align/polaralign.cpp
 
             # Guide
             ekos/guide/guide.cpp

kstars/auxiliary/dms.h

@@ -172,12 +172,18 @@ class dms
          */
     int msecond() const;
 
-    /** @return angle in hours expressed as a double.
+    /** @return angle in hours expressed as a double in the range 0 to 23.999...
          * @note an angle can be measured in degrees/arcminutes/arcseconds
          * or hours/minutes/seconds.  An hour is equal to 15 degrees.
          */
     inline double Hours() const { return reduce().Degrees() / 15.0; }
 
+    /** @return angle in hours expressed as a double in the range -11.999 to 0 to 12.0
+         * @note an angle can be measured in degrees/arcminutes/arcseconds
+         * or hours/minutes/seconds.  An hour is equal to 15 degrees.
+         */
+    inline double HoursHa() const { return Hours() <= 12.0 ? Hours() : Hours() - 24.0; }
+
     /** Sets floating-point value of angle, in degrees.
          * @param x new angle (double)
          */

kstars/ekos/align/polaralign.cpp

@@ -0,0 +1,65 @@
+/*  polaralign.cpp  determines the mount polar axis position
+
+    Copyright (C) 2020 Chris Rowland <chris.rowland@cherryfield.me.uk>
+
+    This application is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public
+    License as published by the Free Software Foundation; either
+    version 2 of the License, or (at your option) any later version.
+ */
+
+#include "polaralign.h"
+
+#include <cmath>
+
+QVector3D PolarAlign::dirCos(const dms primary, const dms secondary)
+{
+    return QVector3D(
+                static_cast<float>(secondary.cos() * primary.cos()),
+                static_cast<float>(secondary.cos() * primary.sin()),
+                static_cast<float>(secondary.sin()));
+}
+
+QVector3D PolarAlign::dirCos(const SkyPoint sp)
+{
+    return dirCos(sp.ra(), sp.dec());
+}
+
+dms PolarAlign::primary(QVector3D dirCos)
+{
+    dms p;
+    p.setRadians(static_cast<double>(std::atan2(dirCos.y(), dirCos.x())));
+    return p;
+}
+
+dms PolarAlign::secondary(QVector3D dirCos)
+{
+    dms p;
+    p.setRadians(static_cast<double>(std::asin(dirCos.z())));
+    return p;
+}
+
+SkyPoint PolarAlign::skyPoint(QVector3D dc)
+{
+    return SkyPoint(primary(dc), secondary(dc));
+}
+
+QVector3D PolarAlign::poleAxis(SkyPoint p1, SkyPoint p2, SkyPoint p3)
+{
+    // convert the three positions to vectors, these define the plane
+    // of the Ha axis rotation
+    QVector3D v1 = PolarAlign::dirCos(p1);
+    QVector3D v2 = PolarAlign::dirCos(p2);
+    QVector3D v3 = PolarAlign::dirCos(p3);
+
+    // the Ha axis direction is the normal to the plane
+    QVector3D p = QVector3D::normal(v1, v2, v3);
+
+    // p points to the north or south pole depending on the rotation of the points
+    // the other pole position can be determined by reversing the sign of the Dec and
+    // adding 12hrs to the Ha value.
+    // if we want only the north then this would do it
+    //if (p.z() < 0)
+    //    p = -p;
+    return p;
+}

kstars/ekos/align/polaralign.h

@@ -0,0 +1,76 @@
+/*  polaralign.h  determines the mount polar axis position
+
+    Copyright (C) 2020 Chris Rowland <chris.rowland@cherryfield.me.uk>
+
+    This application is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public
+    License as published by the Free Software Foundation; either
+    version 2 of the License, or (at your option) any later version.
+ */
+
+#ifndef POLARALIGN_H
+#define POLARALIGN_H
+
+#include <QVector3D>
+#include <dms.h>
+#include <skypoint.h>
+
+/**
+ *@class PolarAlign
+ *@short PolarAlign class handles determining the mount Ha axis position given three positions taken with the same mount declination.
+ *
+ *@author Chris Rowland
+ *@version 1.0
+ */
+class PolarAlign
+{
+public:
+    ///
+    /// \brief dirCos converts primary and secondary angles to a directional cosine
+    /// \param primary angle, can be Ra, Ha, Azimuth or the corresponding axis values
+    /// \param secondary angle, can be Dec, Altitude.  90 deg is the pole
+    /// \return QVector3D containing the directional cosine.
+    static QVector3D dirCos(const dms primary, const dms secondary);
+
+    ///
+    /// \brief dirCos converts a SkyPoint to a directional cosine
+    /// \param sp SkyPoint with the position
+    /// \return  QVector3D containing the directional cosine.
+    ///
+    static QVector3D dirCos(const SkyPoint sp);
+
+    ///
+    /// \brief primary returns the primary dms value in the directional cosine
+    /// \param dirCos
+    /// \return primary angle, Ra, Ha, Azimuth etc.
+    ///
+    static dms primary(QVector3D dirCos);
+
+    ///
+    /// \brief secondary returns the secondary dms angle in the directional cosine
+    /// \param dirCos
+    /// \return
+    ///
+    static dms secondary(QVector3D dirCos);
+
+    ///
+    /// \brief skyPoint returns a skypoint derived from the directional cosine vector
+    /// \param dc
+    /// \return
+    ///
+    static SkyPoint skyPoint(QVector3D dc);
+
+    ///
+    /// \brief poleAxis returns the pole axis vector given three SkyPoints with the same mount declination
+    /// \param p1
+    /// \param p2
+    /// \param p3
+    /// \return vector giving the direction of the pole. The rotation between the three points determines which pole
+    /// the other pole can be determined either by reversing the sign of the declination and adding 12 hrs to the Ha or
+    /// by negating the vector
+    ///
+    static QVector3D poleAxis(SkyPoint p1, SkyPoint p2, SkyPoint p3);
+
+};
+
+#endif // POLARALIGN_H