# Author: Simon-Pierre Boucher — contact@spboucher.ai
"""Geospatial helpers."""

import numpy as np

from src.config import EARTH_RADIUS_KM


def haversine_matrix(lat1, lon1, lat2, lon2):
    """
    Vectorised Haversine distance between every pair (i, j) where
    i indexes lat1/lon1 and j indexes lat2/lon2.

    Parameters
    ----------
    lat1, lon1 : np.ndarray, shape (n,)
    lat2, lon2 : np.ndarray, shape (m,)

    Returns
    -------
    dist : np.ndarray, shape (n, m) — distances in kilometres
    """
    lat1 = np.radians(lat1)[:, None]
    lon1 = np.radians(lon1)[:, None]
    lat2 = np.radians(lat2)[None, :]
    lon2 = np.radians(lon2)[None, :]

    dlat = lat2 - lat1
    dlon = lon2 - lon1

    a = np.sin(dlat / 2.0) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2.0) ** 2
    c = 2.0 * np.arcsin(np.sqrt(a))
    return EARTH_RADIUS_KM * c