# -*- coding: utf-8 -*-
"""
Wildlife Sound Locator — Geodesic Computation Engine.

Provides:
    1. Destination point calculation (bearing + distance → target)
    2. Confidence sector polygon generation
    3. Bearing line generation (observer → estimated point)
    4. Two-observer triangulation (intersection of two bearings)

All calculations use WGS-84 ellipsoid with Vincenty-style formulas
for accuracy at typical forest patrol distances (10 m – 5 km).

Author: Lộc Vũ Trung (LVT) / Slow Forest
License: GPL-3.0

Scientific references:
    - Vincenty, T. (1975). Direct / inverse geodetic solutions.
    - Movable Type Scripts (movable-type.co.uk / scripts / latlong.html)
"""

import math

# WGS-84 Earth radius (mean) in meters
_R = 6_371_000.0


def destination_point(lat, lon, bearing_deg, distance_m):
    """Calculate destination point given start, bearing, and distance.

    Uses the spherical law of cosines (Haversine-family).

    Args:
        lat: Latitude of start point in decimal degrees.
        lon: Longitude of start point in decimal degrees.
        bearing_deg: Azimuth from North (0–360°), clockwise.
        distance_m: Distance in meters.

    Returns:
        (lat2, lon2) in decimal degrees.
    """
    lat1 = math.radians(lat)
    lon1 = math.radians(lon)
    brng = math.radians(bearing_deg)
    d_r = distance_m / _R

    lat2 = math.asin(
        math.sin(lat1) * math.cos(d_r) +  # noqa: W504
        math.cos(lat1) * math.sin(d_r) * math.cos(brng)
    )
    lon2 = lon1 + math.atan2(
        math.sin(brng) * math.sin(d_r) * math.cos(lat1),
        math.cos(d_r) - math.sin(lat1) * math.sin(lat2),
    )

    return math.degrees(lat2), math.degrees(lon2)


def bearing_line_coords(lat, lon, bearing_deg, distance_m, segments=2):
    """Generate coordinates for a bearing line (observer → target).

    Args:
        lat, lon: Observer position in decimal degrees.
        bearing_deg: Azimuth from North.
        distance_m: Distance to target.
        segments: Number of line segments (2 = straight line).

    Returns:
        List of (lon, lat) tuples (GIS convention: x, y).
    """
    coords = [(lon, lat)]  # start = observer
    for i in range(1, segments + 1):
        frac = i / segments
        d = distance_m * frac
        pt_lat, pt_lon = destination_point(lat, lon, bearing_deg, d)
        coords.append((pt_lon, pt_lat))
    return coords


def confidence_sector_coords(lat, lon, bearing_deg, distance_m,
                             error_angle_deg=15, error_distance_pct=30,
                             arc_steps=20):
    """Generate a sector polygon representing the confidence zone.

    The sector is defined by:
        - Angular range: [bearing - error_angle, bearing + error_angle]
        - Distance range: [distance * (1 - error_pct%), distance * (1 + error_pct%)]

    Args:
        lat, lon: Observer position in decimal degrees.
        bearing_deg: Central azimuth.
        distance_m: Central distance estimate.
        error_angle_deg: Half-width of angular uncertainty (±degrees).
        error_distance_pct: Distance uncertainty (±percent).
        arc_steps: Number of points along each arc.

    Returns:
        List of (lon, lat) tuples forming a closed polygon.
    """
    b_min = bearing_deg - error_angle_deg
    b_max = bearing_deg + error_angle_deg
    d_min = distance_m * (1 - error_distance_pct / 100.0)
    d_max = distance_m * (1 + error_distance_pct / 100.0)
    d_min = max(d_min, 1.0)  # avoid zero / negative

    coords = []

    # Origin (observer point) — include for visual fan shape
    coords.append((lon, lat))

    # Left edge: from observer to outer arc at b_min
    outer_lat, outer_lon = destination_point(lat, lon, b_min, d_max)
    coords.append((outer_lon, outer_lat))

    # Outer arc: b_min → b_max at d_max
    for i in range(1, arc_steps + 1):
        angle = b_min + (b_max - b_min) * i / arc_steps
        pt_lat, pt_lon = destination_point(lat, lon, angle, d_max)
        coords.append((pt_lon, pt_lat))

    # Close back to observer
    coords.append((lon, lat))

    return coords


def confidence_sector_ring_coords(lat, lon, bearing_deg, distance_m,
                                  error_angle_deg=15, error_distance_pct=30,
                                  arc_steps=20):
    """Generate a ring-sector (donut sector) polygon.

    Same as confidence_sector_coords but with inner / outer arcs,
    showing only the annular sector (not from observer).

    Returns:
        List of (lon, lat) tuples forming a closed polygon.
    """
    b_min = bearing_deg - error_angle_deg
    b_max = bearing_deg + error_angle_deg
    d_min = distance_m * (1 - error_distance_pct / 100.0)
    d_max = distance_m * (1 + error_distance_pct / 100.0)
    d_min = max(d_min, 1.0)

    coords = []

    # Outer arc: b_min → b_max at d_max
    for i in range(arc_steps + 1):
        angle = b_min + (b_max - b_min) * i / arc_steps
        pt_lat, pt_lon = destination_point(lat, lon, angle, d_max)
        coords.append((pt_lon, pt_lat))

    # Inner arc: b_max → b_min at d_min (reversed)
    for i in range(arc_steps, -1, -1):
        angle = b_min + (b_max - b_min) * i / arc_steps
        pt_lat, pt_lon = destination_point(lat, lon, angle, d_min)
        coords.append((pt_lon, pt_lat))

    # Close polygon
    if coords:
        coords.append(coords[0])

    return coords


def triangulate(lat1, lon1, bearing1_deg, lat2, lon2, bearing2_deg):
    """Find intersection of two bearings from two observer points.

    Uses the spherical intersection formula.

    Args:
        lat1, lon1: Observer A position (decimal degrees).
        bearing1_deg: Bearing from Observer A (degrees from North).
        lat2, lon2: Observer B position (decimal degrees).
        bearing2_deg: Bearing from Observer B (degrees from North).

    Returns:
        (lat, lon) of intersection point, or None if lines are
        parallel / do not converge.
    """
    phi1 = math.radians(lat1)
    lam1 = math.radians(lon1)
    phi2 = math.radians(lat2)
    lam2 = math.radians(lon2)
    theta13 = math.radians(bearing1_deg)
    theta23 = math.radians(bearing2_deg)

    d_phi = phi2 - phi1
    d_lam = lam2 - lam1

    # Angular distance p1-p2
    delta12 = 2 * math.asin(math.sqrt(
        math.sin(d_phi / 2) ** 2 +  # noqa: W504
        math.cos(phi1) * math.cos(phi2) * math.sin(d_lam / 2) ** 2
    ))
    if abs(delta12) < 1e-12:
        return None  # coincident points

    # Initial/final bearings between points
    cos_theta_a = (math.sin(phi2) - math.sin(phi1) * math.cos(delta12)) / (
        math.sin(delta12) * math.cos(phi1))
    cos_theta_b = (math.sin(phi1) - math.sin(phi2) * math.cos(delta12)) / (
        math.sin(delta12) * math.cos(phi2))

    cos_theta_a = max(-1.0, min(1.0, cos_theta_a))
    cos_theta_b = max(-1.0, min(1.0, cos_theta_b))
    theta_a = math.acos(cos_theta_a)
    theta_b = math.acos(cos_theta_b)

    if math.sin(d_lam) > 0:
        theta12 = theta_a
        theta21 = 2 * math.pi - theta_b
    else:
        theta12 = 2 * math.pi - theta_a
        theta21 = theta_b

    alpha1 = theta13 - theta12
    alpha2 = theta21 - theta23

    sin_a1 = math.sin(alpha1)
    sin_a2 = math.sin(alpha2)

    # Check for parallel / near-parallel bearings
    bearing_diff = abs(bearing1_deg - bearing2_deg) % 360
    if bearing_diff > 180:
        bearing_diff = 360 - bearing_diff
    if bearing_diff < 2.0 or abs(bearing_diff - 180) < 2.0:
        return None  # parallel or anti-parallel bearings

    if sin_a1 == 0 and sin_a2 == 0:
        return None  # infinite intersections (parallel)
    if sin_a1 * sin_a2 < 0:
        return None  # ambiguous intersection

    cos_a1 = math.cos(alpha1)
    cos_a2 = math.cos(alpha2)

    alpha3 = math.acos(
        -cos_a1 * cos_a2 + sin_a1 * sin_a2 * math.cos(delta12)
    )

    delta13 = math.atan2(
        math.sin(delta12) * sin_a1 * sin_a2,
        cos_a2 + cos_a1 * math.cos(alpha3),
    )

    phi3 = math.asin(
        math.sin(phi1) * math.cos(delta13) +  # noqa: W504
        math.cos(phi1) * math.sin(delta13) * math.cos(theta13)
    )

    d_lam13 = math.atan2(
        math.sin(theta13) * math.sin(delta13) * math.cos(phi1),
        math.cos(delta13) - math.sin(phi1) * math.sin(phi3),
    )

    lam3 = lam1 + d_lam13

    return math.degrees(phi3), math.degrees(lam3)


def sector_area_m2(distance_m, error_angle_deg, error_distance_pct):
    """Approximate area of the confidence sector in m².

    Uses the formula for an annular sector:
        A = (θ/360) × π × (R_max² - R_min²)

    Args:
        distance_m: Central distance.
        error_angle_deg: Half-angular uncertainty (total sector = 2×this).
        error_distance_pct: Distance uncertainty percent.

    Returns:
        Area in square meters.
    """
    d_min = distance_m * (1 - error_distance_pct / 100.0)
    d_max = distance_m * (1 + error_distance_pct / 100.0)
    d_min = max(d_min, 0)
    theta_total = 2 * error_angle_deg  # full sector angle
    return (theta_total / 360.0) * math.pi * (d_max ** 2 - d_min ** 2)


def sector_area_ha(distance_m, error_angle_deg, error_distance_pct):
    """Same as sector_area_m2 but returns hectares."""
    return sector_area_m2(distance_m, error_angle_deg, error_distance_pct) / 10000.0
