monge-rs 0.1.0

//! External and internal homothetic centers of circles.
//!
//! For two circles (c₁, r₁) and (c₂, r₂):
//! - **External center**: divides centers externally in ratio of radii.
//!   Formula: `S₁₂ = (r₂·c₁ − r₁·c₂) / (r₂ − r₁)` for r₁ ≠ r₂.
//!   When r₁ = r₂, the external center is at infinity (returns `None`).
//! - **Internal center**: divides centers internally in ratio of radii.
//!   Formula: `S₁₂^ⱽ = (r₂·c₁ + r₁·c₂) / (r₂ + r₁)`.
//!
//! In Monge's theorem, the three *external* homothetic centers of three
//! circles are collinear, and the three *internal* centers lie on the same
//! line as each pair's corresponding external center.

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

/// Compute the external homothetic center of two circles.
///
/// Returns `None` when r₁ ≈ r₂ (the center is at infinity).
///
/// ## Formula
/// ```text
/// S₁₂ = (r₂·c₁ − r₁·c₂) / (r₂ − r₁)
/// ```
pub fn external_homothetic_center(c1: &Circle, c2: &Circle) -> Option<Vector2<f64>> {
    let dr = c2.radius - c1.radius;
    if dr.abs() < f64::EPSILON * 1e4 {
        // Radii are effectively equal — center at infinity
        return None;
    }
    let result = (c2.radius * c1.center - c1.radius * c2.center) / dr;
    Some(result)
}

/// Compute the internal homothetic center of two circles.
///
/// This always exists (for non-degenerate circles).
///
/// ## Formula
/// ```text
/// S₁₂^ⱽ = (r₂·c₁ + r₁·c₂) / (r₂ + r₁)
/// ```
pub fn internal_homothetic_center(c1: &Circle, c2: &Circle) -> Vector2<f64> {
    let sum_r = c2.radius + c1.radius;
    if sum_r.abs() < f64::EPSILON {
        // Both radii are zero — return midpoint
        return (c1.center + c2.center) / 2.0;
    }
    (c2.radius * c1.center + c1.radius * c2.center) / sum_r
}

/// Compute all three external homothetic centers of three circles.
///
/// Returns `(S₁₂, S₂₃, S₃₁)` — the three centers that Monge's theorem
/// asserts are collinear.
///
/// Returns `None` if any pair has equal radii (degenerate case).
pub fn external_triple(
    c1: &Circle,
    c2: &Circle,
    c3: &Circle,
) -> Option<(Vector2<f64>, Vector2<f64>, Vector2<f64>)> {
    let s12 = external_homothetic_center(c1, c2)?;
    let s23 = external_homothetic_center(c2, c3)?;
    let s31 = external_homothetic_center(c3, c1)?;
    Some((s12, s23, s31))
}

/// Compute all three internal homothetic centers of three circles.
pub fn internal_triple(
    c1: &Circle,
    c2: &Circle,
    c3: &Circle,
) -> (Vector2<f64>, Vector2<f64>, Vector2<f64>) {
    let s12 = internal_homothetic_center(c1, c2);
    let s23 = internal_homothetic_center(c2, c3);
    let s31 = internal_homothetic_center(c3, c1);
    (s12, s23, s31)
}

/// Verify Monge's theorem: the three external homothetic centers are collinear.
///
/// Returns the area of the triangle formed by the three centers.
/// If they are collinear, this area should be zero (within numerical tolerance).
pub fn monge_collinear_area(
    c1: &Circle,
    c2: &Circle,
    c3: &Circle,
) -> Result<f64, &'static str> {
    let (s12, s23, s31) = external_triple(c1, c2, c3)
        .ok_or("Equal radii — external centers at infinity")?;
    Ok(triangle_area(&s12, &s23, &s31))
}

/// Signed area of triangle (p, q, r). Zero iff points are collinear.
pub fn triangle_area(p: &Vector2<f64>, q: &Vector2<f64>, r: &Vector2<f64>) -> f64 {
    0.5 * ((q.x - p.x) * (r.y - p.y) - (q.y - p.y) * (r.x - p.x)).abs()
}

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

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

    #[test]
    fn test_external_homothetic_basic() {
        // Circles at (0,0) r=1 and (4,0) r=2
        let c1 = Circle::new(Vector2::new(0.0, 0.0), 1.0);
        let c2 = Circle::new(Vector2::new(4.0, 0.0), 2.0);
        // S12 = (2*0 - 1*4)/(2-1) = (-4)/1 = (-4, 0)
        let s = external_homothetic_center(&c1, &c2).unwrap();
        assert!((s.x - (-4.0)).abs() < 1e-10);
        assert!((s.y).abs() < 1e-10);
    }

    #[test]
    fn test_internal_homothetic_basic() {
        let c1 = Circle::new(Vector2::new(0.0, 0.0), 1.0);
        let c2 = Circle::new(Vector2::new(4.0, 0.0), 2.0);
        // S12^V = (2*0 + 1*4)/(2+1) = 4/3 ≈ 1.333
        let s = internal_homothetic_center(&c1, &c2);
        assert!((s.x - 4.0 / 3.0).abs() < 1e-10);
        assert!((s.y).abs() < 1e-10);
    }

    #[test]
    fn test_external_equal_radii() {
        let c1 = Circle::new(Vector2::new(0.0, 0.0), 1.0);
        let c2 = Circle::new(Vector2::new(4.0, 0.0), 1.0);
        // Should be None (parallel lines, center at infinity)
        assert!(external_homothetic_center(&c1, &c2).is_none());
    }

    #[test]
    fn test_monge_theorem() {
        // Three circles — all radii distinct
        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 area = monge_collinear_area(&c1, &c2, &c3).unwrap();
        assert!(
            area < 1e-10,
            "Monge centers should be collinear, area = {}",
            area
        );
    }

    #[test]
    fn test_internal_center_collinear_with_circle_centers() {
        // Each internal homothetic center lies on the line connecting
        // its two circle centers (but the three internal centers are
        // NOT necessarily collinear with each other)
        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 s = internal_homothetic_center(&c1, &c2);
        // Internal center lies between the two centers
        let area = triangle_area(&c1.center, &c2.center, &s);
        assert!(
            area < 1e-10,
            "Internal center should be collinear with circle centers, area = {}",
            area
        );
    }

    #[test]
    fn test_equal_radii_edge_case() {
        // When two circles have equal radii, the external center is at infinity
        let c1 = Circle::new(Vector2::new(0.0, 0.0), 2.0);
        let c2 = Circle::new(Vector2::new(4.0, 0.0), 2.0);
        assert!(external_homothetic_center(&c1, &c2).is_none());

        // Internal still works
        let s = internal_homothetic_center(&c1, &c2);
        assert!((s.x - 2.0).abs() < 1e-10);

        // Verify collinearity of internal center with circle centers
        let area = triangle_area(&c1.center, &c2.center, &s);
        assert!(area < 1e-10);
    }
}