fix(libastro): use apparent sidereal time in EquatorialToHorizontal; add accuracy tests

EquatorialToHorizontal called ln_get_hrz_from_equ which uses mean sidereal
time internally, while HorizontalToEquatorial's ln_get_equ_from_hrz uses
apparent sidereal time. The mismatch introduced a fixed ~17" RA offset
(equation of the equinoxes) in every round-trip. Fix: call
ln_get_hrz_from_equ_sidereal_time with ln_get_apparent_sidereal_time(JD)
explicitly, matching the inverse path.

Five tests document the accuracy of the libnova coordinate pipeline:

Reciprocity: J2000toObserved -> ObservedToJ2000 closes to 0.004" for Deneb
at two epochs (2020, 2026).

HorizontalAccuracy_vs_IMCCE: measures accuracy against IMCCE Miriade (INPOP19)
for Vega and Polaris at four sites (Greenwich, Mitaka, Mauna Kea, Siding Spring).
Observed errors:
  Vega:    az 5-12",  alt 3-6"
  Polaris: az 18-31", alt 3-11"
Polaris az errors are larger because polar distance is only 0.74°: a nutation-
in-RA error is amplified ~1/sin(0.74°) ≈ 77x when mapped to azimuth. The
underlying error is the IAU 1980 nutation accuracy floor of libnova, common to
both stars. Per-case errors are logged; the assertion catches gross regressions
which cause degree-level errors.

ObserverLongitudeMatters: Vega altitude spans > 30° across sites covering 316°
of longitude.

ObserverLatitudeMatters: Polaris is above the horizon at Greenwich (51.5 N) and
below it at Siding Spring (31.3 S).

RoundTrip_HorizontalToEquatorial: after the sidereal time fix, the round-trip
closes to sub-picosecond (was ~11" RA before the fix).

Golden data generated by tools/generate_golden_files.py (IMCCE Miriade, INPOP19).

Author: ckemper67

libs/indicore/libastro.cpp

@@ -39,6 +39,7 @@
 #include <libnova/aberration.h>
 #include <libnova/transform.h>
 #include <libnova/nutation.h>
+#include <libnova/sidereal_time.h>
 
 namespace INDI
 {
@@ -145,7 +146,11 @@ void EquatorialToHorizontal(IEquatorialCoordinates *object, IGeographicCoordinat
     // RA Hours --> Degrees
     struct ln_equ_posn libnova_object = {object->rightascension * 15.0, object->declination};
     struct ln_hrz_posn horizontalPos;
-    ln_get_hrz_from_equ(&libnova_object, &libnova_location, JD, &horizontalPos);
+    // Use apparent sidereal time to match ln_get_equ_from_hrz (which already uses apparent).
+    // ln_get_hrz_from_equ uses mean sidereal time — a libnova bug that introduces a ~17" RA
+    // offset (equation of the equinoxes) in the round-trip.
+    double apparentSidereal = ln_get_apparent_sidereal_time(JD);
+    ln_get_hrz_from_equ_sidereal_time(&libnova_object, &libnova_location, apparentSidereal, &horizontalPos);
     position->azimuth = range360(180 + horizontalPos.az);
     position->altitude = horizontalPos.alt;
 }
@@ -161,6 +166,7 @@ void HorizontalToEquatorial(IHorizontalCoordinates *object, IGeographicCoordinat
     // Convert from INDI standard location to libnova standard location
     struct ln_hrz_posn libnova_object = {range360(object->azimuth + 180), object->altitude};
     struct ln_equ_posn equatorialPos;
+    // ln_get_equ_from_hrz uses apparent sidereal time, matching EquatorialToHorizontal above.
     ln_get_equ_from_hrz(&libnova_object, &libnova_location, JD, &equatorialPos);
     // Degrees --> Hours
     position->rightascension = equatorialPos.ra / 15.0;

test/core/CMakeLists.txt

@@ -1,5 +1,3 @@
-
-
 SET (test_base64_SRCS
 	test_base64.cpp
 )
@@ -34,7 +32,7 @@ SET (test_rotator_limits_SRCS
 ADD_EXECUTABLE(test_rotator_limits
     ${test_rotator_limits_SRCS}
 )
-TARGET_LINK_LIBRARIES(test_rotator_limits   
+TARGET_LINK_LIBRARIES(test_rotator_limits
     ${GTEST_BOTH_LIBRARIES}
 	${GMOCK_LIBRARIES}
     ${CMAKE_THREAD_LIBS_INIT}
@@ -52,3 +50,18 @@ TARGET_LINK_LIBRARIES(test_libnova_nutation
     ${CMAKE_THREAD_LIBS_INIT}
 )
 ADD_TEST(test_libnova_nutation test_libnova_nutation)
+
+SET (test_libastro_SRCS
+    test_libastro.cpp
+)
+ADD_EXECUTABLE(test_libastro ${test_libastro_SRCS})
+TARGET_LINK_LIBRARIES(test_libastro
+    indidriver
+    ${GTEST_BOTH_LIBRARIES}
+    ${GMOCK_LIBRARIES}
+    ${CMAKE_THREAD_LIBS_INIT}
+)
+TARGET_COMPILE_DEFINITIONS(test_libastro PRIVATE
+    TEST_DATA_DIR="${CMAKE_SOURCE_DIR}/test/data"
+)
+ADD_TEST(test_libastro test_libastro)

test/core/test_libastro.cpp

@@ -0,0 +1,253 @@
+/*******************************************************************************
+ * Tests for libastro coordinate functions (libnova backend)
+ *
+ * Tests for the libastro coordinate pipeline against external truth (IMCCE Miriade)
+ * and internal round-trip consistency.
+ *
+ * Golden data: test/data/horizontal_golden.json
+ *   8 cases: Vega and Polaris at Greenwich, Mitaka, Mauna Kea, Siding Spring.
+ *   Source: IMCCE Miriade (INPOP19), -tcoor=5 (horizontal coordinates).
+ *   Regenerate with: python3 tools/generate_golden_files.py
+ *******************************************************************************/
+
+#include <gtest/gtest.h>
+#include <libastro.h>
+#include <nlohmann/json.hpp>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <cmath>
+
+// ---------------------------------------------------------------------------
+// Golden data loader
+// ---------------------------------------------------------------------------
+
+struct HorizCase {
+    std::string label;
+    std::string object;
+    double ra_j2000_h;
+    double dec_j2000_deg;
+    double jd;
+    std::string site;
+    double lon_deg;
+    double lat_deg;
+    double elev_m;
+    double az_truth;
+    double alt_truth;
+};
+
+static std::vector<HorizCase> load_golden()
+{
+    std::ifstream f(TEST_DATA_DIR "/horizontal_golden.json");
+    if (!f.is_open()) return {};
+    nlohmann::json j = nlohmann::json::parse(f);
+    std::vector<HorizCase> cases;
+    for (auto &e : j) {
+        if (!e.contains("az_deg")) continue;
+        HorizCase c;
+        c.label         = e["label"];
+        c.object        = e["object"];
+        c.ra_j2000_h    = e["ra_j2000_h"];
+        c.dec_j2000_deg = e["dec_j2000_deg"];
+        c.jd            = e["jd"];
+        c.site          = e["site"];
+        c.lon_deg       = e["lon_deg"];
+        c.lat_deg       = e["lat_deg"];
+        c.elev_m        = e["elev_m"];
+        c.az_truth      = e["az_deg"];
+        c.alt_truth     = e["alt_deg"];
+        cases.push_back(c);
+    }
+    return cases;
+}
+
+static double az_diff_arcsec(double a, double b) {
+    double d = a - b;
+    while (d >  180.0) d -= 360.0;
+    while (d < -180.0) d += 360.0;
+    return std::abs(d) * 3600.0;
+}
+
+// ---------------------------------------------------------------------------
+// J2000toObserved -> ObservedToJ2000 round-trip
+// The libnova forward and inverse paths must close to < 1 arcsecond.
+// A large error here indicates a broken nutation or precession inversion.
+// ---------------------------------------------------------------------------
+TEST(Libastro, Reciprocity)
+{
+    // Use two different epochs to exercise different precession amounts.
+    struct { double jd; const char *label; } epochs[] = {
+        { 2459019.833333, "2020-06-19" },
+        { 2461112.5,      "2026-01-01" },
+    };
+    INDI::IEquatorialCoordinates j2000_in = { 20.69053168, 45.28033881 }; // Deneb
+
+    for (auto &ep : epochs) {
+        INDI::IEquatorialCoordinates jnow, j2000_out;
+        INDI::J2000toObserved(&j2000_in, ep.jd, &jnow);
+        INDI::ObservedToJ2000(&jnow, ep.jd, &j2000_out);
+
+        double cos_dec = std::cos(j2000_in.declination * M_PI / 180.0);
+        double ra_err  = std::abs(j2000_in.rightascension - j2000_out.rightascension) * 15.0 * 3600.0 * cos_dec;
+        double dec_err = std::abs(j2000_in.declination    - j2000_out.declination)    * 3600.0;
+
+        GTEST_LOG_(INFO) << ep.label << " round-trip: RA=" << ra_err << "\"  Dec=" << dec_err << "\"";
+        EXPECT_LT(ra_err,  1.0) << "Round-trip RA error at " << ep.label;
+        EXPECT_LT(dec_err, 1.0) << "Round-trip Dec error at " << ep.label;
+    }
+}
+
+// ---------------------------------------------------------------------------
+// EquatorialToHorizontal accuracy vs IMCCE Miriade (INPOP19) truth.
+// Per-case errors are logged; the assertion catches gross regressions
+// (wrong observer, wrong sidereal time, etc.) which cause degree-level errors.
+// ---------------------------------------------------------------------------
+TEST(Libastro, HorizontalAccuracy_vs_IMCCE)
+{
+    static constexpr double TOLERANCE_ARCSEC = 60.0;
+
+    auto cases = load_golden();
+    ASSERT_GT(cases.size(), 0u) << "Could not load horizontal_golden.json";
+
+    double max_az_err = 0, max_alt_err = 0;
+    int n = 0;
+    for (auto &c : cases) {
+        INDI::IEquatorialCoordinates j2000 = { c.ra_j2000_h, c.dec_j2000_deg };
+        INDI::IEquatorialCoordinates jnow;
+        INDI::J2000toObserved(&j2000, c.jd, &jnow);
+
+        INDI::IGeographicCoordinates obs = { c.lon_deg, c.lat_deg, c.elev_m };
+        INDI::IHorizontalCoordinates hrz;
+        INDI::EquatorialToHorizontal(&jnow, &obs, c.jd, &hrz);
+
+        double az_err  = az_diff_arcsec(hrz.azimuth,  c.az_truth);
+        double alt_err = std::abs(hrz.altitude - c.alt_truth) * 3600.0;
+        max_az_err  = std::max(max_az_err,  az_err);
+        max_alt_err = std::max(max_alt_err, alt_err);
+
+        GTEST_LOG_(INFO) << c.object << " at " << c.site
+                         << "  az_err=" << az_err << "\"  alt_err=" << alt_err << "\"";
+        EXPECT_LT(az_err,  TOLERANCE_ARCSEC)
+            << "Az error for " << c.object << " at " << c.site
+            << ": " << az_err << "\" (truth=" << c.az_truth << " got=" << hrz.azimuth << ")";
+        EXPECT_LT(alt_err, TOLERANCE_ARCSEC)
+            << "Alt error for " << c.object << " at " << c.site
+            << ": " << alt_err << "\" (truth=" << c.alt_truth << " got=" << hrz.altitude << ")";
+        n++;
+    }
+    GTEST_LOG_(INFO) << "Horizontal accuracy (" << n << " cases):"
+                     << " max az_err=" << max_az_err << "\""
+                     << " max alt_err=" << max_alt_err << "\"";
+}
+
+// ---------------------------------------------------------------------------
+// Structural: observer longitude must change altitude significantly.
+// If EquatorialToHorizontal ignores the observer position, Vega has the same
+// altitude at all sites. Sites span ~316 deg of longitude so the spread must
+// exceed 30 deg.
+// ---------------------------------------------------------------------------
+TEST(Libastro, ObserverLongitudeMatters)
+{
+    auto cases = load_golden();
+    ASSERT_GT(cases.size(), 0u);
+
+    double alt_min = 999.0, alt_max = -999.0;
+    int vega_count = 0;
+    for (auto &c : cases) {
+        if (c.object != "Vega") continue;
+        INDI::IEquatorialCoordinates j2000 = { c.ra_j2000_h, c.dec_j2000_deg };
+        INDI::IEquatorialCoordinates jnow;
+        INDI::J2000toObserved(&j2000, c.jd, &jnow);
+        INDI::IGeographicCoordinates obs = { c.lon_deg, c.lat_deg, c.elev_m };
+        INDI::IHorizontalCoordinates hrz;
+        INDI::EquatorialToHorizontal(&jnow, &obs, c.jd, &hrz);
+        alt_min = std::min(alt_min, hrz.altitude);
+        alt_max = std::max(alt_max, hrz.altitude);
+        vega_count++;
+    }
+    ASSERT_GE(vega_count, 3) << "Need at least 3 Vega cases";
+    EXPECT_GT(alt_max - alt_min, 30.0)
+        << "Vega altitude spread=" << (alt_max - alt_min)
+        << " deg — expected >30 deg. Observer longitude may be ignored.";
+    GTEST_LOG_(INFO) << "Vega altitude spread: " << alt_min << " to " << alt_max
+                     << " deg (spread=" << (alt_max - alt_min) << " deg)";
+}
+
+// ---------------------------------------------------------------------------
+// Structural: observer latitude must determine whether Polaris is above the
+// horizon. Polaris is circumpolar from Greenwich (lat 51.5 N) and below the
+// horizon from Siding Spring (lat 31.3 S). If latitude is ignored, the sign
+// of the altitude is wrong at one or both sites.
+// ---------------------------------------------------------------------------
+TEST(Libastro, ObserverLatitudeMatters)
+{
+    auto cases = load_golden();
+    ASSERT_GT(cases.size(), 0u);
+
+    double alt_greenwich = 999.0, alt_siding = 999.0;
+    for (auto &c : cases) {
+        if (c.object != "Polaris") continue;
+        INDI::IEquatorialCoordinates j2000 = { c.ra_j2000_h, c.dec_j2000_deg };
+        INDI::IEquatorialCoordinates jnow;
+        INDI::J2000toObserved(&j2000, c.jd, &jnow);
+        INDI::IGeographicCoordinates obs = { c.lon_deg, c.lat_deg, c.elev_m };
+        INDI::IHorizontalCoordinates hrz;
+        INDI::EquatorialToHorizontal(&jnow, &obs, c.jd, &hrz);
+        if (c.site == "Greenwich")     alt_greenwich = hrz.altitude;
+        if (c.site == "Siding Spring") alt_siding    = hrz.altitude;
+    }
+    ASSERT_NE(alt_greenwich, 999.0) << "Greenwich Polaris case missing";
+    ASSERT_NE(alt_siding,    999.0) << "Siding Spring Polaris case missing";
+
+    EXPECT_GT(alt_greenwich, 0.0)
+        << "Polaris should be above horizon at Greenwich (got " << alt_greenwich << " deg)";
+    EXPECT_LT(alt_siding, 0.0)
+        << "Polaris should be below horizon at Siding Spring (got " << alt_siding << " deg)";
+    GTEST_LOG_(INFO) << "Polaris: Greenwich=" << alt_greenwich
+                     << " deg, Siding Spring=" << alt_siding << " deg";
+}
+
+// ---------------------------------------------------------------------------
+// Round-trip: EquatorialToHorizontal -> HorizontalToEquatorial must recover
+// the original JNow coordinates to floating-point precision (sub-arcsecond).
+// Skips cases below 5 deg altitude where spherical trig is less stable.
+// ---------------------------------------------------------------------------
+TEST(Libastro, RoundTrip_HorizontalToEquatorial)
+{
+    auto cases = load_golden();
+    ASSERT_GT(cases.size(), 0u);
+
+    double max_ra_err = 0, max_dec_err = 0;
+    for (auto &c : cases) {
+        INDI::IEquatorialCoordinates j2000 = { c.ra_j2000_h, c.dec_j2000_deg };
+        INDI::IEquatorialCoordinates jnow_in;
+        INDI::J2000toObserved(&j2000, c.jd, &jnow_in);
+
+        INDI::IGeographicCoordinates obs = { c.lon_deg, c.lat_deg, c.elev_m };
+        INDI::IHorizontalCoordinates hrz;
+        INDI::EquatorialToHorizontal(&jnow_in, &obs, c.jd, &hrz);
+
+        if (hrz.altitude < 5.0) continue;
+
+        INDI::IEquatorialCoordinates jnow_out;
+        INDI::HorizontalToEquatorial(&hrz, &obs, c.jd, &jnow_out);
+
+        double cos_dec = std::cos(jnow_in.declination * M_PI / 180.0);
+        double ra_err  = std::abs(jnow_in.rightascension - jnow_out.rightascension) * 15.0 * 3600.0 * cos_dec;
+        double dec_err = std::abs(jnow_in.declination    - jnow_out.declination)    * 3600.0;
+        max_ra_err  = std::max(max_ra_err,  ra_err);
+        max_dec_err = std::max(max_dec_err, dec_err);
+
+        GTEST_LOG_(INFO) << c.object << " at " << c.site
+                         << "  ra_err=" << ra_err << "\"  dec_err=" << dec_err << "\"";
+        EXPECT_LT(ra_err,  1.0) << "Round-trip RA error for "  << c.object << " at " << c.site;
+        EXPECT_LT(dec_err, 1.0) << "Round-trip Dec error for " << c.object << " at " << c.site;
+    }
+    GTEST_LOG_(INFO) << "Round-trip max: RA=" << max_ra_err << "\" Dec=" << max_dec_err << "\"";
+}
+
+int main(int argc, char **argv)
+{
+    ::testing::InitGoogleTest(&argc, argv);
+    return RUN_ALL_TESTS();
+}