Ephemeriden mit Skyfield optimieren und dokumentieren

public/py/ephemeriden_api.py

@@ -12,6 +12,8 @@ if SCRIPT_DIR not in sys.path:
     sys.path.insert(0, SCRIPT_DIR)
 
 import astronomy
+from skyfield import almanac as skyfield_almanac
+from astronomical_conversions import skyfield_context, skyfield_wgs84
 
 
 EPHEMERIS_BODIES = {
@@ -26,6 +28,21 @@ EPHEMERIS_BODIES = {
     "Neptune": ("Neptun", astronomy.Body.Neptune),
 }
 
+# DE421 enthält die äußeren Planeten als Systemschwerpunkte, nicht als
+# separate Planetenzentren. Für Ephemeriden ist das die vorgesehene
+# Zielkörperdefinition dieser Datei.
+SKYFIELD_EPHEMERIS_TARGETS = {
+    "Sun": "sun",
+    "Moon": "moon",
+    "Mercury": "mercury",
+    "Venus": "venus",
+    "Mars": "mars",
+    "Jupiter": 5,
+    "Saturn": 6,
+    "Uranus": 7,
+    "Neptune": 8,
+}
+
 J2000_OBLIQUITY_DEG = 23.439279444444445
 GAUSSIAN_GRAVITATIONAL_CONSTANT = 0.01720209895
 PARABOLIC_ECCENTRICITY_TOLERANCE = 1.0e-6
@@ -118,6 +135,45 @@ def search_event(
     return result
 
 
+def skyfield_observer(latitude: float, longitude: float, elevation: float):
+    _timescale, planets = skyfield_context()
+    return planets["earth"] + skyfield_wgs84.latlon(latitude, longitude, elevation_m=elevation)
+
+
+def skyfield_planet_state(body_name: str, observer, dt_utc: datetime, timescale=None, planets=None):
+    if timescale is None or planets is None:
+        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)
+    altitude, azimuth, _distance = apparent.altaz("standard")
+    return float(ra.hours), float(dec.degrees), float(azimuth.degrees), float(altitude.degrees)
+
+
+def search_skyfield_rise_set(observer, body_name: str, local_day_start: datetime, timescale, planets) -> tuple[datetime | None, datetime | None]:
+    """Sucht Auf- und Untergang gemeinsam über Skyfields diskrete Ereignissuche."""
+    target = planets[SKYFIELD_EPHEMERIS_TARGETS[body_name]]
+    event_function = skyfield_almanac.risings_and_settings(
+        planets,
+        target,
+        observer - planets["earth"],
+    )
+    local_day_end = local_day_start + timedelta(days=1)
+    start_time = timescale.from_datetime(local_day_start.astimezone(timezone.utc))
+    end_time = timescale.from_datetime(local_day_end.astimezone(timezone.utc))
+    event_times, event_states = skyfield_almanac.find_discrete(start_time, end_time, event_function)
+
+    rise = None
+    set_ = None
+    for event_time, state in zip(event_times.utc_datetime(), event_states):
+        event_dt = event_time.astimezone(timezone.utc)
+        if bool(state) and rise is None:
+            rise = event_dt
+        elif not bool(state) and set_ is None:
+            set_ = event_dt
+    return rise, set_
+
+
 def format_ra_hours(ra_hours: float) -> str:
     total_seconds = int(round(float(ra_hours) * 3600.0))
     total_seconds %= 24 * 3600
@@ -522,21 +578,28 @@ def action_planet_ephemeris(args: list[str]) -> dict:
 
     label, body = EPHEMERIS_BODIES[body_name]
     observer = astronomy.Observer(latitude, longitude, elevation)
+    skyfield_timescale, skyfield_planets = skyfield_context()
+    skyfield_topos = skyfield_observer(latitude, longitude, elevation)
 
     def row_builder(current_local: datetime, time_value: astronomy.Time, rise_label: str | None, set_label: str | None) -> dict:
-        eq = astronomy.Equator(body, time_value, observer, False, True)
-        horizontal = astronomy.Horizon(time_value, observer, float(eq.ra), float(eq.dec), astronomy.Refraction.Normal)
+        ra_hours, dec_deg, azimuth_deg, altitude_deg = skyfield_planet_state(
+            body_name,
+            skyfield_topos,
+            current_local.astimezone(timezone.utc),
+            skyfield_timescale,
+            skyfield_planets,
+        )
         return {
             "object_name": label,
             "date_local": current_local.strftime("%d.%m.%Y"),
             "time_local": current_local.strftime("%H:%M"),
             "local_iso": current_local.isoformat(),
-            "ra": format_ra_hours(float(eq.ra)),
-            "ra_decimal_hours": round(float(eq.ra), 8),
-            "dec": format_dec_deg(float(eq.dec)),
-            "dec_decimal_deg": round(float(eq.dec), 8),
-            "azimuth_deg": round(float(horizontal.azimuth), 6),
-            "altitude_deg": round(float(horizontal.altitude), 6),
+            "ra": format_ra_hours(ra_hours),
+            "ra_decimal_hours": round(ra_hours, 8),
+            "dec": format_dec_deg(dec_deg),
+            "dec_decimal_deg": round(dec_deg, 8),
+            "azimuth_deg": round(azimuth_deg, 6),
+            "altitude_deg": round(altitude_deg, 6),
             "rise": rise_label,
             "set": set_label,
         }
@@ -547,13 +610,15 @@ def action_planet_ephemeris(args: list[str]) -> dict:
         except Exception as exc:
             fail("Zeitzone ist ungueltig.", extra={"details": str(exc)})
 
-        local_day_end = local_day_start + timedelta(days=1)
-        day_start_time = dt_to_time(local_day_start.astimezone(timezone.utc))
-        day_end_utc = local_day_end.astimezone(timezone.utc)
-        rise = search_event(body, astronomy.Direction.Rise, observer, day_start_time, day_end_utc)
-        set_ = search_event(body, astronomy.Direction.Set, observer, day_start_time, day_end_utc)
-        rise_label = serialize_event("Aufgang", rise, tz).get("local_time") if rise is not None else None
-        set_label = serialize_event("Untergang", set_, tz).get("local_time") if set_ is not None else None
+        rise, set_ = search_skyfield_rise_set(
+            skyfield_topos,
+            body_name,
+            local_day_start,
+            skyfield_timescale,
+            skyfield_planets,
+        )
+        rise_label = rise.astimezone(tz).strftime("%H:%M") if rise is not None else None
+        set_label = set_.astimezone(tz).strftime("%H:%M") if set_ is not None else None
         return rise_label, set_label
 
     return build_rows_response(