Bezugssystem-Fehler in Skyfield-Ephemeridenrechnung beheben
apparent.radec(epoch=timescale.J2000) lieferte faelschlich das wahre Aequinoktium/Aequator von J2000.0 statt ICRF-Koordinaten und rechnete damit die Nutation von J2000.0 (~14,5") in die Merkur/Planeten-Ausgabe ein, obwohl die CSV-Spalten "RA (J2000)"/"Dek (J2000)" ICRF-Koordinaten meinen sollen. radec() ohne epoch-Argument liefert das korrekte Bezugssystem. Regressionstest ergaenzt, der bei erneutem epoch=-Argument fehlschlaegt. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
This commit is contained in:
@@ -145,7 +145,7 @@ def skyfield_planet_state(body_name: str, observer, dt_utc: datetime, timescale=
|
|||||||
timescale, planets = skyfield_context()
|
timescale, planets = skyfield_context()
|
||||||
time_value = timescale.from_datetime(dt_utc.astimezone(timezone.utc))
|
time_value = timescale.from_datetime(dt_utc.astimezone(timezone.utc))
|
||||||
apparent = observer.at(time_value).observe(planets[SKYFIELD_EPHEMERIS_TARGETS[body_name]]).apparent()
|
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")
|
altitude, azimuth, _distance = apparent.altaz("standard")
|
||||||
return float(ra.hours), float(dec.degrees), float(azimuth.degrees), float(altitude.degrees)
|
return float(ra.hours), float(dec.degrees), float(azimuth.degrees), float(altitude.degrees)
|
||||||
|
|
||||||
|
|||||||
@@ -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()
|
||||||
Reference in New Issue
Block a user