feat: Koch, Regiomontanus, Campanus, Topocentric, Alcabitius, Morinus…

README.md

@@ -19,7 +19,8 @@ on the server (Python bindings), in the desktop and mobile apps (native) and in
 
 > Status: the complete chart pipeline is ported and matches the reference implementation
 > on every golden chart: planets, lunar nodes, Chiron, Lilith, houses (Placidus, Whole Sign,
-> Equal, Porphyry), aspects and orbs, fixed stars, Arabic parts, temperament, lunar status,
+> Equal, Porphyry; Koch, Regiomontanus, Campanus, Topocentric, Alcabitius, Morinus and Vehlow
+> are checked against Swiss Ephemeris), aspects and orbs, fixed stars, Arabic parts, temperament, lunar status,
 > horary charts (planetary hours from computed sunrise and sunset), transits, synastry and
 > derived charts, with Python and WebAssembly bindings.
 

crates/astroceleste-engine/src/houses.rs

@@ -1,5 +1,9 @@
-//! House cusps, Ascendant and Midheaven (`SkyfieldEngine.calculate_houses`):
-//! Placidus (default), Whole Sign, Equal and Porphyry.
+//! House cusps, Ascendant and Midheaven: Placidus (default), Koch, Regiomontanus, Campanus,
+//! Topocentric (Polich-Page), Alcabitius, Morinus, Porphyry, Equal, Vehlow and Whole Sign.
+//!
+//! Placidus, Whole Sign, Equal and Porphyry follow the reference implementation
+//! (`SkyfieldEngine.calculate_houses`); the other systems use the standard formulas, as in
+//! Swiss Ephemeris' `swehouse.c`, against which the application tests them.
 
 use std::f64::consts::PI;
 
@@ -11,18 +15,33 @@ use crate::time::Time;
 #[non_exhaustive]
 pub enum HouseSystem {
     Placidus,
-    WholeSign,
-    Equal,
+    Koch,
+    Regiomontanus,
+    Campanus,
+    Topocentric,
+    Alcabitius,
+    Morinus,
     Porphyry,
+    Equal,
+    Vehlow,
+    WholeSign,
 }
 
 impl HouseSystem {
-    /// From the first letter of the code, like the reference: anything unknown is Placidus.
+    /// From the first letter of the code (the Swiss Ephemeris letters), like the reference:
+    /// anything unknown is Placidus.
     pub fn from_code(code: &str) -> Self {
         match code.chars().next().map(|c| c.to_ascii_uppercase()) {
-            Some('W') => HouseSystem::WholeSign,
-            Some('E') => HouseSystem::Equal,
+            Some('K') => HouseSystem::Koch,
+            Some('R') => HouseSystem::Regiomontanus,
+            Some('C') => HouseSystem::Campanus,
+            Some('T') => HouseSystem::Topocentric,
+            Some('B') => HouseSystem::Alcabitius,
+            Some('M') => HouseSystem::Morinus,
             Some('O') => HouseSystem::Porphyry,
+            Some('E') => HouseSystem::Equal,
+            Some('V') => HouseSystem::Vehlow,
+            Some('W') => HouseSystem::WholeSign,
             _ => HouseSystem::Placidus,
         }
     }
@@ -71,20 +90,96 @@ pub fn calculate_houses(jd: f64, lat: f64, lon: f64, system: HouseSystem, shift:
                 *cusp = pyfloat::rem(asc_deg + (i as f64 - 1.0) * 30.0, 360.0);
             }
         }
-        HouseSystem::Porphyry => {
+        HouseSystem::Vehlow => {
+            // Equal houses with the Ascendant in the middle of the first house.
+            for (i, cusp) in c.iter_mut().enumerate().skip(1) {
+                *cusp = pyfloat::rem(asc_deg - 15.0 + (i as f64 - 1.0) * 30.0 + 360.0, 360.0);
+            }
+        }
+        HouseSystem::Morinus => {
+            // The equator divided into twelve from the RAMC, projected onto the ecliptic
+            // along ecliptic meridians: the angles are not cusps.
+            for i in 1..=12 {
+                let ra = (ramc + 30.0 * i as f64).to_radians();
+                c[(i + 9) % 12 + 1] = norm((ra.sin() * eps_r.cos()).atan2(ra.cos()).to_degrees());
+            }
+        }
+        HouseSystem::Koch if lat.abs() >= 90.0 - eps_deg => {
+            // Koch is undefined inside the polar circles: Porphyry, as Swiss Ephemeris does.
+            porphyry(&mut c, asc_deg, mc_deg);
+        }
+        HouseSystem::Koch
+        | HouseSystem::Regiomontanus
+        | HouseSystem::Campanus
+        | HouseSystem::Topocentric
+        | HouseSystem::Alcabitius => {
             set_angles(&mut c, asc_deg, mc_deg);
-            let q1 = pyfloat::rem(asc_deg - mc_deg, 360.0);
-            c[11] = pyfloat::rem(mc_deg + q1 / 3.0, 360.0);
-            c[12] = pyfloat::rem(mc_deg + (2.0 * q1) / 3.0, 360.0);
-            let q2 = pyfloat::rem(c[4] - asc_deg, 360.0);
-            c[2] = pyfloat::rem(asc_deg + q2 / 3.0, 360.0);
-            c[3] = pyfloat::rem(asc_deg + (2.0 * q2) / 3.0, 360.0);
-            // Houses 5, 6 oppose 11, 12; houses 8, 9 oppose 2, 3.
-            c[5] = pyfloat::rem(c[11] + 180.0, 360.0);
-            c[6] = pyfloat::rem(c[12] + 180.0, 360.0);
-            c[8] = pyfloat::rem(c[2] + 180.0, 360.0);
-            c[9] = pyfloat::rem(c[3] + 180.0, 360.0);
+            // (right ascension of the house circle's east point, its pole height) for
+            // cusps 11, 12, 2 and 3; the others are opposite.
+            let circles: [(f64, f64); 4] = match system {
+                HouseSystem::Koch => {
+                    let sin_a =
+                        (mc_deg.to_radians().sin() * eps_r.sin() / lat_r.cos()).clamp(-1.0, 1.0);
+                    let cos_a = (1.0 - sin_a * sin_a).sqrt();
+                    let pole = (lat_r.tan() / cos_a).atan();
+                    let ad3 = (pole.sin() * sin_a).asin().to_degrees() / 3.0;
+                    [
+                        (ramc + 30.0 - 2.0 * ad3, lat),
+                        (ramc + 60.0 - ad3, lat),
+                        (ramc + 120.0 + ad3, lat),
+                        (ramc + 150.0 + 2.0 * ad3, lat),
+                    ]
+                }
+                HouseSystem::Regiomontanus | HouseSystem::Topocentric => {
+                    let (k1, k2) = if system == HouseSystem::Regiomontanus {
+                        (0.5, 30.0_f64.to_radians().cos())
+                    } else {
+                        (1.0 / 3.0, 2.0 / 3.0)
+                    };
+                    let p1 = (lat_r.tan() * k1).atan().to_degrees();
+                    let p2 = (lat_r.tan() * k2).atan().to_degrees();
+                    [
+                        (ramc + 30.0, p1),
+                        (ramc + 60.0, p2),
+                        (ramc + 120.0, p2),
+                        (ramc + 150.0, p1),
+                    ]
+                }
+                HouseSystem::Campanus => {
+                    let p1 = (lat_r.sin() / 2.0).asin().to_degrees();
+                    let p2 = (3.0_f64.sqrt() / 2.0 * lat_r.sin()).asin().to_degrees();
+                    let x1 = (3.0_f64.sqrt() / lat_r.cos()).atan().to_degrees();
+                    let x2 = (1.0 / 3.0_f64.sqrt() / lat_r.cos()).atan().to_degrees();
+                    [
+                        (ramc + 90.0 - x1, p1),
+                        (ramc + 90.0 - x2, p2),
+                        (ramc + 90.0 + x2, p2),
+                        (ramc + 90.0 + x1, p1),
+                    ]
+                }
+                _ => {
+                    // Alcabitius: the Ascendant's diurnal and nocturnal semi-arcs trisected
+                    // on the equator, projected along hour circles.
+                    let dec = (asc_deg.to_radians().sin() * eps_r.sin()).asin();
+                    let sda = (-lat_r.tan() * dec.tan())
+                        .clamp(-1.0, 1.0)
+                        .acos()
+                        .to_degrees();
+                    let sna = 180.0 - sda;
+                    [
+                        (ramc + sda / 3.0, 0.0),
+                        (ramc + 2.0 * sda / 3.0, 0.0),
+                        (ramc + 180.0 - 2.0 * sna / 3.0, 0.0),
+                        (ramc + 180.0 - sna / 3.0, 0.0),
+                    ]
+                }
+            };
+            for (house, (ra, pole)) in [11, 12, 2, 3].into_iter().zip(circles) {
+                c[house] = ecliptic_intersection(ra, pole, eps_r);
+            }
+            opposite_cusps(&mut c);
         }
+        HouseSystem::Porphyry => porphyry(&mut c, asc_deg, mc_deg),
         HouseSystem::Placidus => {
             set_angles(&mut c, asc_deg, mc_deg);
             let solve = |offset: f64, semi_factor: f64| {
@@ -131,9 +226,148 @@ pub fn calculate_houses(jd: f64, lat: f64, lon: f64, system: HouseSystem, shift:
     }
 }
 
+fn norm(deg: f64) -> f64 {
+    pyfloat::rem(pyfloat::rem(deg, 360.0) + 360.0, 360.0)
+}
+
+/// Longitude where the ecliptic meets the house circle whose east point has right ascension
+/// `ra` and whose pole height is `pole` (degrees). With `ra` = RAMC + 90° and `pole` = the
+/// latitude it is the Ascendant; with `pole` = 0 it projects `ra` along its hour circle.
+fn ecliptic_intersection(ra: f64, pole: f64, eps_r: f64) -> f64 {
+    let ra = ra.to_radians();
+    let x = ra.cos() * eps_r.cos() - pole.to_radians().tan() * eps_r.sin();
+    norm(ra.sin().atan2(x).to_degrees())
+}
+
+/// Porphyry: each quadrant between the angles trisected in longitude.
+fn porphyry(c: &mut [f64; 13], asc_deg: f64, mc_deg: f64) {
+    set_angles(c, asc_deg, mc_deg);
+    let q1 = pyfloat::rem(asc_deg - mc_deg, 360.0);
+    c[11] = pyfloat::rem(mc_deg + q1 / 3.0, 360.0);
+    c[12] = pyfloat::rem(mc_deg + (2.0 * q1) / 3.0, 360.0);
+    let q2 = pyfloat::rem(c[4] - asc_deg, 360.0);
+    c[2] = pyfloat::rem(asc_deg + q2 / 3.0, 360.0);
+    c[3] = pyfloat::rem(asc_deg + (2.0 * q2) / 3.0, 360.0);
+    opposite_cusps(c);
+}
+
+/// Houses 5, 6 oppose 11, 12; houses 8, 9 oppose 2, 3.
+fn opposite_cusps(c: &mut [f64; 13]) {
+    c[5] = pyfloat::rem(c[11] + 180.0, 360.0);
+    c[6] = pyfloat::rem(c[12] + 180.0, 360.0);
+    c[8] = pyfloat::rem(c[2] + 180.0, 360.0);
+    c[9] = pyfloat::rem(c[3] + 180.0, 360.0);
+}
+
 fn set_angles(c: &mut [f64; 13], asc: f64, mc: f64) {
     c[1] = asc;
     c[10] = mc;
     c[7] = pyfloat::rem(asc + 180.0, 360.0);
     c[4] = pyfloat::rem(mc + 180.0, 360.0);
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    const JD: f64 = 2_448_027.104_166_666_5; // 1990-05-15 14:30 UT, Rome (41.9° N, 12.5° E)
+
+    fn gap(a: f64, b: f64) -> f64 {
+        let d = (a - b).rem_euclid(360.0);
+        d.min(360.0 - d)
+    }
+
+    /// Cusps within 1′ of Swiss Ephemeris 2.10 (`swe.houses`) for the same moment.
+    #[test]
+    fn systems_match_swiss_ephemeris() {
+        let reference: [(&str, [f64; 12]); 7] = [
+            (
+                "K",
+                [
+                    190.2408, 219.3528, 248.6584, 282.0137, 311.8995, 340.9197, 10.2408, 39.3528,
+                    68.6584, 102.0137, 131.8995, 160.9197,
+                ],
+            ),
+            (
+                "R",
+                [
+                    190.2408, 214.8788, 245.0128, 282.0137, 317.7672, 346.2050, 10.2408, 34.8788,
+                    65.0128, 102.0137, 137.7672, 166.2050,
+                ],
+            ),
+            (
+                "C",
+                [
+                    190.2408, 222.0070, 252.8709, 282.0137, 310.3369, 339.4310, 10.2408, 42.0070,
+                    72.8709, 102.0137, 130.3369, 159.4310,
+                ],
+            ),
+            (
+                "T",
+                [
+                    190.2408, 216.9311, 248.0089, 282.0137, 315.5716, 345.3618, 10.2408, 36.9311,
+                    68.0089, 102.0137, 135.5716, 165.3618,
+                ],
+            ),
+            (
+                "B",
+                [
+                    190.2408, 223.0794, 253.2543, 282.0137, 309.4047, 339.0304, 10.2408, 43.0794,
+                    73.2543, 102.0137, 129.4047, 159.0304,
+                ],
+            ),
+            (
+                "M",
+                [
+                    192.0137, 220.6069, 251.6333, 284.1881, 315.5252, 344.3861, 12.0137, 40.6069,
+                    71.6333, 104.1881, 135.5252, 164.3861,
+                ],
+            ),
+            (
+                "V",
+                [
+                    175.2408, 205.2408, 235.2408, 265.2408, 295.2408, 325.2408, 355.2408, 25.2408,
+                    55.2408, 85.2408, 115.2408, 145.2408,
+                ],
+            ),
+        ];
+        for (code, expected) in reference {
+            let houses = calculate_houses(JD, 41.9, 12.5, HouseSystem::from_code(code), 0.0);
+            for (i, (ours, theirs)) in houses.cusps.iter().zip(expected).enumerate() {
+                assert!(
+                    gap(*ours, theirs) < 1.0 / 60.0,
+                    "{code} cusp {}: {ours:.4} vs {theirs:.4}",
+                    i + 1
+                );
+            }
+        }
+    }
+
+    #[test]
+    fn codes_select_their_systems() {
+        for (code, system) in [
+            ("P", HouseSystem::Placidus),
+            ("koch", HouseSystem::Koch),
+            ("R", HouseSystem::Regiomontanus),
+            ("C", HouseSystem::Campanus),
+            ("T", HouseSystem::Topocentric),
+            ("B", HouseSystem::Alcabitius),
+            ("M", HouseSystem::Morinus),
+            ("O", HouseSystem::Porphyry),
+            ("E", HouseSystem::Equal),
+            ("V", HouseSystem::Vehlow),
+            ("W", HouseSystem::WholeSign),
+            ("?", HouseSystem::Placidus),
+            ("", HouseSystem::Placidus),
+        ] {
+            assert_eq!(HouseSystem::from_code(code), system, "{code}");
+        }
+    }
+
+    #[test]
+    fn koch_falls_back_to_porphyry_inside_the_polar_circles() {
+        let koch = calculate_houses(JD, 70.0, 25.0, HouseSystem::Koch, 0.0);
+        let porphyry = calculate_houses(JD, 70.0, 25.0, HouseSystem::Porphyry, 0.0);
+        assert_eq!(koch, porphyry);
+    }
+}

docs/accuracy.md

@@ -16,6 +16,10 @@ every key. Three independent references check that.
 | `transits.json`, `synastry.json`, `derived.json` | 27 each |
 | `errors.json` | dates before and after the kernel (1700, 2200) |
 
+The reference implementation only had Placidus, Whole Sign, Equal and Porphyry. The other
+house systems are unit-tested against Swiss Ephemeris (`src/houses.rs`, within 1′), and the
+application cross-checks them at several latitudes.
+
 The engine must reproduce them with:
 
 - floats within 1e-9 (relative to magnitude, absolute below 1);

docs/api.md

@@ -22,14 +22,15 @@ derived charts work on already computed charts and do not.
 | moment | UTC instant: `UtcInstant` in Rust, `datetime` or ISO string in Python, ISO string `utc` in JS | required |
 | latitude | geographic latitude in degrees, north positive | required |
 | longitude | geographic longitude in degrees, east positive | required |
-| house system | `P` Placidus, `W` Whole Sign, `E` Equal, `O` Porphyry | `P` |
+| house system | `P` Placidus, `K` Koch, `R` Regiomontanus, `C` Campanus, `T` Topocentric (Polich-Page), `B` Alcabitius, `M` Morinus, `O` Porphyry, `E` Equal, `V` Vehlow, `W` Whole Sign (the Swiss Ephemeris letters) | `P` |
 | zodiac type | `tropical` or `sidereal` | `tropical` |
 | ayanamsa | sidereal reference, see below (ignored for tropical charts) | `galcent_0sag` |
 | orb settings | overrides merged over the default orbs, see below | none |
 
 Times are always UTC. Convert civil time (with its time zone and daylight saving time) to
 UTC before calling the engine. Only the first letter of the house system is used, and
-unknown codes fall back to Placidus. Zodiac types other than `sidereal` are tropical.
+unknown codes fall back to Placidus. Koch falls back to Porphyry inside the polar circles,
+where it is undefined. Morinus cusps 1 and 10 are not the Ascendant and Midheaven. Zodiac types other than `sidereal` are tropical.
 
 ```rust
 use astroceleste_engine::ephemeris::{Kernel, KernelSet, Spk};

site/index.html

@@ -89,7 +89,14 @@ <h2>Live demo</h2>
                   <option value="P">Placidus</option>
                   <option value="W">Whole Sign</option>
                   <option value="E">Equal</option>
+                  <option value="K">Koch</option>
+                  <option value="R">Regiomontanus</option>
+                  <option value="C">Campanus</option>
+                  <option value="T">Topocentric</option>
+                  <option value="B">Alcabitius</option>
+                  <option value="M">Morinus</option>
                   <option value="O">Porphyry</option>
+                  <option value="V">Vehlow</option>
                 </select>
               </label>
               <label>Zodiac
@@ -208,7 +215,7 @@ <h3>Planets and points</h3>
           </article>
           <article>
             <h3>Houses</h3>
-            <p>Placidus, Whole Sign, Equal and Porphyry cusps, and the house of every point.</p>
+            <p>Eleven house systems (Placidus, Koch, Regiomontanus, Campanus, Topocentric, Alcabitius, Morinus, Porphyry, Equal, Vehlow, Whole Sign) and the house of every point.</p>
           </article>
           <article>
             <h3>Tropical and sidereal</h3>