diff --git a/public/py/ephemeriden_api.py b/public/py/ephemeriden_api.py index a09acc6..85aa982 100644 --- a/public/py/ephemeriden_api.py +++ b/public/py/ephemeriden_api.py @@ -145,7 +145,7 @@ def skyfield_planet_state(body_name: str, observer, dt_utc: datetime, timescale= timescale, planets = skyfield_context() time_value = timescale.from_datetime(dt_utc.astimezone(timezone.utc)) apparent = observer.at(time_value).observe(planets[SKYFIELD_EPHEMERIS_TARGETS[body_name]]).apparent() - ra, dec, _distance = apparent.radec(epoch=timescale.J2000) + ra, dec, _distance = apparent.radec() altitude, azimuth, _distance = apparent.altaz("standard") return float(ra.hours), float(dec.degrees), float(azimuth.degrees), float(altitude.degrees) diff --git a/public/py/test_ephemeriden_api.py b/public/py/test_ephemeriden_api.py new file mode 100644 index 0000000..20b017c --- /dev/null +++ b/public/py/test_ephemeriden_api.py @@ -0,0 +1,76 @@ +#!/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()