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

"""Systemy zodiaku (LOG-04): tropikalny, syderyczny (ayanamsy), draconic + RA.

Wszystkie pozycje silnika są liczone **tropikalnie of-date** (kontrakt LOG-28).
Zmiana zodiaku to — dla zodiaków ekliptycznych — jednolite przesunięcie długości:

    długość_docelowa = (długość_tropikalna − offset) mod 360

gdzie offset to:
  - **syderyczny**: ayanamsa (kąt między tropikalnym a syderycznym punktem Barana),
  - **draconic**: długość wznoszącego węzła Księżyca (węzeł = 0° draconic),
  - **tropikalny**: 0.

Ponieważ to stałe przesunięcie obiektu ORAZ cusps, **numery domów się nie zmieniają**
(geometria jest niezmiennicza względem obrotu) — przesuwamy tylko etykiety znaków.

Model ayanamsy: `ayan(jd) = ayan0 + B·x + C·x²`, gdzie `x = jd − J2000`. Prędkość
precesji (B, C) jest **wspólna** dla wszystkich ayanams; różni je tylko stała `ayan0`
(wybór syderycznego zera). Stałe skalibrowano do Swiss Ephemeris jako wyroczni —
zgodność do ~0,02" w latach 1900–2100 (patrz tests/test_zodiac.py).

RA (right ascension): konwersja ekliptyka→równik dla przyszłego widoku równikowego.
"""
from __future__ import annotations

import math
from datetime import datetime

from app.engine.formats import norm360

TROPICAL = "tropical"
DRACONIC = "draconic"

# stała ayanamsy w J2000.0 (°) — skalibrowana do swisseph (get_ayanamsa_ut)
_AYAN0 = {
    "lahiri": 23.857092,
    "fagan_bradley": 24.740300,
    "krishnamurti": 23.760240,
}
_J2000 = 2451545.0
_B = 3.824459e-5     # °/dobę — liniowy człon precesji (wspólny)
_C = 2.304e-13       # °/dobę² — drobne przyspieszenie (wspólne)

# nazwy zodiaków akceptowane przez API
SIDEREAL = tuple(f"sidereal_{k}" for k in _AYAN0)          # sidereal_lahiri, ...
SYSTEMS = (TROPICAL, *SIDEREAL, DRACONIC)


def julian_day(dt: datetime) -> float:
    """Julian Day (UT) z momentu UTC — algorytm Meeusa (kalendarz gregoriański)."""
    y, m = dt.year, dt.month
    day = dt.day + (dt.hour + dt.minute / 60.0 + dt.second / 3600.0
                    + dt.microsecond / 3.6e9) / 24.0
    if m <= 2:
        y -= 1
        m += 12
    a = y // 100
    b = 2 - a + a // 4
    return math.floor(365.25 * (y + 4716)) + math.floor(30.6001 * (m + 1)) + day + b - 1524.5


def ayanamsha(name: str, jd: float) -> float:
    """Ayanamsa [°] danej szkoły dla Julian Day (UT)."""
    key = name[len("sidereal_"):] if name.startswith("sidereal_") else name
    if key not in _AYAN0:
        raise ValueError(f"Nieznana ayanamsa: {name!r} (dostępne: {', '.join(_AYAN0)})")
    x = jd - _J2000
    return _AYAN0[key] + _B * x + _C * x * x


def offset(zodiac: str, jd: float, node_lon: float | None = None) -> float:
    """Ile odjąć od długości tropikalnej, by dostać wybrany zodiak.

    `node_lon` (tropikalna długość węzła wznoszącego) wymagana tylko dla draconic.
    """
    if zodiac == TROPICAL:
        return 0.0
    if zodiac == DRACONIC:
        if node_lon is None:
            raise ValueError("draconic wymaga długości węzła (node_lon)")
        return norm360(node_lon)
    if zodiac in SIDEREAL:
        return ayanamsha(zodiac, jd)
    raise ValueError(f"Nieznany zodiak: {zodiac!r} (dostępne: {', '.join(SYSTEMS)})")


def apply(lon: float, off: float) -> float:
    """Długość w docelowym zodiaku."""
    return norm360(lon - off)


def to_equatorial(lon: float, lat: float, eps: float) -> tuple[float, float]:
    """Ekliptyka (λ, β) → równik: (RA, deklinacja) w stopniach. Wszystko w °.

    RA rośnie 0–360°; deklinacja w [−90, 90].
    """
    lam, bet, e = math.radians(lon), math.radians(lat), math.radians(eps)
    sin_dec = math.sin(bet) * math.cos(e) + math.cos(bet) * math.sin(e) * math.sin(lam)
    dec = math.asin(max(-1.0, min(1.0, sin_dec)))
    ra = math.atan2(
        math.sin(lam) * math.cos(e) - math.tan(bet) * math.sin(e),
        math.cos(lam),
    )
    return norm360(math.degrees(ra)), math.degrees(dec)