"""
Distance calculations using the haversine formula.

Supports both scalar and batch (NumPy array) operations.
"""

import numpy as np
from typing import Tuple

EARTH_RADIUS = 6371000  # meters


def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
    """Calculate distance between two points in meters."""
    lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
    c = 2 * np.arcsin(np.sqrt(a))
    return float(EARTH_RADIUS * c)


def haversine_batch(
    lat1: float, lon1: float,
    lats2: np.ndarray, lons2: np.ndarray,
) -> np.ndarray:
    """Distance from one point to many points (meters)."""
    lat1_rad = np.radians(lat1)
    lon1_rad = np.radians(lon1)
    lats2_rad = np.radians(lats2)
    lons2_rad = np.radians(lons2)

    dlat = lats2_rad - lat1_rad
    dlon = lons2_rad - lon1_rad

    a = np.sin(dlat / 2) ** 2 + np.cos(lat1_rad) * np.cos(lats2_rad) * np.sin(dlon / 2) ** 2
    c = 2 * np.arcsin(np.sqrt(a))

    return EARTH_RADIUS * c


def points_to_local_coords(
    ref_lat: float, ref_lon: float,
    lats: np.ndarray, lons: np.ndarray,
) -> Tuple[np.ndarray, np.ndarray]:
    """Convert lat/lon to local X/Y meters (equirectangular projection)."""
    cos_lat = np.cos(np.radians(ref_lat))
    x = (lons - ref_lon) * 111320.0 * cos_lat
    y = (lats - ref_lat) * 110540.0
    return x, y