feat: add kerr_isco and kerr_horizon CPU helpers with tests

src/physics.rs

@@ -52,6 +52,26 @@ pub fn impact_parameter(eye: Vec3, dir: Vec3) -> f32 {
     eye.cross(dir).length()
 }
 
+/// Prograde Kerr ISCO in Rs units (Rs=1, so M=0.5). `chi = a/M ∈ [0,1]`.
+/// Bardeen-Press-Teukolsky (1972) closed form. Returns 6M=3.0 at chi=0,
+/// M=0.5 at chi=1.
+pub fn kerr_isco(chi: f32) -> f32 {
+    let m = 0.5;
+    let cbrt_pos = (1.0 + chi).cbrt();
+    let cbrt_neg = (1.0 - chi).cbrt();
+    let z1 = 1.0 + (1.0 - chi * chi).cbrt() * (cbrt_pos + cbrt_neg);
+    let z2 = (3.0 * chi * chi + z1 * z1).sqrt();
+    m * (3.0 + z2 - ((3.0 - z1) * (3.0 + z1 + 2.0 * z2)).sqrt())
+}
+
+/// Kerr event-horizon radius r+ in Rs units (Rs=1, M=0.5). `chi = a/M ∈ [0,1]`.
+/// Returns Rs=1.0 at chi=0, M=0.5 at chi=1.
+pub fn kerr_horizon(chi: f32) -> f32 {
+    let m = 0.5;
+    let a = chi * m;
+    m + (m * m - a * a).max(0.0).sqrt()
+}
+
 // Silence unused-import warning for Vec4 if not used; kept for future expansion.
 #[allow(dead_code)]
 fn _phantom(_v: Vec4) {}

tests/physics_test.rs

@@ -6,3 +6,49 @@ fn public_bcrt_constant_is_correct() {
     let expected = 1.5 * 3.0_f32.sqrt();
     assert!((physics::BCRIT - expected).abs() < 1e-5);
 }
+
+#[test]
+fn kerr_isco_at_zero_is_schwarzschild() {
+    // spin=0 → ISCO = 6M = 3 Rs (Rs=1).
+    let isco = physics::kerr_isco(0.0);
+    assert!((isco - 3.0).abs() < 1e-3, "spin=0 ISCO should be 3.0, got {}", isco);
+}
+
+#[test]
+fn kerr_isco_at_extremal_is_half_rs() {
+    // spin=1 → ISCO = M = Rs/2 = 0.5.
+    let isco = physics::kerr_isco(1.0);
+    assert!((isco - 0.5).abs() < 1e-3, "spin=1 ISCO should be 0.5, got {}", isco);
+}
+
+#[test]
+fn kerr_isco_is_monotonically_decreasing() {
+    let a = physics::kerr_isco(0.3);
+    let b = physics::kerr_isco(0.6);
+    let c = physics::kerr_isco(0.9);
+    assert!(a > b, "0.3 > 0.6: {} vs {}", a, b);
+    assert!(b > c, "0.6 > 0.9: {} vs {}", b, c);
+}
+
+#[test]
+fn kerr_horizon_at_zero_is_rs() {
+    // spin=0 → r+ = Rs = 1.0.
+    let r = physics::kerr_horizon(0.0);
+    assert!((r - 1.0).abs() < 1e-3, "spin=0 horizon should be 1.0, got {}", r);
+}
+
+#[test]
+fn kerr_horizon_at_extremal_is_half_rs() {
+    // spin=1 → r+ = M = 0.5.
+    let r = physics::kerr_horizon(1.0);
+    assert!((r - 0.5).abs() < 1e-3, "spin=1 horizon should be 0.5, got {}", r);
+}
+
+#[test]
+fn kerr_horizon_is_monotonically_decreasing() {
+    let a = physics::kerr_horizon(0.3);
+    let b = physics::kerr_horizon(0.6);
+    let c = physics::kerr_horizon(0.9);
+    assert!(a > b, "0.3 > 0.6: {} vs {}", a, b);
+    assert!(b > c, "0.6 > 0.9: {} vs {}", b, c);
+}