geodesic-wallpaper 1.6.1

//! Helicoid surface parameterization with analytic Christoffel symbols.
//!
//! The helicoid is a ruled minimal surface generated by a straight line
//! rotating about and translating along an axis simultaneously.

use super::Surface;
use glam::Vec3;
use std::f32::consts::TAU;

/// Helicoid surface with pitch parameter `c`.
///
/// # Parameterization
/// - `u ∈ [0, 2π)` — rotation angle
/// - `v ∈ [-2, 2]` — radial distance from axis
///
/// ```text
/// x = v · cos(u)
/// y = v · sin(u)
/// z = c · u
/// ```
///
/// Metric: `g_uu = v² + c²`, `g_uv = 0`, `g_vv = 1`.
pub struct Helicoid {
    /// Pitch: axial advance per unit of `u`.
    pub c: f32,
}

impl Helicoid {
    /// Construct a helicoid with the given pitch `c`.
    pub fn new(c: f32) -> Self {
        Self { c }
    }

    /// First partial derivative `∂φ/∂u` at `(u, v)`.
    fn d_du(&self, u: f32, v: f32) -> Vec3 {
        Vec3::new(-v * u.sin(), v * u.cos(), self.c)
    }

    /// First partial derivative `∂φ/∂v` at `(u, v)`.
    fn d_dv(&self, u: f32, _v: f32) -> Vec3 {
        Vec3::new(u.cos(), u.sin(), 0.0)
    }
}

impl Surface for Helicoid {
    fn position(&self, u: f32, v: f32) -> Vec3 {
        Vec3::new(v * u.cos(), v * u.sin(), self.c * u)
    }

    fn metric(&self, u: f32, v: f32) -> [[f32; 2]; 2] {
        let e1 = self.d_du(u, v);
        let e2 = self.d_dv(u, v);
        [[e1.dot(e1), e1.dot(e2)], [e1.dot(e2), e2.dot(e2)]]
    }

    fn christoffel(&self, _u: f32, v: f32) -> [[[f32; 2]; 2]; 2] {
        // g_uu = v² + c², g_vv = 1, g_uv = 0
        // g^uu = 1/(v² + c²), g^vv = 1
        //
        // Non-zero Christoffel symbols:
        //   Γ^u_uv = (1/2) g^uu ∂_v g_uu = (1/2)(1/(v²+c²)) · 2v = v/(v²+c²)
        //   Γ^v_uu = -(1/2) g^vv ∂_v g_uu = -(1/2) · 1 · 2v = -v
        let denom = v * v + self.c * self.c;
        let gamma_u_uv = v / denom.max(1e-8);
        let gamma_v_uu = -v;

        [
            // k=0 (u): Γ^u_ij
            [[0.0, gamma_u_uv], [gamma_u_uv, 0.0]],
            // k=1 (v): Γ^v_ij
            [[gamma_v_uu, 0.0], [0.0, 0.0]],
        ]
    }

    fn wrap(&self, u: f32, v: f32) -> (f32, f32) {
        let u = u.rem_euclid(TAU);
        let v = v.clamp(-2.0, 2.0);
        (u, v)
    }

    fn normal(&self, u: f32, v: f32) -> Vec3 {
        let e1 = self.d_du(u, v);
        let e2 = self.d_dv(u, v);
        e1.cross(e2).normalize()
    }

    fn random_position(&self, rng: &mut dyn rand::RngCore) -> (f32, f32) {
        use rand::Rng;
        (rng.gen_range(0.0..TAU), rng.gen_range(-2.0..2.0_f32))
    }

    fn random_tangent(&self, u: f32, v: f32, rng: &mut dyn rand::RngCore) -> (f32, f32) {
        use rand::Rng;
        let angle: f32 = rng.gen_range(0.0..TAU);
        let g = self.metric(u, v);
        // Normalise so g_ij du^i du^j = 1.
        let du = angle.cos() / g[0][0].sqrt().max(1e-8);
        let dv = angle.sin() / g[1][1].sqrt().max(1e-8);
        (du, dv)
    }

    fn mesh_vertices(&self, u_steps: u32, v_steps: u32) -> (Vec<[f32; 3]>, Vec<u32>) {
        let mut verts = Vec::new();
        let mut indices = Vec::new();
        for i in 0..=u_steps {
            for j in 0..=v_steps {
                let u = (i as f32 / u_steps as f32) * TAU;
                let v = -2.0 + (j as f32 / v_steps as f32) * 4.0;
                let p = self.position(u, v);
                verts.push([p.x, p.y, p.z]);
            }
        }
        for i in 0..u_steps {
            for j in 0..v_steps {
                let a = i * (v_steps + 1) + j;
                let b = a + 1;
                let c = (i + 1) * (v_steps + 1) + j;
                let d = c + 1;
                indices.extend_from_slice(&[a, b, c, b, d, c]);
            }
        }
        (verts, indices)
    }
}

// ─── Tests ────────────────────────────────────────────────────────────────────

#[cfg(test)]
#[allow(clippy::needless_range_loop)]
mod tests {
    use super::*;

    #[test]
    fn position_at_known_point() {
        // At u=0, v=1: x = 1·cos(0) = 1, y = 0, z = 0.
        let hel = Helicoid::new(1.0);
        let p = hel.position(0.0, 1.0);
        assert!((p.x - 1.0).abs() < 1e-5, "x={}", p.x);
        assert!(p.y.abs() < 1e-5, "y={}", p.y);
        assert!(p.z.abs() < 1e-5, "z={}", p.z);
    }

    #[test]
    fn metric_is_diagonal() {
        let hel = Helicoid::new(1.0);
        let g = hel.metric(0.5, 1.0);
        assert!((g[0][1]).abs() < 1e-5, "g_01 not zero: {}", g[0][1]);
    }

    #[test]
    fn metric_g_uu_formula() {
        // g_uu = v² + c²
        let hel = Helicoid::new(2.0);
        let v = 1.5_f32;
        let g = hel.metric(0.0, v);
        let expected = v * v + 4.0;
        assert!(
            (g[0][0] - expected).abs() < 1e-4,
            "g_uu={} expected={}",
            g[0][0],
            expected
        );
    }

    #[test]
    fn christoffel_symmetry_ij() {
        let hel = Helicoid::new(1.0);
        let g = hel.christoffel(0.4, 0.8);
        for k in 0..2 {
            assert!(
                (g[k][0][1] - g[k][1][0]).abs() < 1e-6,
                "Γ^{k}_01 != Γ^{k}_10"
            );
        }
    }

    #[test]
    fn wrap_clamps_v_and_wraps_u() {
        let hel = Helicoid::new(1.0);
        let (u, v) = hel.wrap(TAU + 0.3, 5.0);
        assert!((0.0..TAU).contains(&u), "u={u}");
        assert!(v <= 2.0, "v not clamped: {v}");
        assert!(v >= -2.0, "v not clamped: {v}");
    }

    #[test]
    fn normal_is_unit() {
        let hel = Helicoid::new(1.0);
        let n = hel.normal(1.0, 0.5);
        assert!(
            (n.length() - 1.0).abs() < 1e-5,
            "normal length={}",
            n.length()
        );
    }
}