"""Figure for lunar-trajectories-bates.qmd (img/fig-lun-noncoplanar.png): the geometry of a non-coplanar lunar launch.

(a) Declination of the point reached after a total geocentric sweep psi_t, from a site at declination (latitude) 28.5 N, for launch azimuths
    from 40 to 115 degrees (the Eastern Test Range window of the book). At psi_t = 180 deg every azimuth arrives at the same declination, -28.5.
(b) The launch azimuth needed to arrive at a given declination, for four sweep angles; the shaded band is the azimuth window.

Run from the sfm-slides directory:   python3 codes/fig-lun-noncoplanar.py
"""
import numpy as np
import matplotlib.pyplot as plt
from _manstyle import INK, MAROON, BLUE, GREEN, AMBER, PURPLE, despine
from _lunar import intercept_declination, launch_azimuth_for

D2R = np.pi / 180
d0 = 28.5 * D2R
fig, ax = plt.subplots(1, 2, figsize=(13.4, 5.4))
A = ax[0]
ps = np.linspace(0, 360, 721) * D2R
lo = np.degrees(intercept_declination(d0, ps, 40 * D2R)); hi = np.degrees(intercept_declination(d0, ps, 115 * D2R))
A.fill_between(np.degrees(ps), np.minimum(lo, hi), np.maximum(lo, hi), color="#dfe9f3", label="reachable with A between 40° and 115°")
cols = {40: MAROON, 55: AMBER, 70: GREEN, 90: BLUE, 115: PURPLE}
for az, c in cols.items():
    A.plot(np.degrees(ps), np.degrees(intercept_declination(d0, ps, az * D2R)), color=c, lw=1.9, label=f"$A_0={az}^\\circ$")
for y in (28.5, -28.5): A.axhline(y, color=INK, lw=0.9, ls=":")
A.axvline(180, color=INK, lw=0.8, ls=":")
A.scatter([180], [-28.5], color=INK, s=40, zorder=8); A.text(184, -33, "$\\psi_t=180^\\circ$: every azimuth ends at $-28.5^\\circ$", fontsize=10)
A.set_xlim(0, 360); A.set_ylim(-60, 60); A.set_xticks(range(0, 361, 60))
A.set_xlabel(r"total geocentric sweep angle $\psi_t$ [deg]"); A.set_ylabel(r"declination at intercept $\delta_1$ [deg]")
A.legend(loc="upper right", fontsize=9, frameon=False, ncol=2)
A.set_title("(a)  declination reached after a sweep $\\psi_t$", fontsize=11.5, loc="left"); despine(A)
B = ax[1]
d1 = np.linspace(-28.5, 28.5, 400) * D2R
for pt, c in ((120, BLUE), (150, GREEN), (180, AMBER), (210, MAROON)):
    az = np.array([launch_azimuth_for(d0, d, pt * D2R) for d in d1], dtype=object)
    ok = np.array([a is not None for a in az])
    y = np.array([np.degrees(a) if a is not None else np.nan for a in az], dtype=float)
    B.plot(np.degrees(d1), y, color=c, lw=2.2, label=f"$\\psi_t={pt}^\\circ$")
B.axhspan(40, 115, color="#dfe9f3", zorder=0); B.text(-27, 108, "azimuth window", fontsize=10, color=INK)
B.set_xlabel(r"declination at intercept $\delta_1$ [deg]"); B.set_ylabel(r"launch azimuth $A_0$ [deg]"); B.set_ylim(0, 180); B.set_xlim(-30, 30)
B.legend(loc="upper right", fontsize=9.5, frameon=False)
B.set_title("(b)  azimuth needed to arrive at a given declination", fontsize=11.5, loc="left"); despine(B)
fig.tight_layout(); fig.savefig("img/fig-lun-noncoplanar.png", dpi=150, bbox_inches="tight")
print("ok")