#!/usr/bin/env python3
"""The numbers behind "A great storm, not a weaker field, made multi-MeV electrons
inside the belts' barrier" (factchecked.science/inner-belt-electrons/).

What it does, with the Python standard library only:
  1. Reads the IGRF-14 dipole moment every five years, as computed for our guide to
     Earth's magnetic field (magnetic-field-guide-numbers.json, published at
     https://factchecked.science/magnetic-field-guide/), and interpolates it linearly
     to the dates of the storms on the page.
  2. Prints how much weaker the dipole was in May 2024 than at each earlier date.
  3. Prints the equatorial altitude of the L shells on the page's chart, in a dipole:
     altitude = (L - 1) x 6371.2 km, and the L shell of a few quoted altitudes.
  4. Prints the storms' Dst minima, which are constants here: they are the hourly
     minima of Kyoto's Dst index (World Data Center for Geomagnetism, Kyoto,
     doi:10.17593/14515-74000), final to 2020 and provisional after, as cited on the page.

Run:  python3 barrier_numbers.py                 (fetches the guide's file by URL)
      python3 barrier_numbers.py <local copy>    (uses a downloaded copy instead)
"""
import hashlib
import json
import sys
import urllib.request

GUIDE_URL = "https://factchecked.science/magnetic-field-guide/magnetic-field-guide-numbers.json"
# SHA-256 of the guide's file when this page was made (2026-10-29); a later revision
# of the guide may change it, in which case the script says so and carries on.
GUIDE_SHA256 = "cb6cfee622898a281e64b8e725dae1240d36f8d3b3bbafbeeb5fc92f101974fd"
RE_KM = 6371.2  # IGRF reference radius, km

# Dates (UT) as decimal years: middle of the day.
def decimal_year(year, day_of_year, leap):
    return year + (day_of_year - 0.5) / (366 if leap else 365)

DATES = [
    ("24 March 1991: CRRES sees electrons above 15 MeV injected into 2 < L < 3", decimal_year(1991, 83, False)),
    ("30 October 2003: the Halloween storms", decimal_year(2003, 303, False)),
    ("1 July 2013: middle of the 20 months in which the barrier was found", decimal_year(2013, 182, False)),
    ("17 March 2015: the deepest storm of 2012-2015; the barrier held", decimal_year(2015, 76, False)),
    ("24 March 2017: Arase starts full operations; its record begins", decimal_year(2017, 83, False)),
    ("11 May 2024: the May 2024 superstorm", decimal_year(2024, 132, True)),
]

# Kyoto Dst hourly minima, nT (final to 2020, provisional after)
DST = [
    ("March 1991", -298, "25 March 1991, 00-01 UT, final"),
    ("Halloween 2003", -383, "30 October 2003, 22-23 UT, final"),
    ("September 2012 to April 2014, the 20 months in which the barrier was found", -132, "17 March 2013, 20-21 UT, final"),
    ("September 2012 to December 2015", -234, "17 March 2015, 22-23 UT, final"),
    ("1996-2000, the years of the 37 storms behind Tsyganenko and Sitnov's 2005 model", -300, "16 July 2000, 00-01 UT, final"),
    ("24 March 2017 to 9 May 2024, Arase's record before the superstorm", -213, "24 April 2023, 05-06 UT, provisional"),
    ("May 2024", -406, "11 May 2024, 02-03 UT, provisional (-412 nT in the real-time index)"),
    ("October 2024, the next severe storm; Arase's electrons decayed gradually until it", -333, "11 October 2024, 01-02 UT, provisional"),
]

L_SHELLS = [1.5, 1.9, 2.0, 2.5, 2.8, 3.0, 3.5, 4.0]
ALTITUDES_KM = [("the press release's 12,000 km", 12000), ("Baker et al. 2004's 10,000 km", 10000)]


def load(arg):
    if arg:
        raw = open(arg, "rb").read()
        where = arg
    else:
        req = urllib.request.Request(GUIDE_URL, headers={"User-Agent": "barrier_numbers.py"})
        with urllib.request.urlopen(req, timeout=60) as r:
            raw = r.read()
        where = GUIDE_URL
    digest = hashlib.sha256(raw).hexdigest()
    note = "matches" if digest == GUIDE_SHA256 else "differs from"
    print("Input: %s\n  SHA-256 %s (%s the file used for the page)\n" % (where, digest, note))
    return json.loads(raw.decode("utf-8"))


def main():
    data = load(sys.argv[1] if len(sys.argv) > 1 else None)
    series = [(p["epoch"], p["M_1e22_Am2"]) for p in data["field_decline"]["dipole_series"]]

    def moment(t):
        for (e1, m1), (e2, m2) in zip(series, series[1:]):
            if e1 <= t <= e2:
                return m1 + (m2 - m1) * (t - e1) / (e2 - e1)
        raise ValueError("date outside the series: %s" % t)

    m2024 = moment(DATES[-1][1])
    print("IGRF-14 dipole moment, 10^22 A m^2, and how much weaker it was in May 2024")
    for label, t in DATES:
        m = moment(t)
        print("  %-75s %.4f  %5.2f%%" % (label, m, (1 - m2024 / m) * 100))

    print("\nKyoto Dst minima, nT")
    for label, value, when in DST:
        print("  %-75s %5d  (%s)" % (label, value, when))

    print("\nEquatorial altitude of a dipole L shell: (L - 1) x %.1f km" % RE_KM)
    for L in L_SHELLS:
        print("  L = %.1f  %6.0f km" % (L, (L - 1) * RE_KM))
    for label, km in ALTITUDES_KM:
        print("  %-35s L = %.2f" % (label, 1 + km / RE_KM))


if __name__ == "__main__":
    main()