"""Figures for topocentric.qmd.

img/fig-topo-latitude.png : a meridian section of the WGS-84 ellipsoid (flattening exaggerated): the ellipsoid
normal at the site (geodetic latitude phi), the radius from the centre (geocentric latitude phi_gc), and
the height h.
img/fig-topo-sez.png      : the SEZ triad at a site, the line of sight to a satellite and the angles azimuth
(from north, clockwise) and elevation.

Run from the sfm-slides directory:   python3 codes/fig-topo-frames.py
"""
import numpy as np
import matplotlib.pyplot as plt
from _manstyle import INK, MAROON, BLUE, GREEN, AMBER, PURPLE

# ------------------------------------------------------------------ latitude
f_draw = 0.12
a, b = 1.0, 1.0 - f_draw
t = np.linspace(0, 2 * np.pi, 400)
phi = np.radians(50.0)
N = a / np.sqrt(1 - (2 * f_draw - f_draw ** 2) * np.sin(phi) ** 2)
e2 = 2 * f_draw - f_draw ** 2
P = np.array([N * np.cos(phi), N * (1 - e2) * np.sin(phi)])           # surface point
nrm = np.array([np.cos(phi), np.sin(phi)])
hgt = 0.22
Q = P + hgt * nrm                                                      # site above the surface
phi_gc = np.arctan2(Q[1], Q[0])
fig, ax = plt.subplots(figsize=(7.6, 6.4))
ax.plot(a * np.cos(t), b * np.sin(t), color=BLUE, lw=2)
ax.plot(a * np.cos(t), a * np.sin(t), color="#bdb6aa", lw=0.9, ls=":")
ax.plot([-1.25, 1.55], [0, 0], color=INK, lw=0.8); ax.plot([0, 0], [-1.05, 1.5], color=INK, lw=0.8)
ax.text(1.5, -0.08, "equator", fontsize=10); ax.text(0.03, 1.45, "polar axis", fontsize=10)
foot = N * e2 * np.sin(phi)                                            # where the normal meets the axis (below the centre)
axis_pt = np.array([0.0, -foot])
ax.plot([axis_pt[0], Q[0] + 0.2 * nrm[0]], [axis_pt[1], Q[1] + 0.2 * nrm[1]], color=GREEN, lw=1.8)
ax.plot([0, Q[0]], [0, Q[1]], color=MAROON, lw=1.8)
ax.scatter(*Q, color=INK, s=45, zorder=5)
ax.scatter(*P, color=BLUE, s=25, zorder=5)
ax.text(Q[0] + 0.05, Q[1] + 0.05, "site", fontsize=11)
# angles
th = np.linspace(0, phi, 40); ax.plot(0.55 * np.cos(th), 0.55 * np.sin(th) - foot * 0.0, color=MAROON, lw=1.2)
ax.text(0.58, 0.18, r"$\varphi_{gc}$", color=MAROON, fontsize=13)
th2 = np.linspace(0, phi, 40); ax.plot(axis_pt[0] + 0.85 * np.cos(th2), axis_pt[1] + 0.85 * np.sin(th2), color=GREEN, lw=1.2)
ax.text(0.9, -0.12 - foot, r"$\varphi$", color=GREEN, fontsize=13)
ax.annotate("", xy=Q, xytext=P, arrowprops=dict(arrowstyle="-", color=AMBER, lw=3))
ax.text(Q[0] + 0.08, Q[1] - 0.12, "$h$", color=AMBER, fontsize=13)
ax.text(0.02, -0.28 - foot * 0, "$N$", color=GREEN, fontsize=12)
ax.text(-1.2, -1.02, "green: ellipsoid normal, geodetic latitude $\\varphi$", color=GREEN, fontsize=10.5)
ax.text(-1.2, -1.14, "maroon: radius from the centre, geocentric latitude $\\varphi_{gc}$", color=MAROON, fontsize=10.5)
ax.set_aspect("equal"); ax.axis("off"); ax.set_xlim(-1.25, 1.7); ax.set_ylim(-1.2, 1.55)
ax.set_title("A meridian section of the ellipsoid (flattening drawn $\\approx 40\\times$ too large)", fontsize=11.5)
fig.savefig("img/fig-topo-latitude.png", dpi=150, bbox_inches="tight"); plt.close(fig)

# ------------------------------------------------------------------ SEZ
fig = plt.figure(figsize=(8.2, 6.6))
ax = fig.add_subplot(111, projection="3d")
th = np.linspace(0, 2 * np.pi, 200)
ax.plot(np.cos(th), np.sin(th), 0 * th, color="#bdb6aa", lw=0.9)              # horizon circle
ax.plot_surface(*np.meshgrid(np.linspace(-1, 1, 2), np.linspace(-1, 1, 2)), np.zeros((2, 2)), color="#cfdcec", alpha=0.25, lw=0)
# axes: x = east, y = north, z = zenith (a right-handed view); S = -north
E_, N_, Z_ = np.array([1, 0, 0]), np.array([0, 1, 0]), np.array([0, 0, 1])
for v, lab, col in ((E_, r"$\hat{\mathbf{E}}$ (east)", GREEN), (N_, "north $=-\\hat{\\mathbf{S}}$", AMBER), (Z_, r"$\hat{\mathbf{Z}}$ (zenith)", PURPLE)):
    ax.quiver(0, 0, 0, *(1.2 * v), color=col, lw=2.4, arrow_length_ratio=0.1)
    ax.text(*(1.32 * v + np.array([0, 0, 0.05])), lab, color=col, fontsize=11)
az, el = np.radians(55), np.radians(35)
u = np.array([np.cos(el) * np.sin(az), np.cos(el) * np.cos(az), np.sin(el)])  # (east, north, up)
ax.quiver(0, 0, 0, *(1.6 * u), color=MAROON, lw=2.8, arrow_length_ratio=0.08)
ax.text(*(1.7 * u + np.array([0, 0, 0.08])), r"line of sight, range $\rho$", color=MAROON, fontsize=11)
foot = np.array([u[0], u[1], 0]) * 1.6
ax.plot([foot[0], 1.6 * u[0]], [foot[1], 1.6 * u[1]], [0, 1.6 * u[2]], color=MAROON, lw=1, ls=":")
a_ = np.linspace(0, az, 40); ax.plot(0.7 * np.sin(a_), 0.7 * np.cos(a_), 0 * a_, color=AMBER, lw=2)
ax.text(0.45, 0.78, 0, r"$A$", color=AMBER, fontsize=14)
e_ = np.linspace(0, el, 40); ax.plot(0.9 * np.cos(e_) * np.sin(az), 0.9 * np.cos(e_) * np.cos(az), 0.9 * np.sin(e_), color=MAROON, lw=2)
ax.text(*(1.0 * np.array([np.cos(el / 2) * np.sin(az), np.cos(el / 2) * np.cos(az), np.sin(el / 2)]) + np.array([0.12, 0, 0])), r"$\varepsilon$", color=MAROON, fontsize=14)
ax.set_xlim(-1.1, 1.8); ax.set_ylim(-1.1, 1.8); ax.set_zlim(-0.1, 1.5)
ax.set_box_aspect((1, 1, 0.6), zoom=1.5); ax.view_init(26, -40); ax.set_axis_off()
ax.set_title(r"The topocentric frame: azimuth $A$ from north through east, elevation $\varepsilon$ above the horizon", fontsize=11.5)
fig.savefig("img/fig-topo-sez.png", dpi=150, bbox_inches="tight")
print("wrote img/fig-topo-latitude.png, img/fig-topo-sez.png")