monge-rs 0.1.0

//! Radical axis as a 1-cocycle in the Čech cohomology of circle covers.
//!
//! ## Cocycle interpretation
//!
//! For each circle `i`, define the **power function** pᵢ(x) = |x − cᵢ|² − rᵢ².
//! The **radical axis** of circles i and j is the set where pᵢ = pⱼ.
//!
//! In cohomological terms:
//! - Let {Uᵢ} be the open cover by disk interiors.
//! - On overlaps Uᵢ ∩ Uⱼ, define f_{ij}(x) = pᵢ(x) − pⱼ(x).
//! - This is a **1-cocycle**: f_{ij} + f_{jk} = f_{ik} on triple overlaps.
//! - The radical axis is the zero set of f_{ij}.
//!
//! For Monge's theorem, the cocycle condition f_{12} + f_{23} = f_{13}
//! is exactly the cohomological reason that the three external homothetic
//! centers are collinear.

use nalgebra::Vector2;
use crate::geometry::Circle;

/// Represents the radical axis (line) between two circles.
///
/// The axis is stored as coefficients (a, b, c) of the line equation
/// `a·x + b·y + c = 0`.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct RadicalAxis {
    pub a: f64,
    pub b: f64,
    pub c: f64,
}

impl RadicalAxis {
    /// Create a new radical axis from circles `c1` and `c2`.
    ///
    /// The radical axis is the set of points where the power functions
    /// of the two circles are equal.  This is always a line (or empty
    /// for concentric circles with equal radii).
    pub fn new(c1: &Circle, c2: &Circle) -> Self {
        let a = 2.0 * (c2.center.x - c1.center.x);
        let b = 2.0 * (c2.center.y - c1.center.y);
        let c = c1.center.norm_squared()
            - c1.radius.powi(2)
            - c2.center.norm_squared()
            + c2.radius.powi(2);
        Self { a, b, c }
    }

    /// Evaluate the line equation at point `x`.
    pub fn eval(&self, x: &Vector2<f64>) -> f64 {
        self.a * x.x + self.b * x.y + self.c
    }

    /// Check if point `x` lies on the radical axis (within `eps` tolerance).
    pub fn contains(&self, x: &Vector2<f64>, eps: f64) -> bool {
        self.eval(x).abs() < eps
    }

    /// Get the slope of this line (None for vertical lines).
    pub fn slope(&self) -> Option<f64> {
        if self.b.abs() < f64::EPSILON {
            None
        } else {
            Some(-self.a / self.b)
        }
    }

    /// Intersection point of two radical axes.
    /// Returns `None` if they are parallel.
    pub fn intersect(&self, other: &RadicalAxis) -> Option<Vector2<f64>> {
        let det = self.a * other.b - self.b * other.a;
        if det.abs() < f64::EPSILON * 1e4 {
            return None; // Parallel
        }
        Some(Vector2::new(
            (self.b * other.c - self.c * other.b) / det,
            (self.c * other.a - self.a * other.c) / det,
        ))
    }

    /// Check the cocycle condition: f_{ij} + f_{jk} = f_{ik}.
    ///
    /// For the radical axes of three circles, this must hold at all points.
    pub fn check_cocycle(c1: &Circle, c2: &Circle, c3: &Circle, at: &Vector2<f64>) -> bool {
        let r12 = Self::new(c1, c2);
        let r23 = Self::new(c2, c3);
        let r13 = Self::new(c1, c3);
        let lhs = r12.eval(at) + r23.eval(at);
        let rhs = r13.eval(at);
        (lhs - rhs).abs() < 1e-10
    }
}

/// The radical center of three circles — the unique point where all three
/// radical axes intersect.  Returns `None` if two are parallel.
pub fn radical_center(
    c1: &Circle,
    c2: &Circle,
    c3: &Circle,
) -> Option<Vector2<f64>> {
    let r12 = RadicalAxis::new(c1, c2);
    let r23 = RadicalAxis::new(c2, c3);
    // The three axes intersect at a common point (the radical center)
    r12.intersect(&r23)
}

/// Compute the 1-cocycle condition for three circles at point `x`.
///
/// Returns (f_{12}(x), f_{23}(x), f_{13}(x)) where f_{ij} = p_i - p_j.
/// The cocycle condition is f_{12} + f_{23} = f_{13}.
pub fn cocycle_values(
    c1: &Circle,
    c2: &Circle,
    c3: &Circle,
    x: &Vector2<f64>,
) -> (f64, f64, f64) {
    let f12 = c1.power(x) - c2.power(x);
    let f23 = c2.power(x) - c3.power(x);
    let f13 = c1.power(x) - c3.power(x);
    (f12, f23, f13)
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

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

    #[test]
    fn test_radical_axis_basic() {
        // Two circles centered on x-axis
        let c1 = Circle::new(Vector2::new(0.0, 0.0), 1.0);
        let c2 = Circle::new(Vector2::new(4.0, 0.0), 2.0);
        let axis = RadicalAxis::new(&c1, &c2);
        // Line: 8x - 13 = 0 => x = 13/8 = 1.625
        let test_pt = Vector2::new(13.0 / 8.0, 0.0);
        assert!(axis.contains(&test_pt, 1e-10));
    }

    #[test]
    fn test_concentric_equal() {
        // Concentric, equal radii => degenerate line (0=0 everywhere)
        let c1 = Circle::new(Vector2::new(0.0, 0.0), 1.0);
        let c2 = Circle::new(Vector2::new(0.0, 0.0), 1.0);
        let axis = RadicalAxis::new(&c1, &c2);
        // a=0, b=0, c=0 — every point satisfies
        assert!((axis.a).abs() < 1e-12);
        assert!((axis.b).abs() < 1e-12);
        assert!((axis.c).abs() < 1e-12);
    }

    #[test]
    fn test_intersecting_axes() {
        // Three circles arranged so radical axes intersect at a common point
        let c1 = Circle::new(Vector2::new(0.0, 0.0), 1.0);
        let c2 = Circle::new(Vector2::new(4.0, 0.0), 2.0);
        let c3 = Circle::new(Vector2::new(2.0, 3.0), 1.5);

        let center = radical_center(&c1, &c2, &c3).unwrap();
        let r12 = RadicalAxis::new(&c1, &c2);
        let r23 = RadicalAxis::new(&c2, &c3);
        let r31 = RadicalAxis::new(&c3, &c1);

        assert!(r12.contains(&center, 1e-8));
        assert!(r23.contains(&center, 1e-8));
        assert!(r31.contains(&center, 1e-8));
    }

    #[test]
    fn test_cocycle_condition() {
        let c1 = Circle::new(Vector2::new(0.0, 0.0), 1.0);
        let c2 = Circle::new(Vector2::new(4.0, 0.0), 2.0);
        let c3 = Circle::new(Vector2::new(2.0, 3.0), 1.5);

        // Test at several arbitrary points
        for pt in &[
            Vector2::new(0.0, 0.0),
            Vector2::new(1.0, 1.0),
            Vector2::new(-3.0, 7.0),
        ] {
            assert!(
                RadicalAxis::check_cocycle(&c1, &c2, &c3, pt),
                "Cocycle condition should hold at {:?}",
                pt
            );
        }
    }

    #[test]
    fn test_cocycle_values() {
        let c1 = Circle::new(Vector2::new(0.0, 0.0), 1.0);
        let c2 = Circle::new(Vector2::new(4.0, 0.0), 2.0);
        let c3 = Circle::new(Vector2::new(2.0, 3.0), 1.5);
        let x = Vector2::new(1.0, 2.0);

        let (f12, f23, f13) = cocycle_values(&c1, &c2, &c3, &x);
        assert!((f12 + f23 - f13).abs() < 1e-10, "Cocycle f12 + f23 = f13");
    }
}