fix(ephemeris): include nutation in the sidereal frame

恒星黄道改为与 Swiss Ephemeris FLG_SIDEREAL 同一套含章动岁差。产品默认路径不再用平岁差去减含章动的视位置。上升和宫头用同一个岁差。校正身份升到 rectification-v5-matrix-scoring-12。 BUG-1275。未推送，未部署。v5 先验头名三格下降 1.3 个百分点，不调参。
2026-10-08 18:42:45 +08:00
parent 1a8c97f64a
commit e5f2a75183
37 changed files with 52297 additions and 113 deletions
@@ -132,3 +132,20 @@ def temporary_ayanamsa(name: Any = None, swe_module: Any = None) -> Iterator[str
 def sidereal_flags(swe_module: Any, ayanamsa_name: Any = None) -> int:
    apply_ayanamsa(ayanamsa_name, swe_module)
    return getattr(swe_module, 'FLG_SWIEPH', 2) | getattr(swe_module, 'FLG_SIDEREAL', 65536)
+
+
+def sidereal_ayanamsa_degrees(jd_ut: float, swe_module: Any = None, flags: int | None = None) -> float:
+    """Ayanamsa that FLG_SIDEREAL subtracts for these flags.
+
+    ``get_ayanamsa`` and ``get_ayanamsa_ut`` omit nutation. Apparent tropical
+    positions include it, so that subtraction shifts every longitude by Δψ.
+    """
+    module = swe_module or swe
+    if module is None:
+        raise RuntimeError("swisseph is not available")
+    if flags is None:
+        flags = getattr(module, "FLG_SWIEPH", 2)
+    raw = module.get_ayanamsa_ex_ut(float(jd_ut), int(flags))
+    if isinstance(raw, (tuple, list)):
+        return float(raw[1])
+    return float(raw)
@@ -219,7 +219,8 @@ class BhavaChalitCalculator:
        hsys = self._SWE_HSYS.get(house_system, b'P')
        cusps_swe, ascmc = self._swe.houses(jd, lat, lon, hsys)

-        ayanamsa = self._swe.get_ayanamsa(jd)
+        from ayanamsa_utils import sidereal_ayanamsa_degrees
+        ayanamsa = sidereal_ayanamsa_degrees(jd, self._swe)
        result = []
        for i in range(12):
            result.append(_norm(cusps_swe[i] - ayanamsa))
@@ -13,11 +13,11 @@ if str(SCRIPT_DIR) not in sys.path:

 try:
    import swisseph as swe
-    from ayanamsa_utils import apply_ayanamsa, normalize_ayanamsa_name
+    from ayanamsa_utils import apply_ayanamsa, normalize_ayanamsa_name, sidereal_ayanamsa_degrees
    from domain_calculation_service import compute_chart, compute_vimshottari_timeline
 except ModuleNotFoundError:
    import swisseph as swe
-    from scripts.ayanamsa_utils import apply_ayanamsa, normalize_ayanamsa_name
+    from scripts.ayanamsa_utils import apply_ayanamsa, normalize_ayanamsa_name, sidereal_ayanamsa_degrees
    from scripts.domain_calculation_service import compute_chart, compute_vimshottari_timeline

 SIGNS = [
@@ -65,7 +65,7 @@ def _moon_transit(reference_date: str, tz: float, ayanamsa: str) -> dict[str, An
    local_dt = datetime.strptime(reference_date[:10], "%Y-%m-%d").replace(hour=12)
    apply_ayanamsa(normalize_ayanamsa_name(ayanamsa), swe)
    jd = swe.julday(local_dt.year, local_dt.month, local_dt.day, 12.0 - float(tz))
-    ayanamsa_value = swe.get_ayanamsa(jd)
+    ayanamsa_value = sidereal_ayanamsa_degrees(jd, swe)
    position, _flags = swe.calc_ut(jd, swe.MOON)
    longitude = (position[0] - ayanamsa_value) % 360
    sign_index = int(longitude // 30)
@@ -18,6 +18,7 @@ from ayanamsa_utils import (
    apply_ayanamsa,
    is_supported_ayanamsa_name,
    normalize_ayanamsa_name,
+    sidereal_ayanamsa_degrees,
    supported_ayanamsa_names,
 )
 from dasha_analyzer import build_dasha_timeline, lon_to_nakshatra
@@ -300,7 +301,7 @@ def compute_transit_longitude(
            local_dt.day,
            12.0 - float(tz),
        )
-        ayanamsa_value = swe.get_ayanamsa(jd)
+        ayanamsa_value = sidereal_ayanamsa_degrees(jd, swe)
        position, flags = swe.calc_ut(jd, _PLANET_IDS[planet])
    longitude = (position[0] - ayanamsa_value) % 360
    return {
@@ -97,6 +97,7 @@ from ayanamsa_utils import (
    ayanamsa_already_applied,
    ayanamsa_display_name,
    current_ayanamsa_name,
+    sidereal_ayanamsa_degrees,
 )

 # ============================================================================
@@ -816,16 +817,17 @@ def _normalize_position_mode(position_mode):


 def _sidereal_calculation_profile(jd, position_mode):
-    """Return matching Swiss Ephemeris flags and ayanamsa for one frame."""
+    """Return matching Swiss Ephemeris flags and ayanamsa for one frame.
+
+    ``legacy`` used to take ``get_ayanamsa`` (no nutation) against apparent
+    positions. That was the Δψ shift. It now uses the same true ayanamsa as
+    ``apparent``. ``mean`` still drops nutation on both sides. No fourth mode.
+    """
    mode = _normalize_position_mode(position_mode)
    ephemeris_flags = getattr(swe, 'FLG_SWIEPH', 2)
    if mode == 'mean':
        ephemeris_flags |= getattr(swe, 'FLG_NONUT', 64)
-        ayanamsa = swe.get_ayanamsa_ex_ut(jd, ephemeris_flags)[1]
-    elif mode == 'apparent':
-        ayanamsa = swe.get_ayanamsa_ex_ut(jd, ephemeris_flags)[1]
-    else:
-        ayanamsa = swe.get_ayanamsa(jd)
+    ayanamsa = sidereal_ayanamsa_degrees(jd, swe, ephemeris_flags)
    return mode, ephemeris_flags, ayanamsa


@@ -838,8 +840,7 @@ def _calc_sidereal_planets_for_jd(jd, node_mode='mean', include_ketu=True, ayana
        return {}, None
    if ayanamsa_name:
        _apply_ayanamsa(ayanamsa_name)
-    ayanamsa = swe.get_ayanamsa(jd)
-    calc_flags = getattr(swe, 'FLG_SWIEPH', 2)
+    _mode, calc_flags, ayanamsa = _sidereal_calculation_profile(jd, 'legacy')
    data = {}
    for pname, pid in _planet_map_for_node_mode(node_mode).items():
        pos, ret = swe.calc_ut(jd, pid, calc_flags | getattr(swe, 'FLG_SPEED', 256))
@@ -12347,7 +12348,7 @@ def compute_chart_data(year, month, day, hour, minute, lat, lon, tz, node_mode='
    second = int(second or 0)
    hour_decimal = _birth_hour_decimal(hour, minute, second) - tz
    jd = swe.julday(year, month, day, hour_decimal)
-    ayanamsa = swe.get_ayanamsa(jd)
+    _mode, calc_flags, ayanamsa = _sidereal_calculation_profile(jd, 'legacy')

    node_mode = (node_mode or 'mean').lower()
    if node_mode not in ('mean', 'true'):
@@ -12384,7 +12385,6 @@ def compute_chart_data(year, month, day, hour, minute, lat, lon, tz, node_mode='
            "lord": SIGN_LORDS[SIGNS[si]]}

    planets_swe = {**BASE_PLANETS_SWE, 'Rahu': node_pid}
-    calc_flags = getattr(swe, 'FLG_SWIEPH', 2)
    for pname, pid in planets_swe.items():
        try:
            pos, _ret_flags = swe.calc_ut(jd, pid, calc_flags | getattr(swe, 'FLG_SPEED', 256))
@@ -13003,8 +13003,6 @@ def cmd_transit(args):
    except Exception as e:
        return {"error": f"日期计算失败: {e}"}

-    ayanamsa = swe.get_ayanamsa(jd_mid)
-
    # --- 计算行星位置 ---
    node_mode = getattr(args, 'node_mode', 'mean')
    transit_data, ayanamsa = _calc_sidereal_planets_for_jd(
@@ -13368,12 +13366,12 @@ def cmd_double_transit_pac(args):
    transit_year, transit_month, transit_day = map(int, args.date.split('-'))
    transit_hour = 12.0 - args.tz  # 正午 UT
    transit_jd = swe.julday(transit_year, transit_month, transit_day, transit_hour)
-    transit_ayanamsa = swe.get_ayanamsa(transit_jd)
+    _mode, transit_flags, transit_ayanamsa = _sidereal_calculation_profile(transit_jd, 'legacy')

    transit_planets = {}
    for pname, pid in PLANETS_SWE.items():
        try:
-            pos, _ = swe.calc_ut(transit_jd, pid)
+            pos, _ = swe.calc_ut(transit_jd, pid, transit_flags)
            lon_p = (pos[0] - transit_ayanamsa) % 360
            transit_planets[pname] = {'lon': lon_p, 'sign': SIGNS[int(lon_p / 30)]}
            if pname == 'Rahu':
@@ -13579,8 +13577,8 @@ def _calc_transit_lon(jd, planet_name):
    pid = pid_map.get(planet_name)
    if pid is None:
        return None
-    pos, _ = swe.calc_ut(jd, pid)
-    aya = swe.get_ayanamsa_ut(jd)
+    _mode, flags, aya = _sidereal_calculation_profile(jd, 'legacy')
+    pos, _ = swe.calc_ut(jd, pid, flags)
    return (pos[0] - aya) % 360


@@ -13617,7 +13615,7 @@ def cmd_transit_ll7l(args):
        'transit_ayanamsa_policy': _transit_ayanamsa_mixed_review_policy(
            'transit-ll7l',
            requested_ayanamsa=getattr(args, 'ayanamsa', None),
-            transit_ayanamsa=swe.get_ayanamsa_ut(transit_jd),
+            transit_ayanamsa=_sidereal_calculation_profile(transit_jd, 'legacy')[2],
        ),
        'lagna_lord': ll_name,
        'seventh_lord': seven_lord,
@@ -13711,7 +13709,7 @@ def cmd_planetary_congregation(args):
    if args.transit_date:
        t_year, t_month, t_day = map(int, args.transit_date.split('-'))
        transit_jd = swe.julday(t_year, t_month, t_day, 12.0 - args.tz)
-        transit_aya = swe.get_ayanamsa(transit_jd)
+        _mode, transit_flags, transit_aya = _sidereal_calculation_profile(transit_jd, 'legacy')
        result['transit_ayanamsa_policy'] = _transit_ayanamsa_mixed_review_policy(
            'planetary-congregation',
            requested_ayanamsa=getattr(args, 'ayanamsa', None),
@@ -13720,7 +13718,7 @@ def cmd_planetary_congregation(args):
        result['transit'] = {str(h): [] for h in range(1, 13)}
        for pname, pid in PLANETS_SWE.items():
            try:
-                pos, _ = swe.calc_ut(transit_jd, pid)
+                pos, _ = swe.calc_ut(transit_jd, pid, transit_flags)
                lon_p = (pos[0] - transit_aya) % 360
                si = int(lon_p / 30)
                house = ((si - asc_idx) % 12) + 1
@@ -13802,7 +13800,7 @@ def cmd_vivah_saham(args):
        result['transit_ayanamsa_policy'] = _transit_ayanamsa_mixed_review_policy(
            'vivah-saham',
            requested_ayanamsa=getattr(args, 'ayanamsa', None),
-            transit_ayanamsa=swe.get_ayanamsa_ut(transit_jd),
+            transit_ayanamsa=_sidereal_calculation_profile(transit_jd, 'legacy')[2],
        )

        result['transit_activation'] = {'jupiter': [], 'saturn': [], 'double_activation': False}
@@ -9,7 +9,7 @@ from __future__ import annotations

 from typing import Any

-from ayanamsa_utils import DEFAULT_AYANAMSA_NAME, apply_ayanamsa, current_ayanamsa_name
+from ayanamsa_utils import DEFAULT_AYANAMSA_NAME, apply_ayanamsa, current_ayanamsa_name, sidereal_ayanamsa_degrees
 from kp_system import KP_LORDS, get_kp_lords
 from scripts.domain_calculation_service import swiss_ephemeris_lock

@@ -59,7 +59,7 @@ def observe_kp_cusps(jd: float | None, lat: float | None, lon: float | None) ->
        try:
            apply_ayanamsa("krishnamurti", swe)
            cusps_raw, _ascmc = swe.houses(float(jd), float(lat), float(lon), b"P")
-            ayanamsa = float(swe.get_ayanamsa(float(jd)))
+            ayanamsa = sidereal_ayanamsa_degrees(float(jd), swe)
            houses: dict[str, dict[str, Any]] = {}
            indices: dict[str, int] = {}
            for house in range(1, 13):
@@ -21,7 +21,7 @@ from scripts.rectification.case_holdout import holdout_event_ids
 # scoring-11 (2026-10-07, BUG-1260): D5/D11 sign mappings change education and
 # finance varga points for the same input. Cached scoring-10 results must not
 # be reused as current. Policy and the dated window contract stay as they are.
-ALGORITHM_VERSION = "rectification-v5-matrix-scoring-11"
+ALGORITHM_VERSION = "rectification-v5-matrix-scoring-12"
 INPUT_CONTRACT_VERSION = "rectification-calculation-spec-v5"
 PRECISION_WEIGHTS = {
    "day": 1.0,
@@ -0,0 +1,102 @@
+#!/usr/bin/env python3
+"""Freeze PyJHora 4.8.7 sidereal longitudes. Offline only.
+
+Run from the repo root with a Python that has PyJHora 4.8.7:
+
+    python scripts/research/freeze_sidereal_nutation_pyjhora_fixture.py
+
+Tests read the JSON and do not import this module.
+"""
+
+from __future__ import annotations
+
+import importlib.metadata
+import io
+import json
+import sys
+from contextlib import redirect_stdout
+from pathlib import Path
+
+ROOT = Path(__file__).resolve().parents[2]
+sys.path[:0] = [str(ROOT), str(ROOT / "scripts")]
+
+import swisseph as swe  # noqa: E402
+from jhora import const  # noqa: E402
+from jhora.panchanga import drik  # noqa: E402
+from jyotish_engine import compute_chart_data  # noqa: E402
+
+OUT = ROOT / "tests" / "fixtures" / "sidereal_nutation_pyjhora_487.json"
+MODES = {"lahiri": "LAHIRI", "raman": "RAMAN", "kp": "KP"}
+BODIES = (
+    ("Sun", const._SUN), ("Moon", const._MOON), ("Mars", const._MARS),
+    ("Mercury", const._MERCURY), ("Jupiter", const._JUPITER), ("Venus", const._VENUS),
+    ("Saturn", const._SATURN), ("Rahu", const._RAHU),
+)
+
+
+def _births() -> dict[str, dict]:
+    public = json.loads((ROOT / "references" / "public_oracle_cases.json").read_text(encoding="utf-8"))
+    births = {
+        case["id"]: case["birth"]
+        for case in public["cases"]
+        if case["id"] in {"steve_jobs_1955_aa", "barack_obama_1961_aa", "einstein_1879_aa"}
+    }
+    special = json.loads((ROOT / "tests" / "fixtures" / "special_lagnas_pyjhora_487.json").read_text(encoding="utf-8"))
+    births["fictional_reader_main"] = special["cases"]["fictional_reader_main"]["birth"]
+    return births
+
+
+def _longitudes(birth: dict, mode: str) -> dict[str, float]:
+    """Same UTC conversion as drik.planetary_positions, via sidereal_longitude.
+
+    PyJHora 4.8.7's planetary_positions calls dict.index on planet_list and
+    raises AttributeError, so the wrapper cannot be the frozen source.
+    """
+    local_hour = birth["hour"] + birth["minute"] / 60.0 + birth.get("second", 0) / 3600.0
+    jd_local = swe.julday(birth["year"], birth["month"], birth["day"], local_hour)
+    jd_ut = jd_local - float(birth["tz"]) / 24.0
+    with redirect_stdout(io.StringIO()):
+        drik.set_ayanamsa_mode(mode)
+        longitudes = {
+            name: float(drik.sidereal_longitude(jd_ut, planet)) % 360.0
+            for name, planet in BODIES
+        }
+        longitudes["Ketu"] = float(drik.ketu(longitudes["Rahu"])) % 360.0
+    return longitudes
+
+
+def _gap(left: float, right: float) -> float:
+    return abs(((left - right + 180.0) % 360.0) - 180.0)
+
+
+def main() -> None:
+    cases = {}
+    worst = 0.0
+    for case_id, birth in _births().items():
+        by_mode = {}
+        for name, mode in MODES.items():
+            frozen = _longitudes(birth, mode)
+            chart, *_rest = compute_chart_data(
+                birth["year"], birth["month"], birth["day"], birth["hour"], birth["minute"],
+                birth["lat"], birth["lon"], birth["tz"], node_mode="true",
+                second=birth.get("second", 0), ayanamsa_name=name,
+            )
+            for planet, longitude in frozen.items():
+                worst = max(worst, _gap(chart["planets"][planet]["degree_raw"], longitude))
+            by_mode[name] = frozen
+        cases[case_id] = {"ayanamsa": by_mode}
+    payload = {
+        "pyjhora_version": importlib.metadata.version("PyJHora"),
+        "generator": "python scripts/research/freeze_sidereal_nutation_pyjhora_fixture.py",
+        "generated_on": "2026-10-08",
+        "node_mode": "true",
+        "max_abs_difference_degrees": worst,
+        "note": "Frozen from drik.sidereal_longitude. PyJHora 4.8.7 planetary_positions raises AttributeError because planet_list is a dict. Rahu body is const._RAHU (Swiss 11). Differences are recorded, not a second runtime.",
+        "cases": cases,
+    }
+    OUT.write_text(json.dumps(payload, ensure_ascii=False, indent=2) + "\n", encoding="utf-8", newline="\n")
+    print(OUT.relative_to(ROOT).as_posix())
+
+
+if __name__ == "__main__":
+    main()
@@ -30,8 +30,8 @@ from scripts.research.sealed_holdout_rerun import (
 OFFSETS = (-30, -20, -15, -10, -8, -5, -3, 0, 3, 5, 8, 10, 15, 20, 30)
 RADII = (15, 30, 60)
 MINUTE_STEP = 1
-FREEZE = ROOT / "docs/research/reported_offset_upstream_sync6_2026_10_07.freeze.json"
-REPORT = ROOT / "docs/research/reported_offset_upstream_sync6_2026_10_07.json"
+FREEZE = ROOT / "docs/research/reported_offset_ephemeris_nutation_2026_10_08.freeze.json"
+REPORT = ROOT / "docs/research/reported_offset_ephemeris_nutation_2026_10_08.json"
 LEGACY_REPORT = ROOT / "docs/research/reported_offset_2026_09_20.json"


@@ -27,8 +27,8 @@ from scripts.minute_rectification_feature_facts_v4 import build_feature_fact_row
 from scripts.minute_rectification_holdout_validator import validate

 DATASET = ROOT / "references/real_case_calibration/minute_rectification_holdout_v3.json"
-FREEZE = ROOT / "docs/research/sealed_holdout_rerun_upstream_sync6_2026_10_07.freeze.json"
-REPORT = ROOT / "docs/research/sealed_holdout_rerun_upstream_sync6_2026_10_07.json"
+FREEZE = ROOT / "docs/research/sealed_holdout_rerun_ephemeris_nutation_2026_10_08.freeze.json"
+REPORT = ROOT / "docs/research/sealed_holdout_rerun_ephemeris_nutation_2026_10_08.json"
 LEGACY_REPORT = ROOT / "docs/research/sealed_holdout_rerun_2026_09_20.json"
 ARCHIVE = ROOT / "docs/research/history/rectification_pre_cross_midnight_2026_09_20"
 PRODUCTION_FILES = [
@@ -32,7 +32,7 @@ _script_dir = os.path.dirname(os.path.abspath(__file__))
 if _script_dir not in sys.path:
    sys.path.insert(0, _script_dir)
 from varga import varga_map, SIGNS as VARGA_SIGNS, SIGN_LORDS as VARGA_SIGN_LORDS
-from ayanamsa_utils import DEFAULT_AYANAMSA_NAME, apply_ayanamsa, normalize_ayanamsa_name
+from ayanamsa_utils import DEFAULT_AYANAMSA_NAME, apply_ayanamsa, normalize_ayanamsa_name, sidereal_ayanamsa_degrees

 # ============================================================================
 # 常量
@@ -173,7 +173,7 @@ def build_shadbala_context(jd_ut: float, lat: float, lon: float, ayanamsa: str =
        "jd_ut": float(jd_ut), "jd_local": local_jd,
        "lat": float(lat), "lon": float(lon), "timezone": float(timezone),
        "year": int(year), "month": int(month), "day": int(day), "local_hour": float(local_hour),
-        "ayanamsa": key, "ayanamsa_degrees": float(swe.get_ayanamsa_ut(float(jd_ut))),
+        "ayanamsa": key, "ayanamsa_degrees": sidereal_ayanamsa_degrees(float(jd_ut), swe),
        "sunrise_hour": sunrise, "sunset_hour": sunset, "previous_sunset_hour": previous_sunset,
        "house_midpoints": [float(value) % 360 for value in cusps],
    }
@@ -23,7 +23,7 @@ import traceback
 from pathlib import Path
 from typing import Dict, List, Optional, Tuple, Set

-from ayanamsa_utils import sidereal_flags
+from ayanamsa_utils import sidereal_ayanamsa_degrees, sidereal_flags

 # ============================================================
 # 路径设置
@@ -103,7 +103,7 @@ def skill_compute_chart(year, month, day, hour, minute, lat, lon, tz,

    # Ayanamsa
    flags = sidereal_flags(swe, 'lahiri')
-    ayanamsa = swe.get_ayanamsa_ut(jd)
+    ayanamsa = sidereal_ayanamsa_degrees(jd, swe)

    # Ascendant
    asc_info = swe.houses(jd, lat, lon, b'A')  # 'A' = equal houses (Vedic)