#!/usr/bin/env python3
"""
Reflected-image compression over water: globe vs topographic plane.

Solves the specular point exactly (Alhazen's problem) for an observer at height
h_o looking at a feature at height h, horizontal surface distance D, on:

  * a sphere of radius R  (optionally R_eff = R/(1-k) for refraction)
  * a plane               (the topographic flat model)

and reports the ratio

    reflected angular span / direct angular span

for a pair of features on the mountain. On a plane that ratio is ~1 at every
observer height. On a sphere the reflecting strip is convex, so the reflection
is demagnified and the ratio falls as the observer climbs.

No dependencies beyond the standard library.
"""

import math
from bisect import bisect

# ---------------------------------------------------------------- geometry ---

def _unit(v):
    n = math.hypot(v[0], v[1])
    return (v[0] / n, v[1] / n)


def _angle_between(a, b):
    ax, ay = _unit(a)
    bx, by = _unit(b)
    dot = max(-1.0, min(1.0, ax * bx + ay * by))
    return math.acos(dot)


def specular_sphere(R, h_o, h, D):
    """Angle (from observer, at Earth's centre) of the specular point.

    Observer sits at angle 0 and radius R+h_o; target at angle D/R and radius
    R+h. Returns theta in [0, D/R] where incidence == reflection about the
    local radial normal, plus the surface distance from the observer.
    """
    th_D = D / R
    O = (0.0, R + h_o)
    T = (math.sin(th_D) * (R + h), math.cos(th_D) * (R + h))

    def imbalance(th):
        P = (math.sin(th) * R, math.cos(th) * R)
        N = (math.sin(th), math.cos(th))
        to_O = (O[0] - P[0], O[1] - P[1])
        to_T = (T[0] - P[0], T[1] - P[1])
        return _angle_between(to_O, N) - _angle_between(to_T, N)

    lo, hi = 1e-15, th_D - 1e-15
    flo = imbalance(lo)
    for _ in range(200):
        mid = 0.5 * (lo + hi)
        fm = imbalance(mid)
        if (fm < 0) == (flo < 0):
            lo, flo = mid, fm
        else:
            hi = mid
    th = 0.5 * (lo + hi)
    return th, th * R


def reflected_dir_sphere(R, h_o, h, D):
    """Unit direction from observer toward the reflected image of the feature."""
    th, _ = specular_sphere(R, h_o, h, D)
    O = (0.0, R + h_o)
    P = (math.sin(th) * R, math.cos(th) * R)
    return (P[0] - O[0], P[1] - O[1])


def direct_dir_sphere(R, h_o, h, D):
    th_D = D / R
    O = (0.0, R + h_o)
    T = (math.sin(th_D) * (R + h), math.cos(th_D) * (R + h))
    return (T[0] - O[0], T[1] - O[1])


def specular_plane(h_o, h, D):
    """Distance from observer to the specular point on a plane."""
    return D * h_o / (h_o + h)


def reflected_dir_plane(h_o, h, D):
    # image of the feature sits at depth -h, same horizontal distance
    return (D, -h - h_o)


def direct_dir_plane(h_o, h, D):
    return (D, h - h_o)


# ------------------------------------------------------------------ ratios ---

def ratio_sphere(R, h_o, h_lo, h_hi, D):
    d1 = direct_dir_sphere(R, h_o, h_lo, D)
    d2 = direct_dir_sphere(R, h_o, h_hi, D)
    r1 = reflected_dir_sphere(R, h_o, h_lo, D)
    r2 = reflected_dir_sphere(R, h_o, h_hi, D)
    return _angle_between(r1, r2) / _angle_between(d1, d2)


def ratio_plane(h_o, h_lo, h_hi, D):
    d1 = direct_dir_plane(h_o, h_lo, D)
    d2 = direct_dir_plane(h_o, h_hi, D)
    r1 = reflected_dir_plane(h_o, h_lo, D)
    r2 = reflected_dir_plane(h_o, h_hi, D)
    return _angle_between(r1, r2) / _angle_between(d1, d2)


# ------------------------------------------------------------------- site ----

R_EARTH = 6371000.0
SITE = dict(D=60000.0, summit=3206.0)   # Aoraki above Lake Pukaki, 60 km


def r_eff(k):
    return R_EARTH / (1.0 - k)


def report(h_o, h_lo, h_hi, k=0.0, D=SITE["D"]):
    R = r_eff(k)
    rs = ratio_sphere(R, h_o, h_lo, h_hi, D)
    rp = ratio_plane(h_o, h_lo, h_hi, D)
    return rs, rp, (1 - rs) * 100.0, (1 - rp) * 100.0


if __name__ == "__main__":
    D = SITE["D"]
    summit = SITE["summit"]

    print("=== strip endpoints, observer 2 m (validating against the page) ===")
    for label, h in (("summit 3206 m", summit), ("1000 m contour", 1000.0)):
        _, ds = specular_sphere(R_EARTH, 2.0, h, D)
        dp = specular_plane(2.0, h, D)
        print(f"  {label:16s} sphere {ds:8.1f} m   plane {dp:8.1f} m")

    print("\n=== compression, 1000 m -> summit band ===")
    for h_o in (2.0, 120.0):
        rs, rp, cs, cp = report(h_o, 1000.0, summit)
        print(f"  h_o={h_o:6.1f} m   sphere {cs:6.3f}%   plane {cp:6.4f}%")

    print("\n=== compression, high pair 2000 m -> summit ===")
    for h_o in (2.0, 120.0):
        rs, rp, cs, cp = report(h_o, 2000.0, summit)
        print(f"  h_o={h_o:6.1f} m   sphere {cs:6.3f}%   plane {cp:6.4f}%")

    print("\n=== refraction sensitivity at 120 m, 1000 m -> summit ===")
    for k in (0.0, 0.13, 0.17, 0.30):
        rs, rp, cs, cp = report(120.0, 1000.0, summit, k=k)
        print(f"  k={k:4.2f}   sphere {cs:6.3f}%")


# ------------------------------------------------------------------ figure ---

def sweep(h_lo, h_hi, k=0.0, hmax=150.0, n=140, D=SITE["D"]):
    R = r_eff(k)
    xs, ys = [], []
    for i in range(n + 1):
        h_o = 1.0 + (hmax - 1.0) * i / n
        xs.append(h_o)
        ys.append((1 - ratio_sphere(R, h_o, h_lo, h_hi, D)) * 100.0)
    return xs, ys


def sweep_plane(h_lo, h_hi, hmax=150.0, n=140, D=SITE["D"]):
    xs, ys = [], []
    for i in range(n + 1):
        h_o = 1.0 + (hmax - 1.0) * i / n
        xs.append(h_o)
        ys.append((1 - ratio_plane(h_o, h_lo, h_hi, D)) * 100.0)
    return xs, ys


def make_figure(path="docs/mirrored-reflections/img/fig-compression-vs-height.svg"):
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    S = SITE["summit"]
    ACCENT, UMBRA, INK = "#4A6FA5", "#c0392b", "#6E675B"

    fig, ax = plt.subplots(figsize=(7.6, 4.5))

    x, y0 = sweep(2500.0, S, k=0.00)
    _, y17 = sweep(2500.0, S, k=0.17)
    ax.fill_between(x, y17, y0, color=UMBRA, alpha=0.16, lw=0)
    ax.plot(x, y0, color=UMBRA, lw=2.2,
            label="Globe — high pair (2,500 m to summit), no refraction")
    ax.plot(x, y17, color=UMBRA, lw=1.1, ls=":",
            label="same pair, standard refraction (k = 0.17)")

    _, yb = sweep(1000.0, S, k=0.00)
    ax.plot(x, yb, color=UMBRA, lw=1.1, ls="--", alpha=.75,
            label="Globe — 1,000 m to summit (includes unobservable base)")

    xp, yp = sweep_plane(2500.0, S)
    ax.plot(xp, yp, color=ACCENT, lw=2.2, label="Topographic plane — same pair")

    ax.axvline(120, color=INK, lw=.8, ls=":")
    ax.annotate("120 m\n(proposed drone ceiling)", xy=(120, 2.4),
                xytext=(122, 2.35), fontsize=8, color=INK, va="top")

    ax.set_xlabel("Observer height above the water (m)")
    ax.set_ylabel("Compression of the reflected image (%)")
    ax.set_xlim(0, 150)
    ax.set_ylim(-0.2, 5.6)
    ax.grid(alpha=.18, lw=.6)
    for s in ("top", "right"):
        ax.spines[s].set_visible(False)
    ax.legend(fontsize=7.6, loc="upper left", frameon=False)
    ax.set_title("Aoraki reflected in Lake Pukaki: what each model predicts\n"
                 "peak 3,206 m above the water at 60 km", fontsize=9.5,
                 color="#2E4057", loc="left")
    fig.tight_layout()
    fig.savefig(path, format="svg", transparent=True)
    print("wrote", path)


def make_tilt_figure(path="docs/mirrored-reflections/img/fig-tilt-budget.svg", h_o=2.0):
    """Horizontal log-scale comparison of the tilt budget that makes a glitter path."""
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    AM = 180 * 60 / math.pi
    dip = math.sqrt(2 * h_o / R_EARTH)
    cm = lambda W: math.sqrt(0.003 + 0.00512 * W) * AM

    rows = [
        ("What curvature contributes\n(half the horizon dip, 2 m eye)", dip / 2 * AM, "#c0392b"),
        ("Gentlest waves that can form a streak\n(about 0.2°)", 12.0, "#6E675B"),
        ("Cox–Munk dead calm\n(RMS slope, open ocean)", cm(0), "#4A6FA5"),
        ("Light breeze, 3 m/s\n(RMS slope)", cm(3), "#4A6FA5"),
    ]
    labels = [r[0] for r in rows][::-1]
    vals = [r[1] for r in rows][::-1]
    cols = [r[2] for r in rows][::-1]

    fig, ax = plt.subplots(figsize=(7.8, 3.3))
    ax.barh(range(len(vals)), vals, color=cols, height=.6)
    for i, v in enumerate(vals):
        ax.text(v * 1.12, i, f"{v:,.1f}'", va="center", fontsize=8.5,
                color="#3A352D")
    ax.set_yticks(range(len(labels)))
    ax.set_yticklabels(labels, fontsize=8)
    ax.set_xscale("log")
    ax.set_xlim(0.8, 3000)
    ax.set_xlabel("Surface tilt available or required (arcminutes, log scale)")
    ax.grid(axis="x", alpha=.18, lw=.6)
    for s in ("top", "right", "left"):
        ax.spines[s].set_visible(False)
    ax.set_title("Why curvature cannot break a glitter path\n"
                 "the tilt budget at an eye 2 m above the water",
                 fontsize=9.5, color="#2E4057", loc="left")
    fig.tight_layout()
    fig.savefig(path, format="svg", transparent=True)
    print("wrote", path)


# --------------------------------------------------- synthetic scene render ---

AM = 180 * 60 / math.pi   # radians -> arcminutes


def _ridge(n=420, base=1000.0, peak=None, halfwidth=1800.0):
    """A synthetic peak profile: horizontal offset (m) against height (m).

    Deliberately synthetic. We define the profile, so the warp applied below is
    computed from the solver rather than illustrated over a photograph whose
    true feature elevations we would not know.
    """
    peak = peak or SITE["summit"]
    xs, hs = [], []
    for i in range(n + 1):
        u = -1.0 + 2.0 * i / n
        x = u * halfwidth * 1.9
        core = math.exp(-(x / halfwidth) ** 2)
        shoulder = 0.34 * math.exp(-((x - halfwidth * 1.15) / (halfwidth * .85)) ** 2)
        shoulder += 0.27 * math.exp(-((x + halfwidth * 1.25) / (halfwidth * .8)) ** 2)
        notch = 0.05 * math.exp(-((x + halfwidth * .45) / (halfwidth * .16)) ** 2)
        h = base + (peak - base) * min(1.0, core + shoulder - notch)
        xs.append(x); hs.append(max(base, h))
    return xs, hs


def scene(h_o, D=SITE["D"], k=0.0):
    """Both complete scenes, registered on the direct image of the summit.

    A photographer cannot align two frames on an absolute horizon: the bulk
    displacement between models is not observable. What is observable is the
    reflection's position relative to identifiable DIRECT features. So each
    model's scene is shifted so its direct summit sits at the same place, and
    the residual in the reflection is what a measurement would actually see.
    """
    R = r_eff(k)
    xs, hs = _ridge()
    top = max(range(len(hs)), key=lambda i: hs[i])

    def one(direct_fn, refl_fn):
        d = [math.atan2(*reversed(direct_fn(h))) * AM for h in hs]
        r = [math.atan2(*reversed(refl_fn(h))) * AM for h in hs]
        shift = -d[top]
        return [v + shift for v in d], [v + shift for v in r]

    dp, rp = one(lambda h: direct_dir_plane(h_o, h, D),
                 lambda h: reflected_dir_plane(h_o, h, D))
    ds, rs = one(lambda h: direct_dir_sphere(R, h_o, h, D),
                 lambda h: reflected_dir_sphere(R, h_o, h, D))
    ang = [math.atan2(x, D) * AM for x in xs]
    return ang, dp, rp, ds, rs


def make_scene_figure(path="docs/mirrored-reflections/img/fig-what-you-would-see.svg",
                      h_o=120.0, exag=25.0, px_tall=4000, fov_deg=16.0):
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    ang, dp, rp, ds, rs = scene(h_o)
    px_per_am = px_tall / (fov_deg * 60.0)
    worst = max(abs(a - b) for a, b in zip(rp, rs))
    worst_direct = max(abs(a - b) for a, b in zip(dp, ds))

    ROCK, PLN, SPH, SKY = "#3A352D", "#4A6FA5", "#c0392b", "#eef1f4"
    fig, axes = plt.subplots(1, 2, figsize=(8.8, 5.2), sharey=True)

    lo = min(rp) - 10
    for panel, E, is_true in zip(axes, (1.0, exag), (True, False)):
        panel.fill_between(ang, dp, max(dp) + 60, color=SKY, lw=0)
        panel.fill_between(ang, lo - 60, dp, color=ROCK, lw=0)
        warp = [p + E * (s - p) for p, s in zip(rp, rs)]
        panel.fill_between(ang, rp, lo - 60, color=PLN, alpha=.30, lw=0)
        panel.plot(ang, rp, color=PLN, lw=1.7, label="topographic plane")
        panel.plot(ang, warp, color=SPH, lw=1.7, label="globe")
        panel.set_xlim(min(ang), max(ang))
        panel.set_xlabel("angle across the frame (arcmin)", fontsize=8)
        panel.tick_params(labelsize=7.5)
        for sp in ("top", "right"):
            panel.spines[sp].set_visible(False)

    axes[0].set_ylim(lo, max(dp) + 16)
    axes[0].set_ylabel("angle relative to the direct summit (arcmin)", fontsize=8)
    axes[0].set_title("True scale — drawn exactly", fontsize=9.6,
                      color="#2E4057", loc="left")
    axes[1].set_title(f"The same difference, multiplied {exag:.0f}×",
                      fontsize=9.6, color=SPH, loc="left")
    axes[0].legend(fontsize=7.4, loc="lower center", frameon=False,
                   ncol=2, bbox_to_anchor=(.5, -.02))
    axes[0].annotate("both are drawn.\nthe globe is under the plane.",
                     xy=(0, max(dp) - 6), fontsize=7.4, color="#6E675B",
                     ha="center", va="top")
    axes[1].annotate("EXAGGERATED", xy=(0, max(dp) - 6), fontsize=8.4,
                     color=SPH, ha="center", va="top", weight="bold")

    fig.suptitle(
        "What a globe would actually do to the reflection\n"
        f"Aoraki from {h_o:.0f} m over Lake Pukaki, registered on the direct summit — "
        f"largest mismatch {worst:.2f}\u2032, about {worst*px_per_am:.0f} px "
        f"on a {px_tall}-px frame at a {fov_deg:.0f}\u00b0 vertical field",
        fontsize=9.6, color="#2E4057", x=.012, ha="left")
    fig.tight_layout(rect=[0, 0, 1, .91])
    fig.savefig(path, format="svg", transparent=True)
    print(f"worst reflection mismatch {worst:.3f}' = {worst*px_per_am:.1f} px;"
          f"  residual on the direct image {worst_direct:.3f}'")


# ----------------------------------------------- the angle test, as figures ---

def gap_arcmin(h, D, k=0.0):
    """Summit-to-reflected-summit angle. Camera height cancels out of it."""
    lift = D * D * k / (2 * R_EARTH)
    return 2 * math.atan((h + lift) / D) * AM


def gap_globe(h, D, k=0.0, h_o=2.0):
    """Exact: solves the specular point rather than approximating the drop.

    The approximation 2*arctan((h - D^2/2R)/D) is good to about 0.13 arcmin
    here, which is enough to make a figure disagree with the page's own text.
    """
    R = r_eff(k)
    a = lambda v: math.atan2(v[1], v[0]) * AM
    return a(direct_dir_sphere(R, h_o, h, D)) - a(reflected_dir_sphere(R, h_o, h, D))


def make_angle_figure(path="docs/mirrored-reflections/img/fig-the-angle.svg",
                      D=60000.0, h=3206.0):
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    f, g = gap_arcmin(h, D), gap_globe(h, D)
    drop = D * D / (2 * R_EARTH)
    PLN, SPH, INK, WATER = "#4A6FA5", "#c0392b", "#3A352D", "#cfe0ee"
    fig, ax = plt.subplots(figsize=(8.4, 4.0))

    # schematic: vertical exaggerated ~6x so the rays are legible
    EX = 6.0
    Dk = D / 1000.0
    ax.axhspan(-h * EX / 1000 * 1.25, 0, color=WATER, alpha=.55, lw=0)
    ax.axhline(0, color="#7d8fa0", lw=1.2)

    for lab, hh, col, ls in (("plane: summit at 3,206 m", h, PLN, "-"),
                             (f"globe: curvature drops it {drop:.0f} m", h - drop, SPH, "-")):
        y = hh * EX / 1000
        ax.plot([Dk, Dk], [0, y], color=col, lw=2.4, solid_capstyle="butt")
        ax.plot([0, Dk], [0, y], color=col, lw=1.0, ls=ls, alpha=.85)      # direct ray
        ax.plot([0, Dk], [0, -y], color=col, lw=1.0, ls=":", alpha=.85)    # reflected
        ax.plot([Dk, Dk], [0, -y], color=col, lw=2.4, alpha=.35,
                solid_capstyle="butt")
        ax.plot([Dk], [y], marker="^", color=col, ms=7, label=lab)

    ax.plot([0], [0], marker="o", color=INK, ms=6)
    ax.annotate("camera, at the water\n(its height cancels)", xy=(0, 0),
                xytext=(Dk * .05, h * EX / 1000 * .42), fontsize=8, color=INK)
    ax.annotate(f"plane predicts  {f:.2f}′\nglobe predicts  {g:.2f}′\n"
                f"difference  {f-g:.2f}′",
                xy=(Dk * .04, -h * EX / 1000 * .80), fontsize=9.2, color=INK,
                bbox=dict(boxstyle="round,pad=.45", fc="#fcfbf8", ec="#e4ddcb"))
    ax.text(Dk * .5, h * EX / 1000 * .06,
            "solid = direct view   dotted = reflected view", fontsize=7.6,
            color="#6E675B", ha="center")

    ax.set_xlim(-Dk * .03, Dk * 1.06)
    ax.set_ylim(-h * EX / 1000 * 1.25, h * EX / 1000 * 1.25)
    ax.set_xlabel(f"distance from the camera (km) — vertical scale "
                  f"exaggerated {EX:.0f}× so the rays are visible", fontsize=8)
    ax.set_yticks([])
    for s in ("top", "right", "left"):
        ax.spines[s].set_visible(False)
    ax.legend(fontsize=8, loc="upper left", frameon=False)
    ax.set_title("The angle from a summit to its own reflection\n"
                 f"Aoraki, 3,206 m above the water at {Dk:.0f} km",
                 fontsize=9.6, color="#2E4057", loc="left")
    fig.tight_layout()
    fig.savefig(path, format="svg", transparent=True)
    print("wrote", path, f"  flat {f:.2f}' globe {g:.2f}' sep {f-g:.2f}'")


def make_refraction_figure(path="docs/mirrored-reflections/img/fig-refraction-cancels.svg",
                           D=60000.0, h=3206.0):
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    ks = [i / 200.0 for i in range(0, 61)]        # k = 0 .. 0.30
    F = [gap_arcmin(h, D, k) for k in ks]
    G = [gap_globe(h, D, k) for k in ks]
    PLN, SPH, INK = "#4A6FA5", "#c0392b", "#6E675B"

    fig, ax = plt.subplots(figsize=(7.6, 4.2))
    ax.fill_between(ks, G, F, color="#d9a441", alpha=.20, lw=0)
    ax.plot(ks, F, color=PLN, lw=2.4, label="what a plane predicts")
    ax.plot(ks, G, color=SPH, lw=2.4, label="what a globe predicts")
    for k in (0.0, 0.17, 0.30):
        f, g = gap_arcmin(h, D, k), gap_globe(h, D, k)
        ax.plot([k, k], [g, f], color=INK, lw=.8, ls=":")
        ax.annotate(f"{f-g:.2f}′", xy=(k, (f + g) / 2), fontsize=8.4,
                    color="#8a6d1f", ha="center",
                    bbox=dict(boxstyle="round,pad=.28", fc="#fffdf7", ec="none"))
    ax.set_xlabel("refraction coefficient k  (0 = straight rays, 0.17 = standard, "
                  "0.30 = strong inversion)", fontsize=8)
    ax.set_ylabel("summit-to-reflection angle (arcmin)", fontsize=8.4)
    ax.grid(alpha=.18, lw=.6)
    for s in ("top", "right"):
        ax.spines[s].set_visible(False)
    ax.legend(fontsize=8.4, loc="center left", frameon=False)
    ax.set_title("Refraction moves both predictions, and not the gap between them\n"
                 "both rise about 10′ across this range; the separation moves "
                 "by six hundredths of one",
                 fontsize=9.5, color="#2E4057", loc="left")
    fig.tight_layout()
    fig.savefig(path, format="svg", transparent=True)
    print("wrote", path)