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 dataclasses import dataclass

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

WHOLE_SIGN = "whole_sign"
EQUAL = "equal"
EQUAL_MC = "equal_mc"                 # równe domy zakotwiczone na MC, nie na Asc
WHOLE_SIGN_ARIES = "whole_sign_aries"  # znaki jako domy, ale dom I to ZAWSZE Baran
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"
PLACIDUS = "placidus"
KOCH = "koch"
# Systemy WYPUSZCZONE — każdy zweryfikowany wobec Swiss Ephemeris
# (tests/oracle). Placidus i Koch jako jedyne mają granicę dziedziny:
# powyżej koła podbiegunowego nie istnieją i podlegają jawnemu fallbackowi.
SYSTEMS = (WHOLE_SIGN, WHOLE_SIGN_ARIES, EQUAL, EQUAL_MC, PORPHYRY, VEHLOW,
           MORINUS, REGIOMONTANUS, CAMPANUS, ALCABITUS, TOPOCENTRIC,
           PLACIDUS, KOCH)


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.
# Polich–Page: dom pośredni to Ascendent policzony pod własnym „biegunem"
# tan(P) = tan(φ)·k/3, dla RAMC przesuniętego o pozycję domu. Rodzina jest CIĄGŁA:
# przy k=0 (biegun 0, przesunięcie −90°) daje MC, przy k=3 (biegun φ, przesunięcie 0)
# Ascendent, a domy 11 i 12 leżą po drodze.
#
# Cała trudność tego systemu siedziała w wyborze gałęzi — dwa koła wielkie
# przecinają się w dwóch punktach antypodycznych. Heurystyki („po której stronie
# MC", „w łuku kwadrantu", „wschodnia połowa horyzontu", śledzenie ciągłości
# krokami) myliły się na 6–11% przypadków powyżej ~70°, bo każda z nich rozstrzyga
# LOKALNIE, a przy dużych szerokościach kolejność domów potrafi się odwrócić.
#
# Rozwiązanie: nie wybierać w ogóle. Iloczyn wektorowy zenitu z biegunem ekliptyki
# jest ciągłą funkcją parametru rodziny i sam niesie właściwy zwrot — dwuznaczność
# wprowadza dopiero atan2. Zostajemy więc w wektorach, a znak ustalamy RAZ, kotwicząc
# rodzinę na MC górującym. Stąd zgodność co do zera na całej dziedzinie, bez iteracji
# i bez zawężania szerokości.

# Dom → (przesunięcie RAMC [°], ułamek bieguna k/3).
_TOPO_STEP = {10: (-60.0, 1 / 3), 11: (-30.0, 2 / 3),    # domy 11, 12
              1: (30.0, 2 / 3), 2: (60.0, 1 / 3)}        # domy 2, 3


def _cusps_topocentric(ramc: float, eps: float, lat: float,
                       asc: float, mc: float) -> list[float]:
    """Ascendenty pod biegunami tan(P) = tan(φ)·k/3, liczone wektorowo.

    Domy 5, 6, 8, 9 bierzemy jako OPOZYCJE domów 11, 12, 2, 3 — to nie skrót,
    lecz własność konstrukcji: przeciwległe domy leżą na tym samym kole wielkim.

    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"."""
    # Cuspy topocentrica są poprawne na CAŁEJ dziedzinie (zgodne z wyrocznią co do
    # zera), ale powyżej koła podbiegunowego przestają DZIELIĆ OKRĄG: domy nachodzą
    # na siebie, bo cusp VII (= I + 180°) wypada przed cuspem VI. Przypisanie planety
    # do domu traci wtedy sens — co potwierdza sama wyrocznia, której swe_house_pos
    # przeczy tam własnym cuspom (100% zgodności do 62°, 83,9% przy 66°, ok. 50%
    # przy 72°). Odmawiamy, zamiast zwracać liczbę bez znaczenia.
    #
    # To INNY rodzaj granicy niż u Placidusa i Kocha: tam nie istnieją same cuspy,
    # tu istnieją, tylko nie tworzą podziału. Próg jest wyprowadzony z warunku
    # „dwanaście cuspów sumuje się do 360°", nie dobrany pod wynik testu — i wypada
    # na kole podbiegunowym (zmierzone: 100% podziałów do 65°, 78% w pasie 66-67°).
    if abs(lat) >= polar_circle(eps):
        raise HouseSystemUndefined(
            f"φ={lat:.4f}° powyżej koła podbiegunowego ({polar_circle(eps):.4f}° dla "
            f"ε={eps:.4f}°): cuspy topocentryczne przestają dzielić okrąg, "
            f"domy nachodzą na siebie")
    er, tan_lat = math.radians(eps), math.tan(math.radians(lat))
    epole = _ecliptic_pole(er)
    # Gdy MC górujące rozjeżdża się z matematycznym (za kołem podbiegunowym),
    # cała rodzina obraca się razem z dziesiątym domem — stąd zwrot iloczynu.
    culminating = _culminating_mc(mc, eps, lat)
    sign = 1.0 if abs(((culminating - mc + 180.0) % 360.0) - 180.0) > 90.0 else -1.0

    out = [0.0] * 12
    for i, (offset, fraction) in _TOPO_STEP.items():
        th = math.radians(ramc + offset)
        pole = math.atan(tan_lat * fraction)
        zenith = (math.cos(pole) * math.cos(th),
                  math.cos(pole) * math.sin(th),
                  math.sin(pole))
        v = _cross(zenith, epole)
        lon = _equatorial_to_lon(tuple(sign * x for x in v), er)
        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


# ── systemy łuku dobowego (Placidus, Koch) ───────────────────────────────
# Różnią się od wszystkich poprzednich tym, że NIE MAJĄ wzoru zamkniętego: cusp
# jest zdefiniowany warunkiem na samego siebie („punkt, który przebył 1/3 swojego
# półłuku"), więc trzeba go znaleźć iteracyjnie. Mają też jako jedyne REALNĄ
# granicę dziedziny — powyżej koła podbiegunowego stopnie ekliptyki bywają
# okołobiegunowe, nie wschodzą ani nie zachodzą, i półłuk po prostu nie istnieje.


class HouseSystemUndefined(ValueError):
    """System domów nie ma definicji dla podanych parametrów (nie: błąd liczenia).

    Podnoszone zamiast zwrócenia liczby, bo cicha podmiana systemu jest gorsza
    niż błąd: wykres wygląda poprawnie, a planety siedzą w innych domach, niż
    astrolog zamawiał. Warstwa aplikacyjna łapie to w cusps_detailed() i robi
    JAWNY fallback."""


_ITER_MAX = 100
_ITER_TOL_DEG = 1e-11


def _declination_of_ecliptic_lon(lon_deg: float, eps_rad: float) -> float:
    """Deklinacja punktu LEŻĄCEGO NA EKLIPTYCE o danej długości."""
    return math.degrees(math.asin(math.sin(eps_rad) * math.sin(math.radians(lon_deg))))


def _ascensional_difference(dec_deg: float, lat_deg: float) -> float:
    """Różnica wschodnia: o ile półłuk dobowy odbiega od 90°.

    sin(AD) = tan(φ)·tan(δ). Gdy |tan(φ)·tan(δ)| ≥ 1, punkt jest okołobiegunowy
    (nigdy nie wschodzi albo nigdy nie zachodzi) i półłuk nie istnieje."""
    v = math.tan(math.radians(lat_deg)) * math.tan(math.radians(dec_deg))
    if abs(v) >= 1.0:
        raise HouseSystemUndefined(
            f"punkt okołobiegunowy (tan φ·tan δ = {v:.6f}): półłuk dobowy nie istnieje")
    return math.degrees(math.asin(v))


# Dom → (ułamek półłuku, czy łuk NOCNY). Domy 11 i 12 dzielą łuk dzienny licząc
# od MC; domy 2 i 3 — łuk nocny, licząc WSTECZ od IC.
_PLACIDUS_STEP = {10: (1 / 3, False), 11: (2 / 3, False),
                  1: (2 / 3, True), 2: (1 / 3, True)}


def _placidus_cusp(ramc: float, eps: float, lat: float,
                   fraction: float, nocturnal: bool) -> float:
    """Punkt ekliptyki, który przebył `fraction` swojego półłuku.

    Warunek jest uwikłany: półłuk zależy od deklinacji, deklinacja od długości,
    a długość od położenia — więc iterujemy po punkcie stałym. Zbieżność jest
    szybka z dala od koła podbiegunowego i psuje się przy nim, dlatego brak
    zbieżności traktujemy jako wyjście poza dziedzinę, a nie jako wynik."""
    eps_rad = math.radians(eps)
    # Start od podziału równomiernego — to Porphyry na równiku, czyli dokładnie
    # ten przypadek, w którym Placidus się do niego sprowadza.
    ra = ramc + 180.0 - 90.0 * fraction if nocturnal else ramc + 90.0 * fraction
    for _ in range(_ITER_MAX):
        dec = _declination_of_ecliptic_lon(_ra_to_ecliptic_lon(ra, eps_rad), eps_rad)
        ad = _ascensional_difference(dec, lat)
        nxt = (ramc + 180.0 - fraction * (90.0 - ad) if nocturnal
               else ramc + fraction * (90.0 + ad))
        if abs(nxt - ra) < _ITER_TOL_DEG:
            return _ra_to_ecliptic_lon(nxt, eps_rad)
        ra = nxt
    raise HouseSystemUndefined(
        f"brak zbieżności po {_ITER_MAX} krokach (φ={lat:.4f}, RAMC={ramc:.4f})")


def _cusps_placidus(ramc: float, eps: float, lat: float,
                    asc: float, mc: float) -> list[float]:
    """Półłuki dobowe i nocne dzielone na trzy — każdy punkt swoim własnym łukiem."""
    if abs(lat) >= polar_circle(eps):
        raise HouseSystemUndefined(
            f"φ={lat:.4f}° poza kołem podbiegunowym ({polar_circle(eps):.4f}° dla ε={eps:.4f}°)")
    inter = {i: _placidus_cusp(ramc, eps, lat, f, noct)
             for i, (f, noct) in _PLACIDUS_STEP.items()}
    inter.update({(i + 6) % 12: norm360(v + 180.0) for i, v in list(inter.items())})
    return _with_exact_angles(inter, asc, mc)


# Koch dzieli CZAS, nie łuk na niebie. Kryterium: ile czasu minęło od wschodu
# tego stopnia zodiaku, który stoi na MC. Ten odcinek (półłuk dobowy stopnia MC)
# dzielimy na trzy i dla punktów podziału liczymy ZWYKŁY Ascendent — stąd nazwa
# „system miejsca urodzenia". Zgodne z definicją Astrodienst (astro.com/astrowiki).
#
# W przeciwieństwie do Placidusa NIE wymaga iteracji: półłuk zależy od deklinacji
# stopnia MC, którą znamy wprost. Granicę dziedziny dzieli natomiast z Placidusem —
# gdy stopień MC jest okołobiegunowy, „moment jego wschodu" nie istnieje.
_KOCH_OFFSET = {10: -2 / 3, 11: -1 / 3, 1: 1 / 3, 2: 2 / 3}


def _cusps_koch(ramc: float, eps: float, lat: float,
                asc: float, mc: float) -> list[float]:
    """Ascendenty dla chwil trójdzielących drogę stopnia MC od wschodu do górowania."""
    if abs(lat) >= polar_circle(eps):
        raise HouseSystemUndefined(
            f"φ={lat:.4f}° poza kołem podbiegunowym ({polar_circle(eps):.4f}° dla ε={eps:.4f}°)")
    dec_mc = _declination_of_ecliptic_lon(mc, math.radians(eps))
    half_arc = 90.0 + _ascensional_difference(dec_mc, lat)
    inter = {i: compute_asc(ramc + f * half_arc, eps, lat)
             for i, f in _KOCH_OFFSET.items()}
    inter.update({(i + 6) % 12: norm360(v + 180.0) for i, v in list(inter.items())})
    return _with_exact_angles(inter, asc, mc)


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 == WHOLE_SIGN_ARIES:
        # Znaki jako domy, ale numeracja rusza od Barana niezależnie od Ascendentu.
        # Wariant spotykany w tradycji indyjskiej i w części szkół hellenistycznych.
        return [norm360(30.0 * i) for i in range(12)]
    if system == EQUAL_MC:
        # Równe domy jak `equal`, ale zakotwiczone na MC: dom X zaczyna się
        # DOKŁADNIE na MC, więc oś południka wypada na granicy domu, a nie w środku.
        return [norm360(mc + 90.0 + 30.0 * i) for i in range(12)]
    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)
    if system == PLACIDUS:
        return _cusps_placidus(ramc, eps, lat, asc, mc)
    if system == KOCH:
        return _cusps_koch(ramc, eps, lat, asc, mc)
    raise ValueError(f"nieznany system domów: {system}")


# ── jawny fallback poza dziedziną ────────────────────────────────────────
# Placidus i Koch jako jedyne mają miejsca, w których po prostu NIE ISTNIEJĄ.
# Astrolog z Tromsø ma dostać wynik, ale musi wiedzieć, że dostał inny system —
# cicha podmiana jest gorsza niż brak wyniku, bo jest nie do wykrycia z wykresu.

FALLBACK_SYSTEM = PORPHYRY


@dataclass(frozen=True)
class CuspSet:
    """12 cuspów + uczciwa informacja, czym naprawdę zostały policzone."""

    cusps: list[float]
    system: str                          # system FAKTYCZNIE użyty
    requested: str                       # o który poproszono
    reason: str | None = None            # dlaczego nie dało się użyć żądanego

    @property
    def is_fallback(self) -> bool:
        return self.system != self.requested

    @property
    def notice(self) -> str | None:
        """Komunikat dla człowieka. Ma trafić na ekran, do raportu i do PDF-a."""
        if not self.is_fallback:
            return None
        return (f"UWAGA: system domów \u201e{self.requested}\u201d nie ma definicji "
                f"dla tego miejsca i czasu \u2014 {self.reason}. Domy policzono "
                f"systemem \u201e{self.system}\u201d. To NIE jest ten sam podzia\u0142: "
                f"pozycje planet s\u0105 poprawne, ale przypisanie ich do dom\u00f3w "
                f"pochodzi z innego systemu.")


def cusps_detailed(ramc: float, eps: float, lat: float, system: str) -> CuspSet:
    """Jak cusps_for, ale zamiast wyjątku poza dziedziną robi JAWNY fallback.

    cusps_for zostaje funkcją czystą i nieustępliwą (to ją porównuje wyrocznia);
    ustępstwo wobec rzeczywistości jest tutaj — i zawsze zostawia ślad."""
    try:
        return CuspSet(cusps_for(ramc, eps, lat, system), system, system)
    except HouseSystemUndefined as e:
        return CuspSet(cusps_for(ramc, eps, lat, FALLBACK_SYSTEM),
                       FALLBACK_SYSTEM, system, str(e))


def _runs_forward(cusp_list: list[float]) -> bool:
    """Czy domy biegną w stronę rosnących długości ekliptycznych.

    Zwykle tak — ale NIE ZAWSZE. Przy dużych szerokościach systemy dzielące koła
    wielkie (regiomontanus, campanus, topocentric) mają kolejność ODWRÓCONĄ:
    przy φ = −84,3° cusp domu I wypada na 174,2°, a domu II na 165,3°. To nie
    jest błąd — wyrocznia zwraca dokładnie te same wartości.

    Rozstrzygamy sumą przeskoków „do przodu": dwanaście cuspów dzieli okrąg, więc
    idąc we WŁAŚCIWĄ stronę zsumują się do 360°. Idąc pod prąd każdy przeskok
    obchodzi koło dookoła i suma wychodzi 11 × 360° = 3960°."""
    total = sum((cusp_list[(i + 1) % 12] - cusp_list[i]) % 360.0 for i in range(12))
    return abs(total - 360.0) < abs(total - 3960.0)


def assign_house(lon: float, cusp_list: list[float]) -> int:
    """Numer domu (1..12), w którym leży dana długość ekliptyczna.

    Kierunek liczenia bierzemy z samych cuspów. Zaszycie „zawsze do przodu"
    dawało przy |φ| powyżej koła podbiegunowego złe domy dla regiomontanusa,
    campanusa i topocentrica — mimo cuspów zgodnych z wyrocznią co do zera.
    Błąd był CICHY: wykres wyglądał poprawnie, tylko planety siedziały gdzie
    indziej. Sprawdzane wobec swe_house_pos (tests/oracle)."""
    lon = norm360(lon)
    forward = _runs_forward(cusp_list)
    for i in range(12):
        start = cusp_list[i]
        end = cusp_list[(i + 1) % 12]
        if forward:
            span, offset = (end - start) % 360.0, (lon - start) % 360.0
        else:
            span, offset = (start - end) % 360.0, (start - lon) % 360.0
        if offset < span:
            return i + 1
    return 12