astrololo/services/logic/app/engine/skyfield_engine.py

"""SkyfieldEngine — własny, permisywny silnik (LOG-01).

Ścieżka A: Skyfield (MIT) + efemerydy JPL (public domain). Liczy geocentryczne
pozycje pozorne (apparent) i rzutuje je na ekliptykę daty → długość tropikalna,
szerokość, prędkość i kierunek. Brak zależności AGPL.

Prędkość liczymy numerycznie (różnica długości po małym kroku czasu) — wystarcza
do kierunku (D/Rx) i do wykrywania stacji w LOG-03.
"""
from __future__ import annotations

import os
from datetime import timedelta
from functools import lru_cache

from app.engine.base import EphemerisEngine
from app.engine.formats import norm360
from app.engine.models import DEFAULT_OBJECTS, ChartMoment, ObjectPosition

# nazwa obiektu -> cel w jądrze efemeryd (de421 ma centra Merkurego/Wenus,
# dla pozostałych planet używamy barycentrów — różnica nieistotna astrologicznie)
_TARGETS = {
    "Sun": "sun",
    "Moon": "moon",
    "Mercury": "mercury",
    "Venus": "venus",
    "Mars": "mars barycenter",
    "Jupiter": "jupiter barycenter",
    "Saturn": "saturn barycenter",
    "Uranus": "uranus barycenter",
    "Neptune": "neptune barycenter",
    "Pluto": "pluto barycenter",
}


@lru_cache(maxsize=4)
def _load(kernel: str, data_dir: str):
    """Wczytuje skalę czasu i jądro efemeryd raz (kosztowne) i cache'uje."""
    from skyfield.api import Loader

    load = Loader(data_dir)
    ts = load.timescale()
    eph = load(kernel)
    return ts, eph, eph["earth"]


class SkyfieldEngine(EphemerisEngine):
    name = "skyfield"

    def __init__(self, kernel: str | None = None, data_dir: str | None = None) -> None:
        self.kernel = kernel or os.getenv("EPHEMERIS_KERNEL", "de421.bsp")
        self.data_dir = data_dir or os.getenv(
            "EPHEMERIS_DIR", os.path.join(os.path.dirname(__file__), "..", "..", ".ephemeris")
        )
        os.makedirs(self.data_dir, exist_ok=True)
        self.ts, self.eph, self.earth = _load(self.kernel, os.path.abspath(self.data_dir))

    def _ecliptic_lon_lat(self, target, t):
        astrometric = self.earth.at(t).observe(target).apparent()
        lat, lon, _dist = astrometric.ecliptic_latlon(epoch="date")
        return lon.degrees, lat.degrees

    def _virtual_point(self, name: str, tt_jd: float) -> ObjectPosition:
        """Punkty analityczne (LOG-02): mean Node (NN/SN) i mean Lilith.

        Liczone wzorami Meeusa, nie z jądra JPL. SN = NN + 180° (ta sama prędkość).
        Punkty leżą na ekliptyce (latitude = 0).
        """
        from app.engine.points import mean_lilith, mean_lunar_node, point_speed

        if name in ("North Node", "South Node"):
            lon = mean_lunar_node(tt_jd)
            if name == "South Node":
                lon = norm360(lon + 180.0)
            speed = point_speed(mean_lunar_node, tt_jd)
        else:  # Lilith
            lon = mean_lilith(tt_jd)
            speed = point_speed(mean_lilith, tt_jd)
        return ObjectPosition(
            name=name, longitude=float(lon), latitude=0.0,
            speed=float(speed), retrograde=bool(speed < 0),
        )

    def positions(
        self, moment: ChartMoment, objects: list[str] | None = None
    ) -> list[ObjectPosition]:
        names = objects or DEFAULT_OBJECTS
        t = self.ts.from_datetime(moment.when_utc)
        dt = timedelta(hours=1)
        t2 = self.ts.from_datetime(moment.when_utc + dt)

        out: list[ObjectPosition] = []
        for name in names:
            if name not in _TARGETS:                 # punkt wirtualny (NN/SN/Lilith)
                out.append(self._virtual_point(name, t.tt))
                continue
            target = self.eph[_TARGETS[name]]
            lon, lat = self._ecliptic_lon_lat(target, t)
            lon2, _ = self._ecliptic_lon_lat(target, t2)
            # prędkość °/dobę z poprawką na przejście przez 0°/360°
            step = ((lon2 - lon + 180.0) % 360.0) - 180.0
            speed = step * 24.0
            # rzutowanie na czysty float — Skyfield zwraca numpy.float64
            out.append(
                ObjectPosition(
                    name=name,
                    longitude=float(norm360(lon)),
                    latitude=float(lat),
                    speed=float(speed),
                    retrograde=bool(speed < 0),
                )
            )
        return out

    def sidereal(self, moment: ChartMoment) -> tuple[float, float]:
        """(RAMC, ε) w stopniach — lokalny apparent sidereal time i nachylenie ekliptyki.

        Materiał wejściowy do osi i domów (LOG-05). RAMC = GAST·15 + długość geo.
        """
        t = self.ts.from_datetime(moment.when_utc)
        ramc = norm360(t.gast * 15.0 + moment.lon)
        return ramc, true_obliquity(t)

    def health(self) -> dict:
        return {"engine": self.name, "status": "ok", "kernel": self.kernel}


def true_obliquity(t) -> float:
    """ε PRAWDZIWE (średnie + nutacja w nachyleniu), w stopniach.

    MUSI być prawdziwe, nie średnie. RAMC liczymy z `t.gast` — czasu gwiazdowego
    POZORNEGO, mierzonego od równonocy PRAWDZIWEJ. Ekliptyka odniesiona do tej
    samej równonocy ma ε z nutacją; sparowanie GAST z ε średnim miesza dwa układy
    odniesienia. Kosztowało to ~3,2″ na cuspach domów i było niewidoczne dla
    frameworka wyroczni, bo ten z założenia podaje to samo ε obu stronom
    (izoluje samą funkcję domów) — patrz tests/oracle/README.md.

    Nutacja z serii IAU 2000A, czyli z tego samego źródła, którego Skyfield
    używa do policzenia GAST — dzięki temu oba są spójne z definicji."""
    from skyfield.nutationlib import iau2000a, mean_obliquity

    mean_arcsec = float(mean_obliquity(t.tdb))
    d_eps_arcsec = float(iau2000a(t.tt)[1]) * 1e-7      # jednostki 0,1 µas
    return (mean_arcsec + d_eps_arcsec) / 3600.0