services/logic/app/engine/houses.py

@@ -14,7 +14,22 @@ from app.engine.formats import SIGN_ABBR, norm360, sign_index  # noqa: F401
 WHOLE_SIGN = "whole_sign"
 EQUAL = "equal"
 PORPHYRY = "porphyry"
-SYSTEMS = (WHOLE_SIGN, EQUAL, 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:
@@ -109,6 +124,213 @@ def polar_circle(eps_deg: float) -> float:
     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.
 
@@ -119,7 +341,22 @@ def cusps_for(ramc: float, eps: float, lat: float, system: str) -> list[float]:
     """
     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:

services/logic/tests/test_houses.py

@@ -75,3 +75,54 @@ def test_polar_circle_moves_with_obliquity():
     a ε zmienia się z datą. Testy brzegowe muszą ją liczyć per data."""
     assert H.polar_circle(23.4393) == pytest.approx(66.5607, abs=1e-4)   # dziś
     assert H.polar_circle(23.747) == pytest.approx(66.253, abs=1e-4)     # 370 p.n.e.
+
+
+# ── systemy egzotyczne o zamkniętym wzorze (Etap 1) ──────────────────────
+# Zgodność z wyrocznią sprawdza tests/oracle; tu pilnujemy niezmienników, które
+# muszą zachodzić także bez swissepha (czyli w każdym środowisku).
+
+EXOTIC = ("vehlow", "morinus", "regiomontanus", "campanus", "alcabitus")
+
+
+@pytest.mark.parametrize("system", EXOTIC)
+def test_exotic_returns_twelve_cusps_in_range(system):
+    out = H.cusps_for(100.0, 23.4393, 50.0, system)
+    assert len(out) == 12
+    assert all(0.0 <= c < 360.0 for c in out)
+
+
+@pytest.mark.parametrize("system", EXOTIC)
+def test_opposite_houses_are_antipodal(system):
+    """Domy przeciwległe leżą na tym samym kole wielkim, więc ich cuspy są
+    dokładnie antypodyczne. Naruszenie tego oznacza błąd w wyborze gałęzi."""
+    out = H.cusps_for(137.0, 23.4393, 42.0, system)
+    for i in range(6):
+        assert abs(((out[i + 6] - out[i]) % 360.0) - 180.0) < 1e-9, f"domy {i+1}/{i+7}"
+
+
+@pytest.mark.parametrize("system", ("regiomontanus", "campanus", "alcabitus"))
+def test_quadrant_systems_anchor_on_ascendant(system):
+    """Systemy kwadrantowe zaczynają dom I na Ascendencie."""
+    ramc, eps, lat = 100.0, 23.4393, 50.0
+    assert H.cusps_for(ramc, eps, lat, system)[0] == pytest.approx(
+        H.compute_asc(ramc, eps, lat), abs=1e-9)
+
+
+def test_vehlow_puts_ascendant_in_the_middle_of_house_one():
+    ramc, eps, lat = 100.0, 23.4393, 50.0
+    asc = H.compute_asc(ramc, eps, lat)
+    assert H.cusps_for(ramc, eps, lat, "vehlow")[0] == pytest.approx(
+        H.norm360(asc - 15.0), abs=1e-9)
+
+
+def test_morinus_ignores_latitude():
+    """Morinus rzutuje równik wprost na ekliptykę, bez horyzontu — jako jedyny
+    nie zależy od szerokości geograficznej."""
+    a = H.cusps_for(100.0, 23.4393, 20.0, "morinus")
+    b = H.cusps_for(100.0, 23.4393, 65.0, "morinus")
+    assert a == pytest.approx(b, abs=1e-12)
+
+
+def test_unknown_system_is_rejected():
+    with pytest.raises(ValueError):
+        H.cusps_for(100.0, 23.4393, 50.0, "nie-ma-takiego")

tests/oracle/README.md

@@ -61,6 +61,34 @@ wykrył **dwa realne błędy** w pierwszym przebiegu:
 Oba mają teraz testy regresji w `services/logic/tests/test_houses.py`, więc są
 łapane także bez swissepha.
 
+## Stan systemów domów
+
+Etap 1 — systemy o **zamkniętym wzorze** (bez iteracji). Każdy poniższy przeszedł
+zarówno zestaw brzegowy (2520 porównań), jak i losowy przemiał 20 000 przypadków
+(240 000 porównań na system):
+
+| System | Konstrukcja | Maks. odchylenie |
+|---|---|---|
+| whole sign | podział ekliptyki | 0,000000000° |
+| equal | podział ekliptyki | 0,000000000° |
+| porphyry | podział kwadrantów po ekliptyce | 0,000000000° |
+| vehlow | equal z Ascendentem w środku domu I | 0,000000000° |
+| morinus | równik rzutowany wprost na ekliptykę | 0,000000000° |
+| regiomontanus | podział równika, koła przez punkty N/S horyzontu | 0,000000000° |
+| campanus | podział wertykału pierwszego | 0,000000000° |
+| alcabitus | podział łuków dobowych po równiku | 0,000000000° |
+
+**Topocentric (Polich–Page) jest zaimplementowany, ale NIE wypuszczony** — nie ma go
+w `houses.SYSTEMS`. Zgadza się z wyrocznią na całej dziedzinie poza otoczeniem
+bieguna: przy |φ| ≈ 89,9° i RAMC 90°/270° kolejność domów się odwraca i żadna reguła
+oparta na łuku kwadrantu nie rozstrzyga wyboru gałęzi. Konstrukcja jest tam z natury
+źle uwarunkowana — „biegun" `atan(tan(φ)·k/3)` dla φ = 89,9° wynosi już 89,7°.
+Nie zawężamy dziedziny po to, żeby testy przeszły; system dołączy, gdy będzie
+poprawny wszędzie.
+
+Placidus i Koch (iteracyjne, z realną granicą dziedziny na kole podbiegunowym) —
+Etap 2.
+
 ## Licencja
 
 Swiss Ephemeris jest na AGPL i jest tu **wyłącznie wyrocznią testową** — nie wchodzi

tests/oracle/run.py

@@ -27,7 +27,8 @@ from harness import (  # noqa: E402
 
 # Systemy do sprawdzenia. Rośnie wraz z implementacją kolejnych (Etap 1 i 2) —
 # dopisanie nazwy tutaj wystarcza, żeby weszła do każdego builda.
-SYSTEMS = ["whole_sign", "equal", "porphyry"]
+SYSTEMS = ["whole_sign", "equal", "porphyry", "vehlow", "morinus",
+           "regiomontanus", "campanus", "alcabitus"]
 
 
 def main() -> int: