test: CPU-mirrored geodesic physics + bcrit/capture unit tests

Cargo.toml

@@ -3,6 +3,10 @@ name = "singularity-rs"
 version = "0.1.0"
 edition = "2024"
 
+[lib]
+name = "singularity_rs"
+path = "src/lib.rs"
+
 [dependencies]
 bevy = "0.19"
 bevy_egui = "0.41"

src/lib.rs

@@ -0,0 +1 @@
+pub mod physics;

src/physics.rs

@@ -0,0 +1,86 @@
+//! CPU mirror of the shader physics, for unit-testing.
+//! Natural units: Rs = 1.
+
+use bevy::math::{Vec3, Vec4};
+
+pub const RS: f32 = 1.0;
+/// Critical impact parameter for a Schwarzschild hole: bcrit = (3*sqrt(3)/2) * Rs.
+/// Written as a literal because `f32::sqrt` is not `const` on stable Rust.
+/// (3 * sqrt(3) / 2 = 2.598076211353316...)
+pub const BCRIT: f32 = 2.5980762;
+
+/// One Euler step of the discretized geodesic (mirrors the shader's `deriv`).
+pub fn bending_accel(pos: Vec3, dir: Vec3) -> Vec3 {
+    let r = pos.length();
+    let h = pos.cross(dir);
+    let h2 = h.dot(h);
+    let r5 = (r * r * r * r * r).max(1e-6);
+    -1.5 * RS * h2 / r5 * pos
+}
+
+/// Classify a ray by integrating it. Returns true if captured (r < Rs).
+/// Uses RK4 like the shader. `dt` is the step size, `steps` the count.
+pub fn is_captured(start_pos: Vec3, start_dir: Vec3, steps: u32, dt: f32) -> bool {
+    let mut pos = start_pos;
+    let mut dir = start_dir;
+    const R_ESCAPE: f32 = 1000.0;
+    for _ in 0..steps {
+        let r = pos.length();
+        if r < RS {
+            return true;
+        }
+        if r > R_ESCAPE {
+            return false;
+        }
+        let k1p = dir;
+        let k1d = bending_accel(pos, dir);
+        let k2p = dir + k1d * dt * 0.5;
+        let k2d = bending_accel(pos + k1p * dt * 0.5, (dir + k1d * dt * 0.5).normalize());
+        let k3p = dir + k2d * dt * 0.5;
+        let k3d = bending_accel(pos + k2p * dt * 0.5, (dir + k2d * dt * 0.5).normalize());
+        let k4p = dir + k3d * dt;
+        let k4d = bending_accel(pos + k3p * dt, (dir + k3d * dt).normalize());
+        pos += (k1p + 2.0 * k2p + 2.0 * k3p + k4p) * dt / 6.0;
+        dir = (dir + (k1d + 2.0 * k2d + 2.0 * k3d + k4d) * dt / 6.0).normalize();
+    }
+    false
+}
+
+/// Compute the impact parameter of a ray given eye position and direction:
+/// b = |eye x dir| / |dir|  (since dir is unit, b = |eye x dir|).
+pub fn impact_parameter(eye: Vec3, dir: Vec3) -> f32 {
+    eye.cross(dir).length()
+}
+
+// Silence unused-import warning for Vec4 if not used; kept for future expansion.
+#[allow(dead_code)]
+fn _phantom(_v: Vec4) {}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn bcrit_value() {
+        assert!((BCRIT - 2.598).abs() < 0.01, "bcrit should be ~2.598, got {}", BCRIT);
+    }
+
+    #[test]
+    fn ray_below_bcrit_is_captured() {
+        // Eye far on the z-axis; aim slightly off-center with b < bcrit.
+        let eye = Vec3::new(0.0, 0.0, 50.0);
+        let dir = Vec3::new(0.0, 2.0, -50.0).normalize(); // b ~ 2.0 < 2.598
+        let b = impact_parameter(eye, dir);
+        assert!(b < BCRIT, "b {} should be < bcrit {}", b, BCRIT);
+        assert!(is_captured(eye, dir, 2000, 0.1), "ray below bcrit should be captured");
+    }
+
+    #[test]
+    fn ray_above_bcrit_escapes() {
+        let eye = Vec3::new(0.0, 0.0, 50.0);
+        let dir = Vec3::new(0.0, 10.0, -50.0).normalize(); // b ~ 9.8 >> bcrit
+        let b = impact_parameter(eye, dir);
+        assert!(b > BCRIT);
+        assert!(!is_captured(eye, dir, 2000, 0.1), "ray above bcrit should escape");
+    }
+}

tests/physics_test.rs

@@ -0,0 +1,8 @@
+use singularity_rs::physics;
+
+#[test]
+fn public_bcrt_constant_is_correct() {
+    // 3*sqrt(3)/2 ≈ 2.598076
+    let expected = 1.5 * 3.0_f32.sqrt();
+    assert!((physics::BCRIT - expected).abs() < 1e-5);
+}