use super::funejson::Value;
use super::geo_point_in_polygon::{geo_point_from_value, GeoPoint};
use super::math_round_float::round_float;
use super::math_sin_cos::sin_cos;

const RADIUS: f64 = 6371008.8;
const RADIANS: f64 = 3.141592653589793 / 180.0;
const ACRE: f64 = 4046.8564224;

/// Area enclosed by a boundary, by the shoelace formula in the Lambert
/// cylindrical equal-area projection (x = R lng, y = R sin lat), which is the
/// spherical polygon formula of Chamberlain and Duquette (2007).
///
/// Longitudes are unwrapped from the first point and sines are taken relative
/// to the first point's, which keeps a one-hectare field accurate when the
/// coordinates are large; sin comes from math.sin-cos so every language adds
/// the same doubles in the same order.
///
/// # Panics
/// Panics on a coordinate out of range or fewer than 3 distinct points.
pub fn field_area(boundary: &[GeoPoint]) -> FieldArea {
    for (i, p) in boundary.iter().enumerate() {
        if !p.lat.is_finite() || p.lat < -90.0 || p.lat > 90.0 {
            panic!("point {}: latitude must be a finite number from -90 to 90 degrees", i + 1);
        }
        if !p.lng.is_finite() || p.lng < -180.0 || p.lng > 180.0 {
            panic!("point {}: longitude must be a finite number from -180 to 180 degrees", i + 1);
        }
    }
    let mut ring: &[GeoPoint] = boundary;
    if ring.len() > 1 && ring[0].lat == ring[ring.len() - 1].lat && ring[0].lng == ring[ring.len() - 1].lng {
        ring = &ring[..ring.len() - 1];
    }
    let mut distinct: Vec<(f64, f64)> = Vec::new();
    for p in ring {
        if !distinct.iter().any(|&(a, b)| a == p.lat && b == p.lng) {
            distinct.push((p.lat, p.lng));
        }
    }
    if distinct.len() < 3 {
        panic!("a boundary needs at least 3 distinct points");
    }

    let n = ring.len();
    let mut x: Vec<f64> = Vec::with_capacity(n);
    let mut y: Vec<f64> = Vec::with_capacity(n);
    let sin0 = sin_cos(ring[0].lat * RADIANS).sin;
    let mut previous = 0.0;
    for p in ring {
        let mut d = p.lng - ring[0].lng;
        // Take the short way round, so a field across the antimeridian works.
        while d - previous > 180.0 {
            d -= 360.0;
        }
        while d - previous < -180.0 {
            d += 360.0;
        }
        previous = d;
        x.push(d * RADIANS);
        y.push(sin_cos(p.lat * RADIANS).sin - sin0);
    }
    let mut total = 0.0;
    for i in 0..n {
        total += (x[(i + 1) % n] - x[(i + n - 1) % n]) * y[i];
    }
    let area = total.abs() * RADIUS * RADIUS / 2.0;
    FieldArea {
        square_metres: round_float(area, 0),
        hectares: round_float(area / 10000.0, 4),
        acres: round_float(area / ACRE, 4),
    }
}

pub fn field_area_to_value(f: &FieldArea) -> Value {
    Value::obj(vec![
        ("squareMetres", Value::Float(f.square_metres)),
        ("hectares", Value::Float(f.hectares)),
        ("acres", Value::Float(f.acres)),
    ])
}

pub fn fune_vector(args: &[Value]) -> Value {
    let boundary: Vec<GeoPoint> = args[0].as_arr().iter().map(geo_point_from_value).collect();
    field_area_to_value(&field_area(&boundary))
}