tracematch 0.0.2

//! # Geographic Utilities
//!
//! Core geographic computation utilities for GPS track analysis.
//!
//! This module provides fundamental geographic operations used throughout the route matching
//! library. All functions are designed to be efficient and accurate for GPS trajectory data.
//!
//! ## Overview
//!
//! | Function | Description |
//! |----------|-------------|
//! | [`haversine_distance`] | Great-circle distance between two GPS points |
//! | [`compute_bounds`] | Bounding box of a GPS track |
//! | [`compute_center`] | Centroid of a GPS track |
//! | [`bounds_overlap`] | Check if two bounding boxes overlap |
//! | [`meters_to_degrees`] | Convert meters to approximate degrees at a latitude |
//!
//! ## Example
//!
//! ```rust
//! use tracematch::{GpsPoint, geo_utils};
//!
//! let track = vec![
//!     GpsPoint::new(51.5074, -0.1278),  // London
//!     GpsPoint::new(51.5080, -0.1290),
//!     GpsPoint::new(51.5090, -0.1300),
//! ];
//!
//! // Get bounding box
//! let bounds = geo_utils::compute_bounds(&track);
//! println!("Bounds: {:.4}N to {:.4}N", bounds.min_lat, bounds.max_lat);
//!
//! // Distance between two points
//! let dist = geo_utils::haversine_distance(&track[0], &track[2]);
//! println!("Start to end: {:.0}m", dist);
//! ```
//!
//! ## Algorithm Notes
//!
//! ### Haversine Formula
//!
//! The haversine formula calculates the great-circle distance between two points on a sphere.
//! It's the standard method for GPS distance calculation, accurate to within 0.3% for most
//! practical applications.
//!
//! Reference: [Haversine formula (Wikipedia)](https://en.wikipedia.org/wiki/Haversine_formula)
//!
//! ### Coordinate System
//!
//! All functions expect WGS84 coordinates (latitude/longitude in degrees), which is the
//! standard used by GPS receivers and mapping services.

use crate::{Bounds, GpsPoint};
use geo::{Distance, Haversine, Point};

// =============================================================================
// Distance Functions
// =============================================================================

/// Calculate the great-circle distance between two GPS points using the Haversine formula.
///
/// Returns the distance in meters along the Earth's surface (assuming a spherical Earth
/// with radius 6,371 km).
///
/// # Arguments
///
/// * `p1` - First GPS point
/// * `p2` - Second GPS point
///
/// # Returns
///
/// Distance in meters between the two points.
///
/// # Example
///
/// ```rust
/// use tracematch::{GpsPoint, geo_utils};
///
/// let london = GpsPoint::new(51.5074, -0.1278);
/// let paris = GpsPoint::new(48.8566, 2.3522);
///
/// let distance = geo_utils::haversine_distance(&london, &paris);
/// assert!((distance - 343_560.0).abs() < 1000.0); // ~344 km
/// ```
///
/// # Performance
///
/// This function is O(1) and involves trigonometric operations. For comparing distances
/// where exact values aren't needed, consider using squared Euclidean distance on
/// projected coordinates for better performance.
#[inline]
pub fn haversine_distance(p1: &GpsPoint, p2: &GpsPoint) -> f64 {
    let point1 = Point::new(p1.longitude, p1.latitude);
    let point2 = Point::new(p2.longitude, p2.latitude);
    Haversine::distance(point1, point2)
}

/// Convert meters to approximate degrees at a given latitude.
///
/// Uses the WGS84 ellipsoid approximation for latitude-dependent conversion.
/// More accurate at the given latitude than a fixed conversion factor.
///
/// # Arguments
///
/// * `meters` - Distance in meters to convert
/// * `latitude` - Reference latitude for the conversion (in degrees)
///
/// # Returns
///
/// Approximate distance in degrees.
///
/// # Notes
///
/// - At the equator, 1 degree ≈ 111,320 meters
/// - At 45°N/S, 1 degree ≈ 78,710 meters (longitude) / 111,132 meters (latitude)
/// - At the poles, longitude degrees become meaningless
///
/// This function returns a single value suitable for bounding box calculations
/// where a square search area is acceptable.
#[inline]
pub fn meters_to_degrees(meters: f64, latitude: f64) -> f64 {
    // At the equator, 1 degree ≈ 111,320 meters
    // This decreases with cos(latitude) for longitude
    // For simplicity, use a conservative (larger) value based on latitude
    let lat_rad = latitude.to_radians();
    let meters_per_degree = 111_320.0 * lat_rad.cos().max(0.1);
    meters / meters_per_degree
}

// =============================================================================
// Bounding Box Functions
// =============================================================================

/// Compute the bounding box of a GPS track.
///
/// Returns a [`Bounds`] struct containing the minimum and maximum latitude/longitude
/// values that enclose all points in the track.
///
/// # Arguments
///
/// * `points` - Slice of GPS points
///
/// # Returns
///
/// A [`Bounds`] struct with the bounding box coordinates. For empty input,
/// returns a bounds with MIN/MAX values that will fail any overlap check.
///
/// # Example
///
/// ```rust
/// use tracematch::{GpsPoint, geo_utils};
///
/// let track = vec![
///     GpsPoint::new(51.5000, -0.1300),
///     GpsPoint::new(51.5100, -0.1200),
///     GpsPoint::new(51.5050, -0.1250),
/// ];
///
/// let bounds = geo_utils::compute_bounds(&track);
/// assert_eq!(bounds.min_lat, 51.5000);
/// assert_eq!(bounds.max_lat, 51.5100);
/// assert_eq!(bounds.min_lng, -0.1300);
/// assert_eq!(bounds.max_lng, -0.1200);
/// ```
pub fn compute_bounds(points: &[GpsPoint]) -> Bounds {
    let mut min_lat = f64::MAX;
    let mut max_lat = f64::MIN;
    let mut min_lng = f64::MAX;
    let mut max_lng = f64::MIN;

    for p in points {
        min_lat = min_lat.min(p.latitude);
        max_lat = max_lat.max(p.latitude);
        min_lng = min_lng.min(p.longitude);
        max_lng = max_lng.max(p.longitude);
    }

    Bounds {
        min_lat,
        max_lat,
        min_lng,
        max_lng,
    }
}

/// Compute the bounding box as a tuple (min_lat, max_lat, min_lng, max_lng).
///
/// This is a convenience function that returns the bounds as a tuple instead
/// of a [`Bounds`] struct. Useful for quick destructuring.
///
/// # Arguments
///
/// * `points` - Slice of GPS points
///
/// # Returns
///
/// Tuple of (min_lat, max_lat, min_lng, max_lng).
#[inline]
pub fn compute_bounds_tuple(points: &[GpsPoint]) -> (f64, f64, f64, f64) {
    let bounds = compute_bounds(points);
    (
        bounds.min_lat,
        bounds.max_lat,
        bounds.min_lng,
        bounds.max_lng,
    )
}

/// Check if two bounding boxes overlap, with an optional buffer.
///
/// Useful for quick spatial filtering before expensive point-by-point comparisons.
/// Two tracks with non-overlapping bounds cannot share any common points.
///
/// # Arguments
///
/// * `a` - First bounding box
/// * `b` - Second bounding box
/// * `buffer_meters` - Buffer distance in meters to expand the overlap check
/// * `reference_lat` - Reference latitude for meter-to-degree conversion
///
/// # Returns
///
/// `true` if the bounding boxes overlap (including buffer), `false` otherwise.
///
/// # Example
///
/// ```rust
/// use tracematch::{Bounds, geo_utils};
///
/// let bounds_a = Bounds {
///     min_lat: 51.50, max_lat: 51.51,
///     min_lng: -0.13, max_lng: -0.12,
/// };
///
/// let bounds_b = Bounds {
///     min_lat: 51.505, max_lat: 51.515,
///     min_lng: -0.125, max_lng: -0.115,
/// };
///
/// // These bounds overlap
/// assert!(geo_utils::bounds_overlap(&bounds_a, &bounds_b, 0.0, 51.5));
///
/// // With a large negative buffer, they might not
/// // (negative buffer shrinks the overlap zone)
/// ```
pub fn bounds_overlap(a: &Bounds, b: &Bounds, buffer_meters: f64, reference_lat: f64) -> bool {
    let buffer_deg = meters_to_degrees(buffer_meters, reference_lat);

    !(a.max_lat + buffer_deg < b.min_lat
        || b.max_lat + buffer_deg < a.min_lat
        || a.max_lng + buffer_deg < b.min_lng
        || b.max_lng + buffer_deg < a.min_lng)
}

// =============================================================================
// Center/Centroid Functions
// =============================================================================

/// Compute the geographic center (centroid) of a GPS track.
///
/// Returns the arithmetic mean of all latitude and longitude values.
/// This is a simple centroid calculation suitable for small geographic areas.
///
/// # Arguments
///
/// * `points` - Slice of GPS points
///
/// # Returns
///
/// A [`GpsPoint`] at the center of the track. Returns (0, 0) for empty input.
///
/// # Notes
///
/// For tracks spanning large areas or crossing the antimeridian (180°/-180° longitude),
/// this simple averaging may produce unexpected results. For such cases, consider
/// using a proper spherical centroid calculation.
///
/// # Example
///
/// ```rust
/// use tracematch::{GpsPoint, geo_utils};
///
/// let track = vec![
///     GpsPoint::new(51.50, -0.10),
///     GpsPoint::new(51.52, -0.12),
/// ];
///
/// let center = geo_utils::compute_center(&track);
/// assert!((center.latitude - 51.51).abs() < 0.001);
/// assert!((center.longitude - (-0.11)).abs() < 0.001);
/// ```
pub fn compute_center(points: &[GpsPoint]) -> GpsPoint {
    if points.is_empty() {
        return GpsPoint::new(0.0, 0.0);
    }

    let sum_lat: f64 = points.iter().map(|p| p.latitude).sum();
    let sum_lng: f64 = points.iter().map(|p| p.longitude).sum();
    let n = points.len() as f64;

    GpsPoint::new(sum_lat / n, sum_lng / n)
}

// =============================================================================
// Bearing and Gradient Utilities
// =============================================================================

/// Calculate the initial bearing (forward azimuth) from p1 to p2.
///
/// Returns bearing in degrees (0-360), where:
/// - 0° = North
/// - 90° = East
/// - 180° = South
/// - 270° = West
///
/// Uses the spherical law of cosines for azimuth calculation.
#[inline]
pub fn calculate_bearing(p1: &GpsPoint, p2: &GpsPoint) -> f64 {
    let lat1 = p1.latitude.to_radians();
    let lat2 = p2.latitude.to_radians();
    let delta_lng = (p2.longitude - p1.longitude).to_radians();

    let x = delta_lng.sin() * lat2.cos();
    let y = lat1.cos() * lat2.sin() - lat1.sin() * lat2.cos() * delta_lng.cos();

    (x.atan2(y).to_degrees() + 360.0) % 360.0
}

/// Calculate the absolute angular difference between two bearings.
///
/// Returns a value in the range [0, 180] degrees, handling the wraparound at 360°.
#[inline]
pub fn bearing_difference(b1: f64, b2: f64) -> f64 {
    let diff = (b1 - b2).abs();
    if diff > 180.0 {
        360.0 - diff
    } else {
        diff
    }
}

/// Calculate the circular mean of a set of bearings.
///
/// Handles the circular nature of bearings (e.g., average of 350° and 10° should be 0°).
/// Returns a bearing in the range [0, 360).
pub fn circular_mean_bearing(bearings: &[f64]) -> f64 {
    if bearings.is_empty() {
        return 0.0;
    }

    let sum_sin: f64 = bearings.iter().map(|b| b.to_radians().sin()).sum();
    let sum_cos: f64 = bearings.iter().map(|b| b.to_radians().cos()).sum();

    (sum_sin.atan2(sum_cos).to_degrees() + 360.0) % 360.0
}

/// Calculate the circular standard deviation of bearings.
///
/// Measures how spread out the bearings are. A value near 0 indicates
/// consistent direction; higher values indicate more variation.
pub fn circular_std_bearing(bearings: &[f64]) -> f64 {
    if bearings.len() < 2 {
        return 0.0;
    }

    let mean = circular_mean_bearing(bearings);
    let sum_sq: f64 = bearings
        .iter()
        .map(|b| bearing_difference(*b, mean).powi(2))
        .sum();

    (sum_sq / bearings.len() as f64).sqrt()
}

/// Calculate gradient (grade %) between two points.
///
/// Returns `None` if either point lacks elevation data.
/// Positive values indicate uphill (ascending), negative indicates downhill.
///
/// Formula: gradient = (elevation_change / horizontal_distance) * 100
pub fn calculate_gradient(p1: &GpsPoint, p2: &GpsPoint) -> Option<f64> {
    let elev1 = p1.elevation?;
    let elev2 = p2.elevation?;

    let horizontal_dist = haversine_distance(p1, p2);
    if horizontal_dist < 1.0 {
        // Too close - avoid division by very small numbers
        return Some(0.0);
    }

    let elevation_change = elev2 - elev1;
    Some((elevation_change / horizontal_dist) * 100.0)
}

/// Compute average gradient over a segment of points.
///
/// Returns `None` if insufficient elevation data is available.
pub fn segment_gradient(points: &[GpsPoint]) -> Option<f64> {
    if points.len() < 2 {
        return None;
    }

    let first = points.first()?;
    let last = points.last()?;

    let elev_start = first.elevation?;
    let elev_end = last.elevation?;

    // Calculate total horizontal distance
    let mut total_dist = 0.0;
    for i in 1..points.len() {
        total_dist += haversine_distance(&points[i - 1], &points[i]);
    }

    if total_dist < 1.0 {
        return Some(0.0);
    }

    Some(((elev_end - elev_start) / total_dist) * 100.0)
}