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

"""Osie i domy — czysta matematyka sferyczna (LOG-05).

Bezstanowe funkcje: z lokalnego czasu gwiazdowego (RAMC), nachylenia ekliptyki
(ε) i szerokości geograficznej (φ) wyliczają Ascendent i MC, a stąd cusps domów
dla prostych systemów (Whole Sign, Equal, Porphyry). Niezależne od silnika —
silnik dostarcza tylko RAMC i ε.
"""
from __future__ import annotations

import math

from app.engine.formats import SIGN_ABBR, norm360, sign_index  # noqa: F401

WHOLE_SIGN = "whole_sign"
EQUAL = "equal"
PORPHYRY = "porphyry"
# Systemy o ZAMKNIĘTYM wzorze (bez iteracji). Placidus i Koch wymagają rozwiązania
# iteracyjnego i dochodzą osobno.
VEHLOW = "vehlow"
MORINUS = "morinus"
REGIOMONTANUS = "regiomontanus"
CAMPANUS = "campanus"
ALCABITUS = "alcabitus"
TOPOCENTRIC = "topocentric"
# Systemy WYPUSZCZONE — każdy zweryfikowany wobec Swiss Ephemeris co do zera
# (tests/oracle). TOPOCENTRIC celowo POZA listą: implementacja jest zgodna
# z wyrocznią wszędzie poza otoczeniem bieguna (|φ| ≈ 89,9° przy RAMC 90°/270°),
# gdzie kolejność domów się odwraca i konstrukcja traci sens — jego „biegun"
# atan(tan(89,9°)/3) to już 89,7°. Nie zawężamy dziedziny wyłącznie po to, by
# testy przeszły; system dołączy, gdy będzie poprawny na całej dziedzinie.
SYSTEMS = (WHOLE_SIGN, EQUAL, PORPHYRY, VEHLOW, MORINUS,
           REGIOMONTANUS, CAMPANUS, ALCABITUS)


def mean_obliquity(tt_jd: float) -> float:
    """Średnie nachylenie ekliptyki [°] dla daty (Julian TT). Wystarcza do domów."""
    t = (tt_jd - 2451545.0) / 36525.0
    arcsec = 84381.448 - 46.8150 * t - 0.00059 * t * t + 0.001813 * t ** 3
    return arcsec / 3600.0


def compute_mc(ramc_deg: float, eps_deg: float) -> float:
    r, e = math.radians(ramc_deg), math.radians(eps_deg)
    mc = math.atan2(math.sin(r), math.cos(r) * math.cos(e))
    return norm360(math.degrees(mc))


def compute_asc(ramc_deg: float, eps_deg: float, lat_deg: float) -> float:
    """Ascendent — punkt ekliptyki wschodzący na horyzoncie.

    KOREKTA GAŁĘZI (błąd wykryty przez porównanie z wyrocznią, tests/oracle):
    ekliptyka przecina horyzont w DWÓCH punktach — wschodzącym (Asc) i zachodzącym
    (Dsc). `atan2` wybiera jeden z nich, ale powyżej koła podbiegunowego potrafi
    wskazać ten NIEWŁAŚCIWY: dla szerokości 67°+ i szerokiego zakresu RAMC
    zwracaliśmy Descendent, czyli Ascendent przesunięty o 180°. Skutek nie był
    subtelny — planety lądowały w PRZECIWNYCH domach dla całej Skandynawii
    północnej (Tromsø, Rovaniemi, Murmańsk).

    Rozstrzyga położenie względem MC: punkt wschodzący leży zawsze w półkolu
    (0°, 180°) na wschód od MC. Reguła zweryfikowana na 46 800 przypadkach wobec
    Swiss Ephemeris — zero rozbieżności.
    """
    r, e, phi = math.radians(ramc_deg), math.radians(eps_deg), math.radians(lat_deg)
    asc = norm360(math.degrees(math.atan2(
        math.cos(r),
        -(math.sin(r) * math.cos(e) + math.tan(phi) * math.sin(e)),
    )))
    mc = compute_mc(ramc_deg, eps_deg)
    return norm360(asc + 180.0) if (asc - mc) % 360.0 > 180.0 else asc


def _trisect(a: float, b: float) -> tuple[float, float]:
    """Dwa punkty dzielące łuk a→b (w kierunku zodiaku) na trzy równe części."""
    span = (b - a) % 360.0
    return norm360(a + span / 3.0), norm360(a + 2.0 * span / 3.0)


# Poniżej tej odległości od granicy znaku traktujemy Ascendent jak leżący DOKŁADNIE
# na niej. 1e-9° to 3,6 mikrosekundy łuku — o rzędy wielkości poniżej jakiejkolwiek
# realnej dokładności danych urodzeniowych, więc nie zmienia to żadnego horoskopu.
_SIGN_SNAP_DEG = 1e-9


def _snap_to_sign_boundary(lon: float) -> float:
    """Przyciąga długość do granicy znaku, gdy jest od niej o włos.

    Whole sign jest NIECIĄGŁY na granicach znaków: różnica 10⁻¹¹° w Ascendencie
    przerzuca cały dom I o 30°. Bez tego przyciągania ten sam horoskop policzony
    na innej maszynie mógłby dać inny wynik (wykryte przez porównanie z wyrocznią:
    nasz Asc = 359,999999999976, swissepha = 1e-10 — ta sama wartość po dwóch
    stronach granicy Ryby/Baran). Determinizm jest tu ważniejszy niż dosłowność
    zmiennoprzecinkowa."""
    nearest = round(lon / 30.0) * 30.0
    return norm360(nearest) if abs(lon - nearest) < _SIGN_SNAP_DEG else lon


def cusps(asc: float, mc: float, system: str) -> list[float]:
    """Zwraca 12 cusps (długości) domów 1..12."""
    if system == WHOLE_SIGN:
        start = sign_index(_snap_to_sign_boundary(asc)) * 30.0
        return [norm360(start + 30.0 * i) for i in range(12)]
    if system == EQUAL:
        return [norm360(asc + 30.0 * i) for i in range(12)]
    if system == PORPHYRY:
        dsc, ic = norm360(asc + 180.0), norm360(mc + 180.0)
        c = [0.0] * 12
        c[0], c[3], c[6], c[9] = asc, ic, dsc, mc
        c[1], c[2] = _trisect(asc, ic)      # domy 2,3
        c[4], c[5] = _trisect(ic, dsc)      # domy 5,6
        c[7], c[8] = _trisect(dsc, mc)      # domy 8,9
        c[10], c[11] = _trisect(mc, asc)    # domy 11,12
        return c
    raise ValueError(f"nieznany system domów: {system}")


def polar_circle(eps_deg: float) -> float:
    """Szerokość koła podbiegunowego [°] dla danego nachylenia ekliptyki.

    NIE jest to stała 66,56°: ε zmienia się z datą (ok. 23,71° w 370 p.n.e.,
    23,44° dziś), więc granica przesuwa się o ~0,3° w zakresie dat programu.
    Powyżej niej stopnie ekliptyki bywają okołobiegunowe — nie wschodzą ani nie
    zachodzą — przez co systemy oparte na łuku dobowym (Placidus, Koch) tracą
    definicję."""
    return 90.0 - abs(eps_deg)


# ── geometria wektorowa dla systemów dzielących koła wielkie ─────────────
# Wzory na te systemy krążą w literaturze w kilku wariantach i łatwo o pomyłkę
# w gałęzi albo znaku. Liczymy więc WPROST z geometrii: budujemy wektory kierunkowe
# w układzie równikowym, przecinamy płaszczyzny i dopiero wynik zamieniamy na
# długość ekliptyczną. Jest to dłuższe, ale jednoznaczne i sprawdzalne.

def _cross(a, b):
    return (a[1] * b[2] - a[2] * b[1],
            a[2] * b[0] - a[0] * b[2],
            a[0] * b[1] - a[1] * b[0])


def _dot(a, b):
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]


def _equatorial_to_lon(v, eps_rad: float) -> float:
    """Wektor w układzie RÓWNIKOWYM → długość ekliptyczna [°]."""
    x, y, z = v
    y_ecl = y * math.cos(eps_rad) + z * math.sin(eps_rad)
    return norm360(math.degrees(math.atan2(y_ecl, x)))


def _ecliptic_pole(eps_rad: float):
    """Biegun ekliptyki (normalna płaszczyzny ekliptyki) w układzie równikowym."""
    return (0.0, -math.sin(eps_rad), math.cos(eps_rad))


def _horizon_north(ramc_rad: float, lat_rad: float):
    """Punkt północny horyzontu: RA = RAMC+180°, deklinacja = 90°−φ."""
    return (-math.sin(lat_rad) * math.cos(ramc_rad),
            -math.sin(lat_rad) * math.sin(ramc_rad),
            math.cos(lat_rad))


# Domy POŚREDNIE (bez osi) i to, po której stronie MC leżą. Domy 11, 12, 2, 3
# są na wschód od MC (przesunięcie 0–180°), domy 5, 6, 8, 9 — na zachód.
_INTERMEDIATE = {1: True, 2: True, 4: False, 5: False,
                 7: False, 8: False, 10: True, 11: True}


def _house_circle_cusp(north, q, eps_rad: float, mc: float, east_of_mc: bool) -> float:
    """Cusp = przecięcie ekliptyki z kołem domu.

    Koło domu przechodzi przez punkty N/S horyzontu oraz przez punkt podziału `q`
    (na równiku dla Regiomontanusa, na pierwszym wertykale dla Campanusa).
    Przecięcie dwóch płaszczyzn daje PROSTĄ, czyli DWA antypodyczne kierunki;
    wybieramy ten po właściwej stronie południka.

    Używane WYŁĄCZNIE dla domów pośrednich. Osie (1, 4, 7, 10) znamy dokładnie
    z Asc i MC — liczenie ich tą drogą było błędem, bo leżą dokładnie na granicy
    „wschód/zachód" (przesunięcie 0° i 180°), gdzie porównanie zmiennoprzecinkowe
    się chwieje i potrafi wybrać przeciwny punkt nieba."""
    normal = _cross(north, q)                 # normalna płaszczyzny koła domu
    line = _cross(normal, _ecliptic_pole(eps_rad))
    lon = _equatorial_to_lon(line, eps_rad)
    return lon if ((lon - mc) % 360.0 < 180.0) == east_of_mc else norm360(lon + 180.0)


def _culminating_mc(mc: float, eps: float, lat: float) -> float:
    """Punkt południka, który dla tej szerokości leży NAD horyzontem.

    Systemy oparte na horyzoncie (Regiomontanus, Campanus, Topocentric) biorą jako
    dziesiąty dom punkt GÓRUJĄCY, a nie matematyczne MC — a za kołem podbiegunowym
    to nie zawsze to samo. Punkt południka o deklinacji δ ma wysokość 90−|φ−δ|,
    więc jest nad horyzontem dokładnie wtedy, gdy |φ−δ| < 90.

    Systemy dzielące samą ekliptykę (porphyry, equal, alcabitus, whole sign) tego
    nie robią — i tak samo zachowuje się wyrocznia."""
    dec = math.degrees(math.asin(math.sin(math.radians(mc)) * math.sin(math.radians(eps))))
    return norm360(mc + 180.0) if abs(lat - dec) > 90.0 else mc


def _with_exact_angles(intermediate, asc: float, mc: float) -> list[float]:
    """Składa 12 cuspów: osie wstawione dokładnie, reszta z geometrii."""
    out = [0.0] * 12
    out[0], out[3] = asc, norm360(mc + 180.0)      # Asc, IC
    out[6], out[9] = norm360(asc + 180.0), mc      # Dsc, MC
    for i, value in intermediate.items():
        out[i] = value
    return out


def _ra_to_ecliptic_lon(ra_deg: float, eps_rad: float) -> float:
    """Punkt ekliptyki o zadanej rektascensji (koło godzinne → ekliptyka)."""
    r = math.radians(ra_deg)
    return norm360(math.degrees(math.atan2(math.sin(r), math.cos(r) * math.cos(eps_rad))))


def _equator_point(ra_deg: float):
    """Kierunek punktu na równiku niebieskim o danej rektascensji."""
    r = math.radians(ra_deg)
    return (math.cos(r), math.sin(r), 0.0)


def _prime_vertical_point(ramc_rad: float, lat_rad: float, angle_deg: float):
    """Punkt pierwszego wertykału, `angle_deg` od punktu wschodu w stronę nadiru.

    Pierwszy wertykał to koło przez wschód, zenit, zachód i nadir — Campanus dzieli
    właśnie je."""
    east = (-math.sin(ramc_rad), math.cos(ramc_rad), 0.0)
    zenith = (math.cos(lat_rad) * math.cos(ramc_rad),
              math.cos(lat_rad) * math.sin(ramc_rad),
              math.sin(lat_rad))
    a = math.radians(angle_deg)
    return tuple(east[i] * math.cos(a) - zenith[i] * math.sin(a) for i in range(3))


def _cusps_regiomontanus(ramc: float, eps: float, lat: float,
                         asc: float, mc: float) -> list[float]:
    """Równik niebieski dzielony na 12 równych łuków, rzut kołami przez N/S horyzontu."""
    er, rr, lr = math.radians(eps), math.radians(ramc), math.radians(lat)
    north = _horizon_north(rr, lr)
    mid = {i: _house_circle_cusp(north, _equator_point(ramc + 90.0 + 30.0 * i), er, mc, e)
           for i, e in _INTERMEDIATE.items()}
    return _with_exact_angles(mid, asc, _culminating_mc(mc, eps, lat))


def _cusps_campanus(ramc: float, eps: float, lat: float,
                    asc: float, mc: float) -> list[float]:
    """Pierwszy wertykał dzielony na 12 równych łuków, rzut tak samo jak wyżej."""
    er, rr, lr = math.radians(eps), math.radians(ramc), math.radians(lat)
    north = _horizon_north(rr, lr)
    mid = {i: _house_circle_cusp(north, _prime_vertical_point(rr, lr, 30.0 * i), er, mc, e)
           for i, e in _INTERMEDIATE.items()}
    return _with_exact_angles(mid, asc, _culminating_mc(mc, eps, lat))


def _cusps_morinus(ramc: float, eps: float) -> list[float]:
    """Równik dzielony od RAMC i rzutowany WPROST na ekliptykę — bez horyzontu.

    Dlatego Morinus jako jedyny nie zależy od szerokości geograficznej, a jego
    dom I nie pokrywa się z Ascendentem. Uwaga: to ZAMIANA WSPÓŁRZĘDNYCH punktu
    równika (RA, dec=0) na ekliptyczne, a nie rzut kołem godzinnym — te dwie
    operacje dają różne wyniki i pomylenie ich kosztowało tu do 5°."""
    er = math.radians(eps)
    return [_equatorial_to_lon(_equator_point(ramc + 90.0 + 30.0 * i), er)
            for i in range(12)]


def _cusps_alcabitus(ramc: float, eps: float, asc: float) -> list[float]:
    """Łuki równika MC→Asc i Asc→IC dzielone na trzy; rzut kołami godzinnymi."""
    er = math.radians(eps)
    a = math.radians(asc)
    ra_asc = norm360(math.degrees(math.atan2(math.sin(a) * math.cos(er), math.cos(a))))
    day = (ra_asc - ramc) % 360.0             # łuk MC → Asc (domy 11, 12)
    night = (ramc + 180.0 - ra_asc) % 360.0   # łuk Asc → IC (domy 2, 3)

    ra = [0.0] * 12
    ra[9] = ramc                              # dom 10 = MC
    ra[10] = ramc + day / 3.0                 # dom 11
    ra[11] = ramc + 2.0 * day / 3.0           # dom 12
    ra[0] = ra_asc                            # dom 1 = Asc
    ra[1] = ra_asc + night / 3.0              # dom 2
    ra[2] = ra_asc + 2.0 * night / 3.0        # dom 3
    for i in range(6):                        # domy 4–9 naprzeciw 10–3
        ra[i + 3] = ra[(i + 9) % 12] + 180.0
    return [_ra_to_ecliptic_lon(x, er) for x in ra]


# Ułamek szerokości geograficznej użyty jako „biegun" koła domu (Polich–Page).
# Domy na południku (10 i 4) mają biegun 0 — ich koło to sam południk.
_TOPO_POLE_FRACTION = (1.0, 2 / 3, 1 / 3, 0.0, 1 / 3, 2 / 3,
                       1.0, 2 / 3, 1 / 3, 0.0, 1 / 3, 2 / 3)


def _asc_under_pole_raw(ramc_deg: float, eps_deg: float, pole_deg: float) -> float:
    """Surowy wzór na ascendent pod zadanym „biegunem", BEZ korekty gałęzi.

    Korektę stosuje wywołujący — względem PRAWDZIWEGO MC horoskopu. `compute_asc`
    poprawia gałąź względem MC dla PRZESUNIĘTEGO RAMC, co dla cuspu domu jest złym
    punktem odniesienia i w okolicach biegunów dawało obrót o 180°."""
    r, e, phi = math.radians(ramc_deg), math.radians(eps_deg), math.radians(pole_deg)
    return norm360(math.degrees(math.atan2(
        math.cos(r), -(math.sin(r) * math.cos(e) + math.tan(phi) * math.sin(e)))))


def _cusps_topocentric(ramc: float, eps: float, lat: float,
                       asc: float, mc: float) -> list[float]:
    """Polich–Page: dom pośredni to ASCENDENT policzony pod własnym „biegunem"
    tan(P) = tan(φ)·k/3, dla RAMC przesuniętego o pozycję domu.

    Kusi, by liczyć to jak Regiomontanusa z podmienioną szerokością — daje wynik
    bliski, ale nie równy (kilka sekund łuku); wyrocznia rozstrzygnęła na rzecz
    konstrukcji „ascendent pod biegunem".

    Liczymy tylko domy 11, 12, 2, 3, a 5, 6, 8, 9 bierzemy jako ich OPOZYCJE —
    to nie skrót, lecz własność tych systemów: przeciwległe domy leżą na tym samym
    kole wielkim, więc ich cuspy są dokładnie antypodyczne."""
    tan_lat = math.tan(math.radians(lat))
    # Gałąź liczymy względem MC GÓRUJĄCEGO, nie matematycznego: gdy za kołem
    # podbiegunowym te dwa się rozjeżdżają, cała czwórka domów pośrednich musi
    # obrócić się razem z dziesiątym domem.
    culminating = _culminating_mc(mc, eps, lat)
    out = [0.0] * 12
    for i in (10, 11, 1, 2):                  # domy 11, 12, 2, 3
        pole = math.degrees(math.atan(tan_lat * _TOPO_POLE_FRACTION[i]))
        lon = _asc_under_pole_raw(ramc + 30.0 * i, eps, pole)
        if ((lon - culminating) % 360.0 < 180.0) != _INTERMEDIATE[i]:
            lon = norm360(lon + 180.0)
        out[i] = lon
        out[(i + 6) % 12] = norm360(lon + 180.0)
    out[0], out[3] = asc, norm360(culminating + 180.0)
    out[6], out[9] = norm360(asc + 180.0), culminating
    return out


def cusps_for(ramc: float, eps: float, lat: float, system: str) -> list[float]:
    """Kanoniczne wejście: (RAMC, ε, φ) → 12 cusps.

    Systemy proste (whole sign / equal / porphyry) potrzebują tylko Asc i MC,
    ale systemy egzotyczne dzielą inne koła wielkie i wymagają pełnego zestawu
    (RAMC, ε, φ). Ta funkcja jest wspólnym punktem wejścia dla obu rodzajów —
    i to ją porównuje z wyrocznią framework testowy (tests/oracle).
    """
    asc = compute_asc(ramc, eps, lat)
    mc = compute_mc(ramc, eps)
    if system in (WHOLE_SIGN, EQUAL, PORPHYRY):
        return cusps(asc, mc, system)
    if system == VEHLOW:
        # equal, ale Ascendent leży w ŚRODKU domu I, nie na jego początku
        return [norm360(asc - 15.0 + 30.0 * i) for i in range(12)]
    if system == MORINUS:
        return _cusps_morinus(ramc, eps)
    if system == REGIOMONTANUS:
        return _cusps_regiomontanus(ramc, eps, lat, asc, mc)
    if system == CAMPANUS:
        return _cusps_campanus(ramc, eps, lat, asc, mc)
    if system == ALCABITUS:
        return _cusps_alcabitus(ramc, eps, asc)
    if system == TOPOCENTRIC:
        return _cusps_topocentric(ramc, eps, lat, asc, mc)
    raise ValueError(f"nieznany system domów: {system}")


def assign_house(lon: float, cusp_list: list[float]) -> int:
    """Numer domu (1..12), w którym leży dana długość ekliptyczna."""
    lon = norm360(lon)
    for i in range(12):
        start = cusp_list[i]
        end = cusp_list[(i + 1) % 12]
        span = (end - start) % 360.0
        offset = (lon - start) % 360.0
        if offset < span:
            return i + 1
    return 12