Goldener henkel berechnung verbessert

public/py/api.py

@@ -494,19 +494,47 @@ def action_golden_handle_for_month(args: list[str]) -> dict:
     else:
         local_end = datetime(year, month + 1, 1, 0, 0, 0, tzinfo=tz)
 
-    # Approximate the Golden Handle geometrically by the effective selenographic
-    # colongitude over Montes Jura/Sinus Iridum (roughly 44 N, 36 W),
-    # corrected by libration. Empirically, the best-fit target shifts with
-    # libration in latitude by several tenths of a degree, i.e. by hours.
-    def target_colongitude(elat_deg: float) -> float:
-        return 33.385316705346064 + 0.11176238230626076 * elat_deg + 0.0734267792065295 * elat_deg * elat_deg
+    # Minimal physical model for the Golden Handle:
+    # solve the time when the floor of Sinus Iridum is still below sunrise
+    # by half of the topographic advance angle of the Jura rim.
+    moon_radius_km = 1737.4
+    floor_latitude_deg = 42.0
+    floor_longitude_west_deg = 35.5
+    relief_height_difference_km = 1.18
 
-    def golden_handle_metric(dt_utc: datetime) -> float:
+    def relief_advance_angle_deg() -> float:
+        return math.degrees(
+            math.acos(moon_radius_km / (moon_radius_km + relief_height_difference_km))
+        )
+
+    def solar_selenographic_coordinates(dt_utc: datetime) -> tuple[float, float]:
         time_value = dt_to_time(dt_utc)
         phase_deg = normalize_degrees(astronomy.MoonPhase(time_value))
         libration = astronomy.Libration(time_value)
-        effective_colongitude = normalize_degrees(phase_deg - 90.0 - float(libration.elon))
-        return normalize_signed_degrees(effective_colongitude - target_colongitude(float(libration.elat)))
+        axis_latitudes = calculate_moon_axis_latitudes(time_value)
+
+        # Approximate selenographic solar longitude from colongitude using
+        # west-positive longitude for the lunar surface model.
+        colongitude_deg = normalize_degrees(phase_deg - 90.0 - float(libration.elon))
+        solar_longitude_west_deg = normalize_signed_degrees(colongitude_deg - 90.0)
+        solar_latitude_deg = float(axis_latitudes["subsolar_latitude"])
+        return solar_longitude_west_deg, solar_latitude_deg
+
+    def floor_solar_altitude_deg(dt_utc: datetime) -> float:
+        solar_longitude_west_deg, solar_latitude_deg = solar_selenographic_coordinates(dt_utc)
+        phi = math.radians(floor_latitude_deg)
+        lam_f = math.radians(floor_longitude_west_deg)
+        lam_s = math.radians(solar_longitude_west_deg)
+        b_s = math.radians(solar_latitude_deg)
+        x = (
+            math.sin(phi) * math.sin(b_s)
+            + math.cos(phi) * math.cos(b_s) * math.cos(lam_f - lam_s)
+        )
+        x = max(-1.0, min(1.0, x))
+        return math.degrees(math.asin(x))
+
+    def golden_handle_metric(dt_utc: datetime) -> float:
+        return floor_solar_altitude_deg(dt_utc) + relief_advance_angle_deg() / 2.0
 
     def refine_peak_time(left_utc: datetime, right_utc: datetime) -> datetime:
         left = left_utc