use nalgebra::*;

pub struct Circle {
  pub center_x: f64,
  pub center_y: f64,
  pub radius: f64,
}

// construct the circle through the points given by the columns of `points`
pub fn circ_thru(points: Matrix2x3<f64>) -> Option<Circle> {
    // build the matrix that maps the circle's coefficient vector to the
    // negative of the linear part of the circle's equation, evaluated at the
    // given points
    let neg_lin_part = stack![2.0*points.transpose(), Vector3::repeat(1.0)];
    
    // find the quadrdatic part of the circle's equation, evaluated at the given
    // points
    let quad_part = Vector3::from_iterator(
      points.column_iter().map(|v| v.dot(&v))
    );
    
    // find the circle's coefficient vector, and from there its center and
    // radius
    match neg_lin_part.lu().solve(&quad_part) {
        None => None,
        Some(coeffs) => {
            let center_x = coeffs[0];
            let center_y = coeffs[1];
            Some(Circle {
                center_x: center_x,
                center_y: center_y,
                radius: (
                    coeffs[2] + center_x*center_x + center_y*center_y
                ).sqrt(),
            })
        }
    }
}