From 62be271e7eccffd8738aa022f8e409ca12c5dd4b Mon Sep 17 00:00:00 2001 From: ffalcinelli Date: Wed, 23 Sep 2026 08:42:38 +0200 Subject: [PATCH] feat: Koch, Regiomontanus, Campanus, Topocentric, Alcabitius, Morinus and Vehlow houses MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Every code other than P/W/E/O used to fall back to Placidus, so an application offering these systems silently served Placidus cusps under their names. They follow the standard formulas (as in Swiss Ephemeris' swehouse.c) and are unit-tested against Swiss Ephemeris within 1′. Koch falls back to Porphyry inside the polar circles, where it is undefined. Co-Authored-By: Claude Opus 5.5 --- README.md | 3 +- crates/astroceleste-engine/src/houses.rs | 272 +++++++++++++++++++++-- docs/accuracy.md | 4 + docs/api.md | 5 +- site/index.html | 9 +- 5 files changed, 270 insertions(+), 23 deletions(-) diff --git a/README.md b/README.md index 10db9a6..a1202f5 100644 --- a/README.md +++ b/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. diff --git a/crates/astroceleste-engine/src/houses.rs b/crates/astroceleste-engine/src/houses.rs index ec8065e..7a84ec1 100644 --- a/crates/astroceleste-engine/src/houses.rs +++ b/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); + } +} diff --git a/docs/accuracy.md b/docs/accuracy.md index f137623..f379754 100644 --- a/docs/accuracy.md +++ b/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); diff --git a/docs/api.md b/docs/api.md index e578715..8d2ef2a 100644 --- a/docs/api.md +++ b/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}; diff --git a/site/index.html b/site/index.html index ebbea46..ebe060f 100644 --- a/site/index.html +++ b/site/index.html @@ -89,7 +89,14 @@

Live demo

+ + + + + + +