geodesic-wallpaper 1.6.1

//! Catenoid surface parameterization with analytic Christoffel symbols.
//!
//! The catenoid is a minimal surface of revolution generated by rotating a
//! catenary. It is conformally equivalent to the plane.

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

/// Catenoid surface with scale parameter `c`.
///
/// # Parameterization
/// - `u ∈ [0, 2π)` — azimuthal angle
/// - `v ∈ [-3c, 3c]` — axial (height) parameter
///
/// ```text
/// x = c · cosh(v/c) · cos(u)
/// y = c · cosh(v/c) · sin(u)
/// z = v
/// ```
///
/// The metric is conformal: `g = λ²(du² + dv²)` where `λ = c·cosh(v/c)`.
pub struct Catenoid {
    /// Scale factor controlling the waist radius and domain extent.
    pub c: f32,
}

impl Catenoid {
    /// Construct a catenoid with the given scale `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 {
        let lambda = self.c * (v / self.c).cosh();
        Vec3::new(-lambda * u.sin(), lambda * u.cos(), 0.0)
    }

    /// First partial derivative `∂φ/∂v` at `(u, v)`.
    fn d_dv(&self, u: f32, v: f32) -> Vec3 {
        // ∂/∂v [c·cosh(v/c)] = sinh(v/c)
        let sh = (v / self.c).sinh();
        Vec3::new(sh * u.cos(), sh * u.sin(), 1.0)
    }
}

impl Surface for Catenoid {
    fn position(&self, u: f32, v: f32) -> Vec3 {
        let lambda = self.c * (v / self.c).cosh();
        Vec3::new(lambda * u.cos(), lambda * u.sin(), v)
    }

    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] {
        // Conformal metric g = λ²(du² + dv²) with λ = c·cosh(v/c).
        // σ = ln λ = ln(c·cosh(v/c))
        // ∂_u σ = 0  (λ independent of u)
        // ∂_v σ = tanh(v/c) / c
        //
        // Christoffel symbols for conformal metric:
        //   Γ^u_uu =  ∂_u σ = 0
        //   Γ^u_uv =  ∂_v σ = tanh(v/c)/c
        //   Γ^u_vv = -∂_u σ = 0
        //   Γ^v_uu = -∂_v σ = -tanh(v/c)/c
        //   Γ^v_uv =  ∂_u σ = 0
        //   Γ^v_vv =  ∂_v σ = tanh(v/c)/c
        let sigma_v = (v / self.c).tanh() / self.c;

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

    fn wrap(&self, u: f32, v: f32) -> (f32, f32) {
        let u = u.rem_euclid(TAU);
        let v_max = 3.0 * self.c;
        let v = v.clamp(-v_max, v_max);
        (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;
        let v_max = 3.0 * self.c;
        (rng.gen_range(0.0..TAU), rng.gen_range(-v_max..v_max))
    }

    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 v_max = 3.0 * self.c;
        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 = -v_max + (j as f32 / v_steps as f32) * 2.0 * v_max;
                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=0: x = c·cosh(0)·cos(0) = c, y = 0, z = 0.
        let cat = Catenoid::new(1.0);
        let p = cat.position(0.0, 0.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_conformal() {
        // For conformal metric g_00 = g_11 and g_01 = 0.
        let cat = Catenoid::new(1.0);
        let g = cat.metric(0.5, 0.3);
        assert!((g[0][1]).abs() < 1e-5, "g_01 not zero: {}", g[0][1]);
        assert!(
            (g[0][0] - g[1][1]).abs() < 1e-4,
            "g_00 != g_11: {} {}",
            g[0][0],
            g[1][1]
        );
    }

    #[test]
    fn christoffel_symmetry_ij() {
        let cat = Catenoid::new(1.5);
        let g = cat.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 cat = Catenoid::new(1.0);
        let (u, v) = cat.wrap(TAU + 0.5, 10.0);
        assert!((0.0..TAU).contains(&u), "u={u}");
        assert!(v <= 3.0, "v not clamped: {v}");
    }

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

    #[test]
    fn catenoid_metric_is_conformal() {
        // The catenoid is conformally parameterized (g_uu = g_vv, g_uv = 0)
        // when c = 1.  For c ≠ 1 the parameterization is still orthogonal
        // (g_uv = 0) but g_uu = c²·cosh²(v/c) ≠ cosh²(v/c) = g_vv.
        let c = Catenoid::new(1.0);
        for v in [0.0f32, 0.5, 1.0, -0.7] {
            let g = c.metric(0.3, v);
            assert!(
                (g[0][0] - g[1][1]).abs() < 1e-4,
                "not conformal at v={v}: g00={} g11={}",
                g[0][0],
                g[1][1]
            );
            assert!(g[0][1].abs() < 1e-5, "g_uv != 0 at v={v}");
        }
    }

    #[test]
    fn catenoid_christoffel_symmetric() {
        let c = Catenoid::new(1.5);
        let g = c.christoffel(0.3, 0.7);
        for k in 0..2 {
            assert!(
                (g[k][0][1] - g[k][1][0]).abs() < 1e-5,
                "Γ^{k} not symmetric: {} vs {}",
                g[k][0][1],
                g[k][1][0]
            );
        }
    }

    #[test]
    fn catenoid_position_at_waist() {
        // At v=0, cosh(0) = 1, so position = (c*cos(u), c*sin(u), 0).
        let c = Catenoid::new(1.5);
        let p = c.position(0.0, 0.0);
        assert!((p.x - 1.5).abs() < 1e-6, "waist x at u=0: {}", p.x);
        assert!(p.y.abs() < 1e-6);
        assert!(p.z.abs() < 1e-6);
    }
}