#!/usr/bin/env python3
"""Two world maps from IGRF-14 for docs/antarctica-magnetic/: declination and
total intensity at the surface, September 2026.

Top: declination, the angle between a compass needle and true north. The two
pinwheels are where every bearing converges on a dip pole; the heavy line is
the agonic line, where the needle happens to point true. His four Antarctic
sites (p. 454) are marked.

Bottom: total intensity. The two northern lobes over Canada and Siberia, the
single southern maximum south of Australia -- the strongest field on Earth --
and the South Atlantic Anomaly, the weakest.

    pip install ppigrf numpy matplotlib global_land_mask
    python3 magnetic_world_maps.py            -> img/fig-declination-intensity.png

Coastlines are the global_land_mask land/sea boundary, the same source the
site's other maps use; no basemap package is needed.
"""
import datetime as dt
import os
import sys

import numpy as np

try:
    from ppigrf import igrf
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from global_land_mask import globe
except ImportError as e:  # pragma: no cover
    sys.exit(f"{e}: pip install ppigrf numpy matplotlib global_land_mask")

DATE = dt.datetime(2026, 9, 1)
OUT = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "docs", "antarctica-magnetic", "img",
                   "fig-declination-intensity.png")
SITES = {"A": (-77.05, 168.15), "B": (-78.13, 110.52), "C": (-70.30, 24.45), "D": (-72.17, -74.02)}

lats = np.arange(-89.5, 90, 1.0)
lons = np.arange(-180, 180.5, 1.0)
LO, LA = np.meshgrid(lons, lats)
Be, Bn, Bu = (x[0] for x in igrf(LO, LA, 0, DATE))
D = np.degrees(np.arctan2(Be, Bn))
H = np.hypot(Be, Bn)
F = np.sqrt(Be ** 2 + Bn ** 2 + Bu ** 2) / 1000.0

# coastline: land mask on a finer grid
flat = np.arange(-89.75, 90, 0.5)
flon = np.arange(-180, 180.25, 0.5)
FLO, FLA = np.meshgrid(flon, flat)
LAND = globe.is_land(FLA, FLO).astype(float)


def coast(ax):
    ax.contour(FLO, FLA, LAND, levels=[0.5], colors="white", linewidths=0.55, alpha=0.9)
    ax.set_xlim(-180, 180); ax.set_ylim(-90, 90)
    ax.set_xticks(range(-180, 181, 60)); ax.set_yticks(range(-90, 91, 30))
    ax.set_xticklabels([f"{abs(x)}°{'W' if x < 0 else 'E' if x > 0 else ''}" for x in range(-180, 181, 60)], fontsize=8)
    ax.set_yticklabels([f"{abs(y)}°{'S' if y < 0 else 'N' if y > 0 else ''}" for y in range(-90, 91, 30)], fontsize=8)
    ax.set_aspect("equal")
    ax.tick_params(length=2, colors="#3a352d")
    for s in ax.spines.values():
        s.set_edgecolor("#bdb6a7")


def dip_poles():
    out = []
    for south in (True, False):
        m = (LA < -50) if south else (LA > 50)
        Hm = np.where(m, H, np.inf)
        i = np.unravel_index(Hm.argmin(), Hm.shape)
        out.append((LA[i], LO[i]))
    return out


fig, (a1, a2) = plt.subplots(2, 1, figsize=(10.5, 11.2), dpi=150, facecolor="white")
plt.subplots_adjust(left=0.06, right=0.9, top=0.95, bottom=0.05, hspace=0.24)

# ---- declination
# levels dense near zero so the 5-15 degree variation across a continent shows,
# coarse toward the poles where the needle swings through the whole circle
levels = [-180, -135, -90, -60, -40, -30, -20, -15, -10, -5, 0, 5, 10, 15, 20, 30, 40, 60, 90, 135, 180]
from matplotlib.colors import BoundaryNorm
cmap1 = plt.get_cmap("RdBu_r", len(levels) - 1)
c1 = a1.contourf(LO, LA, D, levels=levels, cmap=cmap1, norm=BoundaryNorm(levels, cmap1.N), extend="neither")
cl = a1.contour(LO, LA, D, levels=[-30, -20, -10, 10, 20, 30], colors="black", linewidths=0.3, alpha=0.6)
a1.clabel(cl, fmt="%d°", fontsize=6.5, inline=True, colors="black")
a1.contour(LO, LA, D, levels=[0], colors="black", linewidths=1.6)
coast(a1)
sp, npole = dip_poles()
for (la, lo), lab in ((sp, "south dip pole"), (npole, "north dip pole")):
    a1.plot(lo, la, marker="*", ms=13, mfc="#ffe27a", mec="black", mew=0.8, ls="none")
    a1.annotate(lab, (lo, la), xytext=(8, -12 if la < 0 else 6), textcoords="offset points", fontsize=8.5,
                fontweight="bold", color="black", bbox=dict(boxstyle="round,pad=0.2", fc="white", ec="none", alpha=0.85))
for k, (la, lo) in SITES.items():
    a1.plot(lo, la, marker="o", ms=6, mfc="white", mec="black", mew=0.9, ls="none")
    a1.annotate(k, (lo, la), xytext=(5, 4), textcoords="offset points", fontsize=8.5, fontweight="bold",
                bbox=dict(boxstyle="round,pad=0.15", fc="white", ec="none", alpha=0.85))
a1.annotate("agonic line: the needle points true here", (-92, 40), xytext=(-175, 62), textcoords="data", fontsize=8.5,
            arrowprops=dict(arrowstyle="-", lw=0.8, color="black"), bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="#bdb6a7"))
a1.annotate("every bearing converges\non the dip pole", (npole[1] - 40, 78), xytext=(60, 30), textcoords="data", fontsize=8.5,
            ha="center", arrowprops=dict(arrowstyle="-", lw=0.8, color="black"), bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="#bdb6a7"))
cb1 = fig.colorbar(c1, ax=a1, orientation="vertical", fraction=0.025, pad=0.015, ticks=[-180, -90, -40, -20, -10, 0, 10, 20, 40, 90, 180])
cb1.set_label("declination, degrees  (blue: needle west of true north; red: east)", fontsize=8.5)
cb1.ax.tick_params(labelsize=8)
a1.set_title("Where the compass points, relative to true north — IGRF-14, September 2026", fontsize=11, loc="left", color="#2E4057")

# ---- intensity
lev2 = np.arange(22, 68.1, 2)
c2 = a2.contourf(LO, LA, F, levels=lev2, cmap="turbo", extend="both")
cs = a2.contour(LO, LA, F, levels=np.arange(25, 66, 5), colors="black", linewidths=0.3, alpha=0.55)
a2.clabel(cs, fmt="%d", fontsize=6.5, inline=True, colors="black")
coast(a2)
for (la, lo), lab in ((sp, ""), (npole, "")):
    a2.plot(lo, la, marker="*", ms=11, mfc="#ffe27a", mec="black", mew=0.8, ls="none")


def mark(ax, cond, label, xy_text, fmt="max"):
    Fm = np.where(cond, F, -np.inf if fmt == "max" else np.inf)
    i = np.unravel_index(Fm.argmax() if fmt == "max" else Fm.argmin(), Fm.shape)
    ax.plot(LO[i], LA[i], marker="o", ms=5, mfc="white", mec="black", mew=0.9, ls="none")
    ax.annotate(f"{label}\n{F[i]:.1f} µT", (LO[i], LA[i]), xytext=xy_text, textcoords="data", fontsize=8.5, ha="center",
                arrowprops=dict(arrowstyle="-", lw=0.8, color="black"), bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="#bdb6a7"))
    return F[i], LA[i], LO[i]


mark(a2, (LA > 40) & (LO < -40), "Canadian lobe", (-150, 20))
mark(a2, (LA > 40) & (LO > 40), "Siberian lobe", (150, 25))
mark(a2, (LA < -40) & (LO > 90), "southern maximum —\nthe strongest field on Earth", (60, -45))
mark(a2, (LA < 0) & (LA > -50) & (LO < 0) & (LO > -90), "South Atlantic Anomaly —\nthe weakest", (-120, -50), fmt="min")
cb2 = fig.colorbar(c2, ax=a2, orientation="vertical", fraction=0.025, pad=0.015)
cb2.set_label("total intensity, µT", fontsize=8.5)
cb2.ax.tick_params(labelsize=8)
a2.set_title("How strong the field is — IGRF-14, September 2026", fontsize=11, loc="left", color="#2E4057")

os.makedirs(os.path.dirname(OUT), exist_ok=True)
fig.savefig(OUT, dpi=150)
print("wrote", os.path.normpath(OUT))
print("dip poles", sp, npole)