#!/usr/bin/env python3
"""Antarctica's Magnetic North (pp. 415-473): what the compass measures, and
what a ring magnet under a flat plane would have it measure instead.

The book takes four Antarctic sites, reads their magnetic declination off the
NCEI geomagnetic calculator, and asks whether declination is "a concept created
to make our current maps appear accurate" -- whether the compass "was right all
along". This script does three things with that.

  1. Recomputes the four declinations from IGRF-14, the open reference model
     that the NCEI calculator's sibling (WMM) is built the same way from:
     ground observatories, repeat stations and the Swarm satellites. His four
     numbers are reproduced to a tenth of a degree. They are measurements.

  2. Prints what else the same model says at those sites and across the world:
     the dip (the angle the needle tilts out of the horizontal), the total
     intensity, where the two dip poles sit, where the field is strongest, and
     how it weakens with height. These are the quantities a compass, a dip
     needle, a phone magnetometer, an observatory and a satellite each measure.

  3. Builds the flat-plane alternative the map requires -- a magnet whose
     north pole is under the centre of the plane and whose south pole runs
     round the Antarctic rim -- and prints what THAT field would look like at
     the same places, for a range of magnet depths. Every axisymmetric ring
     puts declination at zero everywhere (the compass always points at the
     centre), which is the book's premise and is contradicted by its own
     table; and no depth reproduces the dip, the intensity ratios or the
     fall-off with height together.

    pip install ppigrf numpy
    python3 antarctica_magnetic.py

IGRF-14: Alken et al., Earth Planets Space (2021) for IGRF-13; the 14th
generation coefficients (2024) ship with ppigrf. Sign conventions: declination
positive east of true north; dip positive downward.
"""
import datetime as dt
import sys

import numpy as np

try:
    from ppigrf import igrf
except ImportError:  # pragma: no cover
    sys.exit("pip install ppigrf")

DATE = dt.datetime(2026, 9, 1)
A_EARTH = 6371.2  # km, IGRF reference radius


def dms(d, m):
    return d + m / 60.0


# His table, p. 454: site, latitude, longitude, declination as printed (E positive)
HIS_SITES = {
    "A": (-dms(77, 3), dms(168, 9), +dms(135, 7)),
    "B": (-dms(78, 8), dms(110, 31), -dms(130, 30)),
    "C": (-dms(70, 18), dms(24, 27), -dms(40, 10)),
    "D": (-dms(72, 10), -dms(74, 1), +dms(26, 18)),
}

PLACES = [
    ("Scott Base", -77.85, 166.77), ("Vostok", -78.46, 106.84),
    ("Dumont d'Urville", -66.66, 140.00), ("South Pole", -89.99, 0.0),
    ("Sydney", -33.87, 151.21), ("Cape Town", -33.93, 18.42),
    ("Singapore", 1.35, 103.82), ("Quito", -0.18, -78.47),
    ("Tampa", 27.97, -82.53), ("London", 51.51, -0.13),
    ("Boston", 42.36, -71.06), ("Seattle", 47.61, -122.33),
]


def elements(lon, lat, h_km=0.0, date=DATE):
    """Declination D (deg, +E), dip I (deg, +down), intensity F and horizontal H (nT)."""
    Be, Bn, Bu = (float(np.ravel(x)[0]) for x in igrf(lon, lat, h_km, date))
    H = float(np.hypot(Be, Bn))
    F = float(np.hypot(H, Bu))
    D = float(np.degrees(np.arctan2(Be, Bn)))
    I = float(np.degrees(np.arctan2(-Bu, H)))
    return D, I, F, H


def dip_pole(south, date=DATE):
    """Where the horizontal field vanishes: (H nT, lat, lon). Coarse-to-fine search."""
    lats = np.arange(-89.9, -50, 0.5) if south else np.arange(50, 89.9, 0.5)
    lons = np.arange(-180, 180, 1.0)
    for step in (1.0, 0.1, 0.01):
        LO, LA = np.meshgrid(lons, lats)
        Be, Bn, _ = igrf(LO, LA, 0, date)
        H = np.hypot(Be, Bn)[0]
        i = np.unravel_index(H.argmin(), H.shape)
        la, lo = float(LA[i]), float(LO[i])
        lats = np.arange(la - step, la + step, step / 10)
        lons = np.arange(lo - 2 * step, lo + 2 * step, step / 10)
    return float(H[i]), la, lo


def intensity_extremes(date=DATE):
    lats = np.arange(-89.5, 90, 1.0)
    lons = np.arange(-180, 180, 1.0)
    LO, LA = np.meshgrid(lons, lats)
    Be, Bn, Bu = igrf(LO, LA, 0, date)
    F = np.sqrt(Be ** 2 + Bn ** 2 + Bu ** 2)[0]
    i = np.unravel_index(F.argmax(), F.shape)
    j = np.unravel_index(F.argmin(), F.shape)
    Fn = np.where(LA > 0, F, 0)
    k = np.unravel_index(Fn.argmax(), F.shape)
    return {"max": (F[i], LA[i], LO[i]), "min": (F[j], LA[j], LO[j]), "max_north": (F[k], LA[k], LO[k])}


def zonal(lat, date=DATE):
    """Longitude-mean dip and intensity at a latitude."""
    la = np.sign(lat) * 89.99 if abs(lat) >= 90 else lat
    rows = np.array([elements(lo, la, 0, date) for lo in np.arange(-180, 180, 5.0)])
    return rows[:, 1].mean(), rows[:, 2].mean(), rows[:, 2].min(), rows[:, 2].max()


def gc_km(lat1, lon1, lat2, lon2):
    p1, p2 = np.radians(lat1), np.radians(lat2)
    dl = np.radians(lon2 - lon1)
    c = np.sin(p1) * np.sin(p2) + np.cos(p1) * np.cos(p2) * np.cos(dl)
    return A_EARTH * np.arccos(np.clip(c, -1, 1))


# ---------------------------------------------------------------- ring magnet
RIM_KM = 180 * 111.2  # centre to Antarctic rim on the azimuthal-equidistant plane
N_RING = 3600


def ring_field(rho, z, depth, ring_depth=None):
    """Field at distance rho from the centre, height z above the plane, for a magnet with
    one pole (strength -1: field lines go IN, as at Earth's north) at depth `depth` under
    the centre and the opposite pole (+1 total) spread uniformly round the rim at
    `ring_depth`. Coulomb pole model; units arbitrary. Returns (B_out, B_up)."""
    if ring_depth is None:
        ring_depth = depth
    P = np.array([rho, 0.0, z])
    r = P - np.array([0.0, 0.0, -depth])
    B = -r / np.linalg.norm(r) ** 3
    th = (np.arange(N_RING) + 0.5) * 2 * np.pi / N_RING
    src = np.stack([RIM_KM * np.cos(th), RIM_KM * np.sin(th), -ring_depth * np.ones(N_RING)], axis=1)
    rr = P[None, :] - src
    B = B + (rr / np.linalg.norm(rr, axis=1)[:, None] ** 3).sum(axis=0) / N_RING
    return float(B[0]), float(B[2])


def ring_elements(rho, depth, z=0.0, **kw):
    """Dip (deg, +down) and intensity (arbitrary) on the plane. Declination is identically
    zero: the horizontal component of an axisymmetric field is radial."""
    bo, bu = ring_field(rho, z, depth, **kw)
    F = float(np.hypot(bo, bu))
    I = float(np.degrees(np.arctan2(-bu, abs(bo))))
    return I, F


def ring_summary(depth):
    Ic, Fc = ring_elements(1.0, depth)
    Ie, Fe = ring_elements(90 * 111.2, depth)
    I70, F70 = ring_elements(160 * 111.2, depth)
    Ir, Fr = ring_elements(RIM_KM - 1.0, depth)
    rhos = np.linspace(100, RIM_KM - 100, 600)
    dips = np.array([ring_elements(x, depth)[0] for x in rhos])
    dip0_deg = float(rhos[np.argmin(abs(dips))] / 111.2)
    Ic450 = np.hypot(*ring_field(1.0, 450.0, depth))
    return {"dip0_deg_from_pole": dip0_deg, "F_eq_over_centre": Fe / Fc, "F_rim_over_centre": Fr / Fc,
            "F_70S_over_centre": F70 / Fc, "dip_70S": I70, "dip_rim": Ir, "centre_450km": Ic450 / Fc}


def bar_field(rho, z, depth, length):
    """A bar magnet standing on end under the centre: the ring shrunk to a point `length` below the near pole."""
    P = np.array([rho, 0.0, z])
    r1 = P - np.array([0.0, 0.0, -depth]); r2 = P - np.array([0.0, 0.0, -depth - length])
    B = -r1 / np.linalg.norm(r1) ** 3 + r2 / np.linalg.norm(r2) ** 3
    return float(B[0]), float(B[2])


def bar_summary(depth, length):
    def el(rho):
        bo, bu = bar_field(rho, 0.0, depth, length)
        return float(np.degrees(np.arctan2(-bu, abs(bo)))), float(np.hypot(bo, bu))
    Ic, Fc = el(1.0); Ie, Fe = el(90 * 111.2); I70, F70 = el(160 * 111.2); Ir, Fr = el(RIM_KM - 1.0)
    return {"F_eq_over_centre": Fe / Fc, "F_rim_over_centre": Fr / Fc, "dip_70S": I70, "dip_rim": Ir}


def dipole_450():
    return (A_EARTH / (A_EARTH + 450.0)) ** 3


def main():
    print("== 1. His four Antarctic sites, p. 454: NCEI value vs IGRF-14 for 2026.67")
    for k, (lat, lon, his) in HIS_SITES.items():
        D, I, F, H = elements(lon, lat)
        print(f"  {k}  {lat:8.3f} {lon:9.3f}   his D {his:+8.2f}   IGRF D {D:+8.2f}  dip {I:+6.1f}  F {F/1000:5.1f} uT  H {H/1000:5.1f} uT")

    print("\n== 2. The same model elsewhere")
    for name, lat, lon in PLACES:
        D, I, F, H = elements(lon, lat)
        print(f"  {name:18s} D {D:+7.1f}  dip {I:+6.1f}  F {F/1000:5.1f} uT  H {H/1000:5.1f} uT")

    sp, npole = dip_pole(True), dip_pole(False)
    print(f"\n  south dip pole  {sp[1]:.2f} {sp[2]:.2f}  (H = {sp[0]:.0f} nT)   F there {elements(sp[2], sp[1])[2]/1000:.1f} uT")
    print(f"  north dip pole  {npole[1]:.2f} {npole[2]:.2f}  (H = {npole[0]:.0f} nT)   F there {elements(npole[2], npole[1])[2]/1000:.1f} uT")
    print(f"  north dip pole 1831 (Ross, Cape Adelaide 70.09N 96.77W) to 2026: {gc_km(70.09, -96.77, npole[1], npole[2]):.0f} km")
    print(f"  south dip pole 1909 (David/Mawson/Mackay 72.42S 155.27E) to 2026: {gc_km(-72.42, 155.27, sp[1], sp[2]):.0f} km")
    ex = intensity_extremes()
    print(f"  strongest surface field   {ex['max'][0]/1000:.1f} uT at {ex['max'][1]:.1f} {ex['max'][2]:.1f}")
    print(f"  strongest in the north    {ex['max_north'][0]/1000:.1f} uT at {ex['max_north'][1]:.1f} {ex['max_north'][2]:.1f}")
    print(f"  weakest surface field     {ex['min'][0]/1000:.1f} uT at {ex['min'][1]:.1f} {ex['min'][2]:.1f}  (South Atlantic Anomaly)")

    print("\n  longitude-mean dip and intensity by latitude")
    for lat in (90, 60, 30, 0, -30, -60, -70, -80, -90):
        I, F, lo, hi = zonal(lat)
        print(f"    {lat:+3d}  dip {I:+6.1f}   F {F/1000:5.1f} uT  (range {lo/1000:.1f}-{hi/1000:.1f})")

    for y in (1902, 2026):
        D, I, F, H = elements(166.64, -77.85, 0, dt.datetime(y, 1, 1))
        print(f"  Hut Point, IGRF epoch {y}: D {D:+.1f}  dip {I:+.1f}  F {F/1000:.1f} uT")
    print(f"  Scott Base to his site A: {gc_km(-77.85, 166.77, HIS_SITES['A'][0], HIS_SITES['A'][1]):.0f} km")

    print("\n  fall-off with height, 0 -> 450 km (Swarm's orbit)")
    for name, lat, lon in (("Hartland", 51.0, -4.5), ("Scott Base", -77.85, 166.77), ("Gulf of Guinea", 0.0, 0.0)):
        F0, F450 = elements(lon, lat, 0)[2], elements(lon, lat, 450)[2]
        print(f"    {name:15s} {F0/1000:5.1f} -> {F450/1000:5.1f} uT   ratio {F450/F0:.3f}")
    print(f"    pure dipole (a/(a+h))^3 = {dipole_450():.3f}")

    print("\n== 3. A ring magnet under a flat plane: pole under the centre, opposite pole round the rim")
    print("   declination = 0 everywhere, for any depth: the horizontal field is radial")
    print("   depth   dip=0 at   F(eq)/F(N)  F(rim)/F(N)  F(70S)/F(N)  dip(70S)  dip(rim)  450km/surface at centre")
    for depth in (1000, 2000, 4178, 6371, 10000):
        s = ring_summary(depth)
        print(f"   {depth:5d}   {s['dip0_deg_from_pole']:5.1f}deg    {s['F_eq_over_centre']:.3f}       {s['F_rim_over_centre']:.3f}        "
              f"{s['F_70S_over_centre']:.3f}      {s['dip_70S']:+6.1f}   {s['dip_rim']:+6.1f}      {s['centre_450km']:.3f}")
    b = bar_summary(4178, 10000)
    print(f"   bar magnet on end, top at 4178 km, 10000 km long: F(eq)/F(N) {b['F_eq_over_centre']:.3f}  F(rim)/F(N) {b['F_rim_over_centre']:.3f}  dip(70S) {b['dip_70S']:+.1f}")
    print("   measured (IGRF): dip=0 near 0-10 deg S of the equator; F(eq)/F(N) about 0.55-0.6;")
    print("   F(south max)/F(north max) above 1; dip at 70S about -70; 450 km ratio about 0.8")


if __name__ == "__main__":
    main()