# SPDX-License-Identifier: MIT
"""YS-EN-0002 · Plan your own parallax-navigation experiment.

Predicts how well a setup can locate the camera, before you build it, using
the same line-of-position solver as the note (parallax_fix.fix). Edit the
SETUPS below for your own "stars" and camera, then run:

    python plan_experiment.py

For each setup it prints the 1-sigma error ellipsoid, then simulates noisy
observations and solves them, so you can see the scatter you should expect.
Units are whatever you use for positions (metres on a tabletop, au in space);
the camera's angular precision is in arcseconds.
"""
from __future__ import annotations

import math

import numpy as np

from parallax_fix import GAIA, fix, observe, perturb, star_position

ARCSEC = math.pi / 648_000


def tabletop():
    """Two LEDs 3-4 m from a camera; units are metres."""
    leds = np.array([[3.0, 1.0, 0.3], [0.8, 3.5, -0.2]])
    camera = np.array([0.0, 0.0, 0.0])
    # A phone camera: ~70° across 4000 px is ~63"/px; centroiding to ~0.1 px
    # gives ~6". Lens distortion is usually worse: calibrate it, or keep the
    # LEDs near the centre of the frame. 30" is a cautious starting point.
    return "Tabletop: two LEDs, phone camera", leds, camera, 30.0, "m"


def earth_orbit():
    """Barnard's Star + Lalande 21185 seen from Earth; units are au."""
    stars = np.array([star_position(*GAIA[n]) for n in ("Barnard's Star", "Lalande 21185")])
    earth = np.array([1.0, 0.0, 0.0])  # any point on Earth's orbit; the answer barely depends on it
    # 0.05" per star assumes many frames averaged over several nights.
    return "Earth's orbit: Barnard's Star + Lalande 21185", stars, earth, 0.05, "au"


def run(name, stars, truth, sigma_arcsec, unit, trials=2000, seed=1):
    sigma = sigma_arcsec * ARCSEC
    dirs = observe(stars, truth)
    # fix() weights by |s|; for a tabletop the camera is not tiny compared
    # with the star distances, so weight by the true ranges instead.
    rel = stars - truth
    _, cov = fix(rel, dirs, sigma)
    axes = np.sqrt(np.linalg.eigvalsh(cov))[::-1]
    rng = np.random.default_rng(seed)
    errs = []
    for _ in range(trials):
        r, _ = fix(rel, perturb(dirs, sigma, rng), sigma)
        errs.append(np.linalg.norm(r))
    errs = np.array(errs)
    print(f"\n{name}")
    print(f"  pointing precision    {sigma_arcsec:g}\" per star")
    print(f"  1-sigma ellipsoid     {axes[0]:.3g} / {axes[1]:.3g} / {axes[2]:.3g} {unit}")
    print(f"  simulated fixes       median error {np.median(errs):.3g} {unit}, 95% within {np.quantile(errs, 0.95):.3g} {unit}")


if __name__ == "__main__":
    for setup in (tabletop, earth_orbit):
        run(*setup())