#!/usr/bin/env python3
"""A two-lamp inverse-optics toy for the Shadow Watch. Standard library only.

All distances are metres, all times seconds. The bird moves in one vertical
plane toward the lighthouse at x=0. Two coded lamps at different heights cast
two distinguishable shadows onto level sand. This is a model, not bird data.
"""

from __future__ import annotations

import argparse
import math
import random
from dataclasses import dataclass


LAMP_A = 40.0
LAMP_B = 30.0
DELTA_T = 2.0


@dataclass(frozen=True)
class Bird:
    range_m: float
    altitude_m: float


def ground_shadow(lamp_height_m: float, bird: Bird) -> float:
    """Ground intercept of the lamp→bird ray, measured from the tower."""
    if lamp_height_m <= 0 or bird.range_m < 0 or not 0 <= bird.altitude_m < lamp_height_m:
        raise ValueError("bird must be between a positive-height lamp and level ground")
    return lamp_height_m * bird.range_m / (lamp_height_m - bird.altitude_m)


def reconstruct(height_a_m: float, shadow_a_m: float,
                height_b_m: float, shadow_b_m: float) -> Bird:
    """Intersect the two lamp→shadow rays without ever seeing the bird."""
    if height_a_m <= 0 or height_b_m <= 0 or height_a_m == height_b_m:
        raise ValueError("lamps need distinct positive heights")
    denominator = shadow_a_m / height_a_m - shadow_b_m / height_b_m
    if abs(denominator) < 1e-12:
        raise ValueError("rays do not provide a finite intersection")
    altitude = (shadow_a_m - shadow_b_m) / denominator
    distance = shadow_a_m * (1 - altitude / height_a_m)
    if not (0 <= altitude < min(height_a_m, height_b_m)) or distance < 0:
        raise ValueError("observations do not locate a bird below both lamps")
    return Bird(distance, altitude)


def arrival_time(first: Bird, second: Bird, dt_s: float = DELTA_T) -> float:
    """Time from the first flash until x=0, assuming constant horizontal speed."""
    if dt_s <= 0:
        raise ValueError("flash separation must be positive")
    approach = first.range_m - second.range_m
    return math.inf if approach <= 0 else dt_s * first.range_m / approach


def single_lamp_arrival(first_shadow_m: float, second_shadow_m: float,
                        dt_s: float = DELTA_T) -> float:
    """Naive shadow-distance / shadow-speed clock; assumes fixed altitude."""
    approach = first_shadow_m - second_shadow_m
    return math.inf if approach <= 0 else dt_s * first_shadow_m / approach


def demo(climb_mps: float = 0.0) -> None:
    first = Bird(152.0, 20.0)
    second = Bird(136.0, 20.0 + climb_mps * DELTA_T)
    readings = [(ground_shadow(LAMP_A, b), ground_shadow(LAMP_B, b)) for b in (first, second)]
    recovered = [reconstruct(LAMP_A, a, LAMP_B, b) for a, b in readings]
    print(f"coded lamps: A={LAMP_A:g} m, B={LAMP_B:g} m; flash separation={DELTA_T:g} s")
    print(f"shadow A: {readings[0][0]:.2f} → {readings[1][0]:.2f} m")
    print(f"shadow B: {readings[0][1]:.2f} → {readings[1][1]:.2f} m")
    print(f"reconstructed bird: {recovered[0]} → {recovered[1]}")
    print(f"one-lamp ETA: {single_lamp_arrival(readings[0][0], readings[1][0]):.2f} s")
    print(f"two-lamp ETA: {arrival_time(*recovered):.2f} s")


def noise_demo(sigma_m: float, trials: int = 20_000) -> None:
    """Seeded Gaussian centroid-error experiment, not a field uncertainty model."""
    if sigma_m <= 0:
        raise ValueError("noise sigma must be positive")
    rng = random.Random(427)
    estimates = []
    for _ in range(trials):
        recovered = []
        for bird in (Bird(152.0, 20.0), Bird(136.0, 20.0)):
            a = ground_shadow(LAMP_A, bird) + rng.gauss(0, sigma_m)
            b = ground_shadow(LAMP_B, bird) + rng.gauss(0, sigma_m)
            try:
                recovered.append(reconstruct(LAMP_A, a, LAMP_B, b))
            except ValueError:
                break
        if len(recovered) == 2:
            eta = arrival_time(*recovered)
            if math.isfinite(eta):
                estimates.append(eta)
    estimates.sort()
    n = len(estimates)
    if not n:
        raise ValueError("all noisy reconstructions were invalid")
    print(f"sigma={sigma_m:g} m; {n}/{trials} finite reconstructions")
    print(f"arrival ETA central 95%: {estimates[int(.025*n)]:.2f}–{estimates[int(.975*n)]:.2f} s")


def selftest() -> None:
    rng = random.Random(427)
    for _ in range(2000):
        x = rng.uniform(40, 240)
        z = rng.uniform(3, 23)
        speed = rng.uniform(2, 15)
        climb = rng.uniform(-1, 1)
        first = Bird(x, z)
        second = Bird(x - speed * DELTA_T, z + climb * DELTA_T)
        for state in (first, second):
            a = ground_shadow(LAMP_A, state)
            b = ground_shadow(LAMP_B, state)
            inferred = reconstruct(LAMP_A, a, LAMP_B, b)
            assert math.isclose(inferred.range_m, state.range_m, rel_tol=1e-11, abs_tol=1e-10)
            assert math.isclose(inferred.altitude_m, state.altitude_m, rel_tol=1e-11, abs_tol=1e-10)
        a0, a1 = ground_shadow(LAMP_A, first), ground_shadow(LAMP_A, second)
        inferred_eta = arrival_time(first, second)
        assert math.isclose(inferred_eta, x / speed, rel_tol=1e-11)
        if abs(climb) < 1e-12:
            assert math.isclose(single_lamp_arrival(a0, a1), inferred_eta, rel_tol=1e-11)
    assert math.isclose(single_lamp_arrival(304, 272), 19)
    assert math.isclose(arrival_time(Bird(152, 20), Bird(136, 20)), 19)
    print("SHADOW_CLOCK_SELFTEST_OK: 2000 moving birds + 19-second calibration")


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--selftest", action="store_true")
    parser.add_argument("--climb", type=float, default=0.0, help="bird climb rate in m/s for the demo")
    parser.add_argument("--noise", type=float, help="run 20,000 seeded trials with this shadow-centroid noise sigma, metres")
    args = parser.parse_args()
    if args.selftest:
        selftest()
    elif args.noise is not None:
        noise_demo(args.noise)
    else:
        demo(args.climb)


if __name__ == "__main__":
    main()