#!/usr/bin/env python3
"""Regressionstest fuer den Bezugssystem-Fehler in skyfield_planet_state().

apparent.radec() muss ohne epoch-Argument aufgerufen werden, damit die
Ausgabe echte ICRF/J2000-Koordinaten liefert. Ein wieder eingefuehrtes
epoch=timescale.J2000 wuerde die Nutation von J2000.0 (~14,5" bei diesem
Merkur-Test) mit einrechnen und diesen Test zum Scheitern bringen.

Aufruf: python test_ephemeriden_api.py
"""
import math
import os
import sys
import unittest
from datetime import datetime, timezone

SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
if SCRIPT_DIR not in sys.path:
    sys.path.insert(0, SCRIPT_DIR)

import ephemeriden_api as api
from astronomical_conversions import skyfield_wgs84

# Sternwarte Sonneberg
TEST_LATITUDE = 50.3773963
TEST_LONGITUDE = 11.1898541
TEST_ELEVATION_M = 640.0
TEST_INSTANT_UTC = datetime(2026, 6, 15, 12, 0, 0, tzinfo=timezone.utc)

# Erwartete topozentrische, scheinbare ICRF-Koordinaten (radec() ohne epoch=),
# berechnet mit DE421 fuer obigen Ort/Zeitpunkt. Geocentric-astrometrisches
# Gegenstueck stimmt mit einer unabhaengigen VSOP87-Rechnung (PyEphem) auf
# rund 0,3" ueberein.
EXPECTED_RA_HOURS = 7.348643358806561
EXPECTED_DEC_DEG = 23.20555540423118
TOLERANCE_ARCSEC = 1.0


def angular_separation_arcsec(ra1_hours: float, dec1_deg: float, ra2_hours: float, dec2_deg: float) -> float:
    ra1 = math.radians(ra1_hours * 15.0)
    ra2 = math.radians(ra2_hours * 15.0)
    dec1 = math.radians(dec1_deg)
    dec2 = math.radians(dec2_deg)
    cos_sep = (
        math.sin(dec1) * math.sin(dec2)
        + math.cos(dec1) * math.cos(dec2) * math.cos(ra1 - ra2)
    )
    cos_sep = max(-1.0, min(1.0, cos_sep))
    return math.degrees(math.acos(cos_sep)) * 3600.0


class SkyfieldReferenceFrameTest(unittest.TestCase):
    def test_mercury_position_matches_icrf_reference(self):
        timescale, planets = api.skyfield_context()
        observer = planets["earth"] + skyfield_wgs84.latlon(
            TEST_LATITUDE, TEST_LONGITUDE, elevation_m=TEST_ELEVATION_M
        )

        ra_hours, dec_deg, _azimuth_deg, _altitude_deg = api.skyfield_planet_state(
            "Mercury", observer, TEST_INSTANT_UTC, timescale, planets
        )

        separation_arcsec = angular_separation_arcsec(
            ra_hours, dec_deg, EXPECTED_RA_HOURS, EXPECTED_DEC_DEG
        )
        self.assertLess(
            separation_arcsec,
            TOLERANCE_ARCSEC,
            f"Position weicht um {separation_arcsec:.3f}\" vom Sollwert ab "
            f"(Toleranz {TOLERANCE_ARCSEC}\"). Wurde epoch= bei radec() wieder "
            "eingefuehrt?",
        )


if __name__ == "__main__":
    unittest.main()