# SPDX-License-Identifier: MIT
"""Build the figures for YS-EN-0002 from results.json and write them into the
note between <!-- fig:name --> markers.

Usage: python build_figures.py   (after parallax_fix.py)
"""
from __future__ import annotations

import json
import math
import sys
from pathlib import Path

HERE = Path(__file__).resolve().parent
sys.path.insert(0, str(HERE.parent))
import svgchart as sc  # noqa: E402

PAGE = HERE.parent.parent / "ys-en-0002-navigating-by-parallax.html"
R = json.loads((HERE / "results.json").read_text())
AU_KM = 149_597_870.7

# colour follows the star in every figure
STAR = {"Proxima Centauri": 1, "Barnard's Star": 2, "Wolf 359": 3}


def arcsec(v):
    return f"{v:.3g}″" if v < 100 else f"{v:,.0f}″"


# ---------------------------------------------------------------- figure 1
def fig_shift():
    s = R["shift"]
    xs = s["baseline_au"]
    series = [{"slot": STAR[n], "name": n, "values": s[n]} for n in STAR]
    fr = sc.Frame(720, 360, "Maximum parallax shift of Proxima Centauri, Barnard's Star and Wolf 359 versus baseline")
    fr.x_scale(1, 1000, log=True)
    fr.y_scale(0.3, 1000, log=True)
    fr.grid_y([1, 10, 100, 1000], fmt=lambda v: f"{v:,.0f}″")
    fr.grid_x([1, 10, 100, 1000], fmt=lambda v: f"{v:,.0f} au", lines=False)
    fr.baseline()
    for x, lab, dy in ((1, "Earth's orbit", 0), (47.1176, "New Horizons, 2020", 0), (171.93, "Voyager 1, 2026", 14)):
        xx = fr.x(x)
        fr.out.append(f'<line class="grid" x1="{xx:.1f}" x2="{xx:.1f}" y1="{fr.mt}" y2="{fr.mt + fr.ph}"/>')
        fr.note(xx + 5, fr.mt + 11 + dy, lab)
    fr.titles("Baseline perpendicular to the line of sight (log)", "Apparent shift (log)")
    fr.cursor()
    sc.draw_lines(fr, xs, series)
    # the measured shifts, Earth–New Horizons separation 46.85 au
    for name, v in (("Proxima Centauri", 32.36), ("Wolf 359", 15.72)):
        cx, cy = fr.x(46.85), fr.y(v)
        fr.out.append(f'<rect class="s{STAR[name]}" x="{cx - 5:.1f}" y="{cy - 5:.1f}" width="10" height="10" '
                      f'transform="rotate(45 {cx:.1f} {cy:.1f})" stroke="var(--panel)" stroke-width="2"/>')
    fr.note(fr.x(46.85) - 9, fr.y(32.36) - 9, "measured, 2020", anchor="end")
    sc.column_hits(fr, xs, series, lambda i, x: f"Baseline {x:,.3g} au", arcsec)
    rows = [[f"{x:,.4g} au"] + [arcsec(s[n][i]) for n in STAR] for i, x in enumerate(xs) if i % 5 == 0]
    rows.append(["46.85 au, measured 2020", "32.36″", "—", "15.72″"])
    return sc.figure(
        code="FIG. 1 — PARALLAX SHIFT VS BASELINE",
        fid="fig-shift",
        plot=fr.svg(),
        legend_html=sc.legend(series) .replace("</div>", '<span class="k0"><i style="width:9px;height:9px;transform:rotate(45deg);background:var(--mist)"></i>Measured by New Horizons (Lauer et al.)</span></div>'),
        caption=("<b>Figure 1.</b> The largest possible shift, when the baseline is perpendicular to the line of "
                 "sight, from Gaia DR3 distances. The diamonds are the two shifts measured between Earth and New "
                 "Horizons; they sit below the lines because the real baseline was at an angle to each star."),
        data_html=sc.table("Maximum shift by baseline", ["Baseline", *STAR], rows),
    )


# ---------------------------------------------------------------- figure 2 (diagram)
def fig_lines():
    W, H = 720, 360
    sun = (120, 300)
    craft = (430, 250)
    stars = {"Proxima Centauri": (60, 40, "1.30 pc"), "Wolf 359": (700, 30, "2.41 pc")}
    o = [f'<svg class="dg" viewBox="0 0 {W} {H}" role="img" aria-label="Two lines of position from the stars cross at the spacecraft; each line can slide sideways by distance times angular error">']

    def band(sx, sy, half_w, cls):
        # a wedge from the star widening toward and past the craft
        dx, dy = craft[0] - sx, craft[1] - sy
        L = math.hypot(dx, dy)
        ux, uy = dx / L, dy / L
        nx, ny = -uy, ux
        ex, ey = sx + ux * (L + 90), sy + uy * (L + 90)
        k = (L + 90) / L
        pts = [(sx, sy), (ex + nx * half_w * k, ey + ny * half_w * k), (ex - nx * half_w * k, ey - ny * half_w * k)]
        o.append(f'<path class="{cls}" opacity="0.16" d="M{pts[0][0]},{pts[0][1]} L{pts[1][0]:.1f},{pts[1][1]:.1f} L{pts[2][0]:.1f},{pts[2][1]:.1f} Z"/>')
        o.append(f'<line class="msg {cls.replace("f", "s")}" x1="{sx}" y1="{sy}" x2="{ex:.1f}" y2="{ey:.1f}"/>')

    band(*stars["Proxima Centauri"][:2], 14, "f1")
    band(*stars["Wolf 359"][:2], 26, "f3")
    for name, (x, y, dist) in stars.items():
        o.append(f'<circle class="f{STAR[name]}" cx="{x}" cy="{y}" r="6"/>')
        anchor = "start" if x < W / 2 else "end"
        dx = 14 if x < W / 2 else -14
        o.append(f'<text class="lbl" x="{x + dx}" y="{y + 4}" text-anchor="{anchor}">{name}</text>')
        o.append(f'<text x="{x + dx}" y="{y + 19}" text-anchor="{anchor}">{dist} · not to scale</text>')
    o.append(f'<circle class="f1" cx="{sun[0]}" cy="{sun[1]}" r="9"/>')
    o.append(f'<text class="lbl" x="{sun[0]}" y="{sun[1] + 28}" text-anchor="middle">Sun</text>')
    o.append(f'<line class="rule" x1="{sun[0]}" y1="{sun[1]}" x2="{craft[0]}" y2="{craft[1]}"/>')
    o.append(f'<text x="{(sun[0] + craft[0]) / 2}" y="{(sun[1] + craft[1]) / 2 + 22}" text-anchor="middle">47.1 au</text>')
    o.append(f'<rect class="f0" x="{craft[0] - 5}" y="{craft[1] - 5}" width="10" height="10" transform="rotate(45 {craft[0]} {craft[1]})"/>')
    o.append(f'<text class="lbl" x="{craft[0] + 22}" y="{craft[1] + 38}">New Horizons</text>')
    o.append(f'<text x="{craft[0] + 22}" y="{craft[1] + 53}">where the lines cross</text>')
    o.append(f'<text x="250" y="150" text-anchor="middle">line of position</text>')
    o.append(f'<text x="600" y="190" text-anchor="middle">band width = d · σ</text>')
    o.append(f'<text x="600" y="205" text-anchor="middle">(wider: Wolf 359 is farther)</text>')
    o.append("</svg>")
    return sc.figure(
        code="FIG. 2 — LINES OF POSITION",
        fid="fig-lines",
        plot="\n".join(o),
        caption=("<b>Figure 2.</b> Each measured direction, drawn back through the star's catalogue position, is "
                 "a line the spacecraft must lie on. Pointing error σ turns each line into a band whose width at "
                 "the spacecraft is the star's distance times σ, so the farther star gives the wider band. "
                 "Schematic, not to scale: the stars are thousands of times farther away than drawn."),
    )


# ---------------------------------------------------------------- figure 3
def fig_scatter():
    mc = R["monte_carlo"]
    a = R["nh_fix"]["axes_au"]
    fx, fy = mc["lorri_fix_in_basis"][:2]
    span_x, span_y = 1.5, 0.85
    pw = 720 - 58 - 132
    ph = round(pw * (2 * span_y) / (2 * span_x))  # equal scale on both axes
    fr = sc.Frame(720, ph + 72, "Monte Carlo cloud of 4000 simulated fixes, with the actual LORRI fix and the true position",
                  margins=(58, 132, 26, 46))
    fr.x_scale(-span_x, span_x)
    fr.y_scale(-span_y, span_y)
    fr.grid_y([-0.5, 0, 0.5], fmt=lambda v: f"{v:+.1f}" if v else "0")
    fr.grid_x([-1, -0.5, 0, 0.5, 1], fmt=lambda v: f"{v:+.1f}" if v else "0")
    fr.titles("Along the ellipsoid's long axis, roughly toward Proxima (au)", "Middle axis (au)")
    for x, y in mc["sample"]:
        fr.out.append(f'<circle class="s2" cx="{fr.x(x):.1f}" cy="{fr.y(y):.1f}" r="1.6" stroke="none" opacity="0.55"/>')
    for k, op in ((1, 0.9), (2, 0.5)):
        rx = (fr.x(a[0] * k) - fr.x(0))
        ry = (fr.y(0) - fr.y(a[1] * k))
        fr.out.append(f'<ellipse cx="{fr.x(0):.1f}" cy="{fr.y(0):.1f}" rx="{rx:.1f}" ry="{ry:.1f}" fill="none" '
                      f'stroke="var(--ivory)" stroke-opacity="{op}" stroke-width="1"/>')
        fr.note(fr.x(0) + rx * 0.72 + 4, fr.y(0) - ry * 0.72 - 4, f"{k}σ")
    tx, ty = fr.x(0), fr.y(0)
    fr.out.append(f'<circle class="dot s3" cx="{tx:.1f}" cy="{ty:.1f}" r="5"/>')
    fr.note(tx + 10, ty + 16, "JPL / DSN (truth)", cls="ink")
    px, py = fr.x(fx), fr.y(fy)
    fr.out.append(f'<circle class="dot s1" cx="{px:.1f}" cy="{py:.1f}" r="5.5"/>')
    fr.note(px + 10, py - 8, "LORRI fix", cls="ink")
    e = R["nh_fix"]["error_norm_au"]
    sc.point_hit(fr, tx, ty, 0, "True position (JPL)", "3|Offset|0 au")
    sc.point_hit(fr, px, py, 1, "Fix from LORRI images",
                 f"1|Long axis|{fx:+.3f} au;1|Middle axis|{fy:+.3f} au;1|Total error (3D)|{e:.3f} au")
    std = mc["std_along_axes_au"]
    rows = [["Long axis", f"{a[0]:.3f} au", f"{std[0]:.3f} au", f"{fx:+.3f} au"],
            ["Middle axis", f"{a[1]:.3f} au", f"{std[1]:.3f} au", f"{fy:+.3f} au"],
            ["Short axis", f"{a[2]:.3f} au", f"{std[2]:.3f} au", f"{mc['lorri_fix_in_basis'][2]:+.3f} au"]]
    leg = ('<div class="legend"><span class="k2"><i></i>Simulated fixes (σ = 0.18″)</span>'
           '<span class="k1"><i></i>Actual LORRI fix</span><span class="k3"><i></i>True position</span></div>')
    return sc.figure(
        code="FIG. 3 — 4,000 FIXES THAT COULD HAVE BEEN",
        fid="fig-scatter",
        plot=fr.svg(),
        legend_html=leg,
        caption=(f"<b>Figure 3.</b> We re-observed the true directions {mc['n']:,} times with 0.18″ of random "
                 "error per star and solved each time; 700 of the fixes are plotted, projected onto the "
                 "ellipsoid's two longest axes (equal scale). Their spread matches the analytic covariance, "
                 f"and the real fix sits inside the 1σ ellipse in this view ({R['nh_fix']['error_sigma']:.2f}σ in 3D)."),
        data_html=sc.table("Error ellipsoid: analytic vs Monte Carlo",
                           ["Axis", "1σ (analytic)", "1σ (Monte Carlo)", "LORRI fix offset"], rows),
    )


# ---------------------------------------------------------------- figure 4
def fig_dop():
    d = R["dop"]
    xs = d["sep_deg"]
    series = [
        {"slot": 2, "name": "2nd star at 1.83 pc (Barnard's)", "end": "at 1.83 pc", "values": d["with_barnard_distance"]},
        {"slot": 3, "name": "2nd star at 2.41 pc (Wolf 359's)", "end": "at 2.41 pc", "values": d["with_wolf_distance"]},
    ]
    fr = sc.Frame(720, 330, "Long axis of the error ellipsoid versus angle between the two stars")
    fr.x_scale(0, 180)
    fr.y_scale(1, 30, log=True)
    fr.grid_y([1, 2, 3, 5, 10, 20, 30], fmt=lambda v: f"{v:g}")
    fr.grid_x([0, 30, 60, 90, 120, 150, 180], fmt=lambda v: f"{v:.0f}°", lines=False)
    fr.baseline()
    fr.titles("Angle between the two stars, seen from the Sun", "Long semi-axis, au per arcsecond (log)")
    fr.cursor()
    sc.draw_lines(fr, xs, series)
    pairs = {(p["a"], p["b"]): p for p in R["pairs"]}
    for key, slot, lab, dy in ((("Proxima Centauri", "Barnard's Star"), 2, "Proxima + Barnard's", 16),
                               (("Proxima Centauri", "Wolf 359"), 3, "Proxima + Wolf 359 (used)", -10)):
        p = pairs[key]
        cx, cy = fr.x(p["sep_deg"]), fr.y(p["long"])
        fr.out.append(f'<circle class="dot s{slot}" cx="{cx:.1f}" cy="{cy:.1f}" r="5"/>')
        fr.note(cx, cy + dy, f'{lab}: {p["long"]:.2f}', anchor="middle", cls="ink")
    sc.column_hits(fr, xs, series, lambda i, x: f"Stars {x}° apart", lambda v: f"{v:.2f} au/″")
    rows = [[f"{x}°", f"{a:.2f}", f"{b:.2f}"] for x, a, b in zip(xs, d["with_barnard_distance"], d["with_wolf_distance"]) if x % 15 == 0]
    return sc.figure(
        code="FIG. 4 — GEOMETRY: WHY THE ANGLE MATTERS",
        fid="fig-dop",
        plot=fr.svg(),
        legend_html=sc.legend(series),
        caption=("<b>Figure 4.</b> Proxima Centauri plus a second star at a fixed distance, swept around the sky. "
                 "The error blows up as the two lines of sight become parallel (0° or 180°). It is smallest "
                 "at 90°, where the long axis equals the farther star's distance in parsecs. The two real pairs are marked at "
                 "their Gaia DR3 separations."),
        data_html=sc.table("Long semi-axis (au per arcsec) by separation", ["Separation", "2nd at 1.83 pc", "2nd at 2.41 pc"], rows),
    )


# ---------------------------------------------------------------- figure 5 (table)
def fig_pairs():
    paper = {("Proxima Centauri", "Barnard's Star"): "1.91 / 1.06", ("Proxima Centauri", "Wolf 359"): "2.45 / 1.15",
             ("Proxima Centauri", "Lalande 21185"): "2.70 / 1.16", ("Proxima Centauri", "Ross 154"): "3.64 / 1.19",
             ("Barnard's Star", "Wolf 359"): "2.58 / 1.46", ("Barnard's Star", "Lalande 21185"): "2.60 / 1.49",
             ("Barnard's Star", "Ross 154"): "6.54 / 1.56", ("Wolf 359", "Lalande 21185"): "7.00 / 1.75",
             ("Wolf 359", "Ross 154"): "3.80 / 1.87", ("Lalande 21185", "Ross 154"): "4.24 / 1.93"}
    rows = []
    for p in sorted(R["pairs"], key=lambda p: p["long"]):
        rows.append([f'{p["a"]} + {p["b"]}', f'{p["sep_deg"]:.1f}°', f'{p["long"]:.2f} / {p["short"]:.2f}', paper[(p["a"], p["b"])]])
    tbl = sc.table("Error semi-axes for every pair, au per arcsecond of pointing error",
                   ["Pair", "Separation", "Long / short (this note)", "Lauer et al. Table 4"], rows)
    return (f'<figure class="fig" id="fig-pairs"><div class="media__bar"><span class="lamp"></span><span class="lamp">'
            f'</span><span class="lamp"></span><span>FIG. 5 — STAR PAIRS, RANKED</span></div>'
            f'<div class="fig__plot">{tbl}</div><figcaption><b>Figure 5.</b> Every pair of the five candidate stars within '
            f"3 pc (Lauer et al. Table 4), best first. Multiply by the pointing error in arcseconds to get the error "
            f"in au. Lalande 21185 is HD 95735 and Ross 154 is CD−23 14742 in the paper.</figcaption></figure>")


# ---------------------------------------------------------------- figure 6
def fig_ladder():
    y = R["ysnav"]
    items = [  # label, km, kind (1 measured, 2 projected/simulated, 3 derived/specified), note
        ("YS·NAV 7 course hold: ±0.02 ly", y["fix_au"] * AU_KM, 3, "our specification, whole crossing"),
        ("1″ camera, best pair", y["fix_at_1arcsec_au"] * AU_KM, 3, "this note, pair solve"),
        ("LORRI 2020, 1σ long axis", 0.44 * AU_KM, 1, "Lauer et al., formal 1σ"),
        ("LORRI 2020, actual error", 0.351 * AU_KM, 1, "Lauer et al."),
        ("LORRI, 4× images + Barnard's", 0.17 * AU_KM, 2, "Lauer et al., projected"),
        ("Few-mas camera, nearest pair", 0.01 * AU_KM, 2, "Lauer et al., projected"),
        ("YS·NAV 7 live fix: ±0.01 au", 0.01 * AU_KM, 3, "our specification"),
        ("20 stars at 1 mas, light-years out", 0.003 * AU_KM, 2, "Bailer-Jones, simulated"),
        ("DSN ΔDOR, plane of sky, 47 au", 0.281 * 47.1, 3, "DSN handbook, 281 m at 1 au, scaled"),
        ("X-ray pulsars, low Earth orbit", 10, 1, "Mitchell et al., < 10 km"),
        ("DSN ranging, line of sight", 0.001, 3, "≈ 1 m, handbook figure"),
    ]
    row_h = 30
    fr = sc.Frame(720, 50 + row_h * len(items) + 40, "Position accuracy of navigation methods, log scale from 1 metre to 0.02 light-years",
                  margins=(232, 24, 30, 46))
    fr.x_scale(1e-3, 1e12, log=True)
    fr.y_scale(0, 1)
    ticks = [1e-3, 1, 1e3, 1e6, 1e9, 1e12]
    names = {1e-3: "1 m", 1: "1 km", 1e3: "1,000 km", 1e6: "10⁶ km", 1e9: "10⁹ km", 1e12: "10¹² km"}
    fr.out.append('<g class="grid">')
    for t in ticks:
        fr.out.append(f'<line x1="{fr.x(t):.1f}" x2="{fr.x(t):.1f}" y1="{fr.mt}" y2="{fr.mt + fr.ph}"/>')
    fr.out.append("</g>")
    for t in ticks:
        fr.out.append(f'<text x="{fr.x(t):.1f}" y="{fr.mt + fr.ph + 18}" text-anchor="middle">{names[t]}</text>')
    au = fr.x(AU_KM)
    fr.out.append(f'<line class="axis" x1="{au:.1f}" x2="{au:.1f}" y1="{fr.mt}" y2="{fr.mt + fr.ph}"/>')
    fr.note(au + 5, fr.mt - 8, "1 au")
    fr.titles("Position uncertainty (log scale)")
    for i, (lab, km, kind, note) in enumerate(items):
        cy = fr.mt + row_h * (i + 0.5) + 6
        fr.out.append(f'<text class="{"ink" if kind else ""}" x="{fr.ml - 12}" y="{cy + 4:.1f}" text-anchor="end">{sc.esc(lab)}</text>')
        fr.out.append(f'<line class="grid" x1="{fr.ml}" x2="{fr.ml + fr.pw}" y1="{cy:.1f}" y2="{cy:.1f}"/>')
        cx = fr.x(km)
        fr.out.append(f'<circle class="dot s{kind}" cx="{cx:.1f}" cy="{cy:.1f}" r="5"/>')
        val = f"{km / AU_KM:.3g} au" if km > 1e5 else (f"{km * 1000:.0f} m" if km < 1 else f"{km:,.0f} km")
        sc.point_hit(fr, cx, cy, i, lab, f"{kind}|{note}|{val}", r=13)
    rows = []
    for lab, km, kind, note in items:
        val = f"{km / AU_KM:.3g} au" if km > 1e5 else (f"{km * 1000:.0f} m" if km < 1 else f"{km:,.0f} km")
        rows.append([lab, val, f"{km:.3g} km", {1: "Measured", 2: "Projected / simulated", 3: "Derived / specified"}[kind], note])
    leg = ('<div class="legend"><span class="k1"><i></i>Measured</span><span class="k2"><i></i>Projected or simulated</span>'
           '<span class="k3"><i></i>Derived or specified</span></div>')
    return sc.figure(
        code="FIG. 6 — ACCURACY LADDER",
        fid="fig-ladder",
        plot=fr.svg(),
        legend_html=leg,
        caption=("<b>Figure 6.</b> About fourteen orders of magnitude of position accuracy. Stellar parallax sits around an "
                 "au today and could reach a hundredth of an au; radio tracking is millions of times better while a "
                 "ground network can see the craft. Methods differ in what they measure (range, plane of sky, full "
                 "3D) and where they were tested; read the table view before comparing rows."),
        data_html=sc.table("Position accuracy by method", ["Method", "Accuracy", "In km", "Kind", "Source"], rows, numeric_from=1),
    )


if __name__ == "__main__":
    for name, fn in [("shift", fig_shift), ("lines", fig_lines), ("scatter", fig_scatter), ("dop", fig_dop),
                     ("pairs", fig_pairs), ("ladder", fig_ladder)]:
        sc.inject(PAGE, name, fn())
        print("built", name)