#!/usr/bin/env python3
"""
The decline forecast in "When is the next magnetic excursion due?"
(Fact-Checked Science, https://factchecked.science/next-excursion/, section 3).

Earth's axial dipole, the main part of its magnetic field, is measured by the Gauss coefficient g10.
This script takes g10 from the International Geomagnetic Reference Field, 14th generation (IGRF-14),
and extends its decline five ways. Each projection is followed until the dipole reaches 5% of its 2025 strength,
Laschamp's low, and the script also gives when each one halves the dipole, the level at which past excursions began:
  - straight lines at 24, 16 and 12.6 nT a year: the fastest five-year pace since 1945 (1970-75),
    the 1840-2010 average, and the IGRF-14 forecast pace for 2025-30
  - exponential decay that loses 5% of its value each century
  - an exponential speed-up that starts at 12.6 nT a year and doubles its pace every 70 years
It also gives the loss since 1840, every five-year pace since 1900, how long 16 nT a year would take to
reach zero, what halving the dipole or a Laschamp-level low by the 2040s would need, and the series drawn in the
page's chart.

Run it
  python3 decline.py
  Needs Python 3 and nothing else. It reads igrf14coeffs.txt from the folder it runs in. If the file
  isn't there, it downloads it from NOAA, which hosts the official IGRF-14 coefficients (1900-2025, plus
  the secular variation for 2025-30):
      https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt
  It prints the results and writes decline.json. To use other file names:
      python3 decline.py <IGRF-14 coefficient file> <output JSON>

Numbers taken from papers (numbered as in the page's source list)
  16 nT a year  the axial dipole's mean decay rate from 1840 to 2010 (Finlay, Aubert and Gillet 2016 [1], p. 2)
  2 nT a year   the 1590-1840 trend, rounded. Gubbins, Jones and Finlay (2006 [3]) found 2.28 +/- 2.72 nT a
                year, which is statistically indistinguishable from zero, so the 70-year doubling time
                derived from it is our own construct and only loosely constrained
  50%           about half today's strength: excursion directions began when the dipole had fallen to about
                half its normal value (Laj and Kissel 2015 [29]), and field models flag excursions below about
                half the present dipole (Korte et al. 2019 [18], section 4.4)
  5%            Laschamp's low relative to today's dipole, in the LSMOD.2 model (Korte et al. 2019 [18]) as
                Larsen et al. (2026 [6]) compute it

Conventions
  The 1840 value is the IGRF-14 value for 2010 plus 170 years at 16 nT a year; the 1590 value adds 250
  years at 2 nT a year. Projections start from |g10| in 2025 (29,350 nT) and stop at 5% of it.
  Positive paces mean weakening.
"""
import hashlib
import json
import math
import os
import sys
import urllib.request

IGRF_URL = "https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt"
# SHA-256 of the file we used on October 11, 2026. A different hash means NOAA has changed the file.
IGRF_SHA256 = "8f8d88403028fc4ee92c4f38d97b46e0a87e2cfc496045b43c9e26c1d6b0903c"


def fetch_igrf(path):
    """Download the IGRF-14 coefficient file if it isn't already here, then check it is the one we used."""
    if not os.path.exists(path):
        print(f"Downloading {IGRF_URL} ...")
        try:
            with urllib.request.urlopen(IGRF_URL, timeout=60) as resp:
                data = resp.read()
        except OSError as err:
            sys.exit(f"Couldn't download the IGRF-14 file ({err}).\n"
                     f"Save {IGRF_URL} as {path} in this folder and run the script again.")
        with open(path, "wb") as fh:
            fh.write(data)
    with open(path, "rb") as fh:
        digest = hashlib.sha256(fh.read()).hexdigest()
    if digest != IGRF_SHA256:
        print(f"Note: {path} differs from the file we used (SHA-256 {IGRF_SHA256[:16]}...), so results may differ.")


def read_igrf(path):
    """Return the epochs (1900.0 ... 2025.0) and every coefficient row, keyed (g/h, degree, order)."""
    years = None
    rows = {}
    with open(path) as fh:
        for line in fh:
            if line.startswith('g/h'):
                parts = line.split()
                years = parts[3:]          # '1900.0' ... '2025.0', '2025-30'
            elif line[:1] in 'gh' and not line.startswith('g/h'):
                p = line.split()
                rows[(p[0], int(p[1]), int(p[2]))] = [float(x) for x in p[3:]]
    epochs = [float(y) for y in years[:-1]]
    return epochs, rows


def main(igrf_path, out_path):
    fetch_igrf(igrf_path)
    epochs, rows = read_igrf(igrf_path)
    g10 = rows[('g', 1, 0)]
    sv = g10[-1]                     # 2025-30 secular variation, nT/yr
    G = {int(e): abs(v) for e, v in zip(epochs, g10[:-1])}
    g2025 = G[2025]
    out = {'source_file': os.path.basename(igrf_path), 'g10_abs_nT': G, 'sv_g10_2025_2030_nT_per_yr': sv}

    # five-year decay rates (positive = weakening)
    rates = []
    for a, b in zip(epochs[:-1], epochs[1:]):
        rates.append({'from': int(a), 'to': int(b), 'nT_per_yr': round((G[int(a)] - G[int(b)]) / (b - a), 2)})
    out['five_year_rates'] = rates
    post45 = [r for r in rates if r['from'] >= 1945]
    out['since_1945_min'] = min(post45, key=lambda r: r['nT_per_yr'])
    out['since_1945_max'] = max(post45, key=lambda r: r['nT_per_yr'])
    out['pre_1945_max'] = max([r for r in rates if r['from'] < 1945], key=lambda r: r['nT_per_yr'])
    out['mean_rate_1900_2025'] = round((G[1900] - g2025) / 125, 2)
    out['mean_rate_1900_2010'] = round((G[1900] - G[2010]) / 110, 2)
    r1980 = [r for r in rates if r['from'] in (1975, 1980)]
    out['rate_around_1980'] = r1980
    starts = [r['from'] for r in rates]
    out['ratio_2020_25_to_1980_85'] = round(rates[starts.index(2020)]['nT_per_yr'] / rates[starts.index(1980)]['nT_per_yr'], 3)
    out['ratio_sv2025_to_1980_85'] = round(sv / rates[starts.index(1980)]['nT_per_yr'], 3)
    out['ratio_sv2025_to_1975_80'] = round(sv / rates[starts.index(1975)]['nT_per_yr'], 3)

    # the 1840 and 1590 values, built from the published average rates
    g1840 = G[2010] + 170 * 16
    g1590 = g1840 + 250 * 2
    out['g1840_from_16nT'] = g1840
    out['g1590_from_2nT'] = g1590
    out['pct_of_2025_in_1840'] = round(100 * g1840 / g2025, 2)
    out['pct_of_2025_in_1590'] = round(100 * g1590 / g2025, 2)
    loss = 1 - g2025 / g1840
    out['loss_since_1840_pct'] = round(100 * loss, 2)
    out['loss_per_century_linear_pct'] = round(100 * loss / 1.85, 2)
    out['loss_per_century_compound_pct'] = round(100 * (1 - (g2025 / g1840) ** (100 / 185)), 2)
    out['exp_5pct_per_century_loss_over_185yr_pct'] = round(100 * (1 - 0.95 ** 1.85), 2)

    # zero crossing at 16 nT/yr
    out['years_to_zero_at_16_from_2025'] = round(g2025 / 16, 1)
    out['year_zero_at_16_from_2025'] = round(2025 + g2025 / 16, 1)
    out['years_to_zero_at_16_from_2010'] = round(G[2010] / 16, 1)
    out['years_to_zero_at_16_from_1840'] = round(g1840 / 16, 1)

    # straight lines to 5% of the 2025 value
    lines = {}
    for label, r in [('fastest_1970s', 24.0), ('mean_1840_2010', 16.0), ('sv_2025_30', sv),
                     ('rate_2020_25', rates[-1]['nT_per_yr']), ('slowest_2015_20', 7.61)]:
        lines[label] = {'rate': r, 'year_5pct': round(2025 + 0.95 * g2025 / r, 1),
                        'year_60pct': round(2025 + 0.40 * g2025 / r, 1),
                        'year_50pct': round(2025 + 0.50 * g2025 / r, 1),
                        'year_zero': round(2025 + g2025 / r, 1)}
    out['straight_lines'] = lines

    # exponential decay, 5% per century
    t5 = 100 * math.log(0.05) / math.log(0.95)
    out['exp_decay'] = {'year_5pct': round(2025 + t5, 1),
                        'pct_at_4500': round(100 * 0.95 ** ((4500 - 2025) / 100), 2),
                        'year_60pct': round(2025 + 100 * math.log(0.6) / math.log(0.95), 1),
                        'year_50pct': round(2025 + 100 * math.log(0.5) / math.log(0.95), 1)}
    # exponential decay at the measured 1840-2025 compound rate
    k = (g2025 / g1840) ** (1 / 185)
    out['exp_decay_measured_rate'] = {'loss_per_century_pct': round(100 * (1 - k ** 100), 2),
                                      'year_5pct': round(2025 + math.log(0.05) / math.log(k), 1)}

    # exponential speed-up: the pace starts at r0 and doubles every D years
    def speedup_year(r0, D, level=0.05):
        return 2025 + D * math.log2(1 + (1 - level) * g2025 / (r0 * D / math.log(2)))
    out['speedup'] = {'r0': sv, 'D': 70, 'year_5pct': round(speedup_year(sv, 70), 1),
                      'year_60pct': round(speedup_year(sv, 70, 0.6), 1),
                      'year_50pct': round(speedup_year(sv, 70, 0.5), 1),
                      'sensitivity': {f'D={D}': round(speedup_year(sv, D), 1) for D in (60, 70, 75, 77, 100, 125)}}

    # where the 70-year doubling comes from: ~2 nT/yr (1590-1840, midpoint 1715) -> 16 nT/yr (1840-2010, midpoint 1925)
    def dbl(r1, r2, span=210):
        return span / math.log2(r2 / r1)
    out['doubling_time'] = {'2_to_16': round(dbl(2, 16), 1),
                            '2.28_to_16': round(dbl(2.28, 16), 1),
                            '2.28_to_15': round(dbl(2.28, 15), 1),
                            '5.0_(2.28+1sd)_to_16': round(dbl(5.0, 16), 1),
                            'note': 'Gubbins et al. 2006 slope is 2.28 +/- 2.72 nT/yr: consistent with zero, so the doubling time is not constrained from below'}

    # what a Laschamp-level low (5% of 2025), or halving the dipole, by the 2040s would need
    need = {}
    for yr in (2040, 2045, 2049):
        r5 = 0.95 * g2025 / (yr - 2025)
        r0 = g2025 / (yr - 2025)
        r50 = 0.50 * g2025 / (yr - 2025)
        need[str(yr)] = {'rate_to_5pct_nT_per_yr': round(r5, 1), 'x_sv2025': round(r5 / sv, 1), 'x_fastest24': round(r5 / 24, 1),
                         'rate_to_zero_nT_per_yr': round(r0, 1),
                         'rate_to_50pct_nT_per_yr': round(r50, 1), 'x_sv2025_50pct': round(r50 / sv, 1),
                         'x_fastest24_50pct': round(r50 / 24, 1)}
    lo, hi = 0.1, 70.0               # doubling time the speed-up would need to reach 5% by 2045 (bisection)
    for _ in range(100):
        mid = (lo + hi) / 2
        if speedup_year(sv, mid) > 2045:
            hi = mid
        else:
            lo = mid
    need['speedup_doubling_time_to_reach_5pct_by_2045_yr'] = round(lo, 2)
    lo, hi = 0.1, 70.0               # and to halve the dipole by 2045
    for _ in range(100):
        mid = (lo + hi) / 2
        if speedup_year(sv, mid, 0.5) > 2045:
            hi = mid
        else:
            lo = mid
    need['speedup_doubling_time_to_reach_50pct_by_2045_yr'] = round(lo, 2)
    out['twenty_forties'] = need

    # chart-ready series (percent of the 2025 axial dipole), as drawn on the page
    pct = lambda gval: round(100 * gval / g2025, 2)
    series = {'reconstructed_1590_1840': [[1590, pct(g1590)], [1840, pct(g1840)]],
              'measured_1840_2025': [[1840, pct(g1840)]] + [[y, pct(G[y])] for y in sorted(G)]}
    for label, r in [('line_24', 24.0), ('line_16', 16.0), ('line_12.6', sv)]:
        end = 2025 + 0.95 * g2025 / r
        series[label] = [[2025, 100.0], [round(end, 1), 5.0]]
    series['exp_decay_5pct_per_century'] = [[y, round(100 * 0.95 ** ((y - 2025) / 100), 2)] for y in range(2025, 4501, 25)]
    # the same curve continued until it reaches 5% (the page's chart draws it to its end)
    T_exp = out['exp_decay']['year_5pct']
    series['exp_decay_to_5pct'] = [[y, round(100 * 0.95 ** ((y - 2025) / 100), 2)] for y in range(2025, int(T_exp), 25)] + [[T_exp, 5.0]]
    T = out['speedup']['year_5pct']
    sp = []
    for y in list(range(2025, 2340, 5)) + [T]:
        loss = sv * (70 / math.log(2)) * (2 ** ((y - 2025) / 70) - 1)
        sp.append([round(y, 1), round(100 * (g2025 - loss) / g2025, 2)])
    series['exp_speedup_double_70yr'] = sp
    series['reference_lines_percent'] = {'long_term_average_Wang2015': 60, 'laschamp_low_Larsen2026': 5,
                                         'half_strength_excursions_began_LajKissel2015_Korte2019': 50}
    out['chart_series'] = series

    with open(out_path, 'w') as fh:
        json.dump(out, fh, indent=1)

    # a readable summary
    sl = out['straight_lines']
    print(f"Axial dipole |g10|: {G[1900]:,.0f} nT in 1900, {g2025:,.0f} nT in 2025; IGRF-14 forecast pace for 2025-30: {sv} nT a year")
    print(f"1840 value (2010 + 170 years at 16 nT a year): {g1840:,.0f} nT = {out['pct_of_2025_in_1840']}% of 2025")
    print(f"Lost since 1840: {out['loss_since_1840_pct']}% ({out['loss_per_century_linear_pct']}% a century)")
    print(f"Five-year paces since 1945: {out['since_1945_min']['nT_per_yr']} ({out['since_1945_min']['from']}-{out['since_1945_min']['to']})"
          f" to {out['since_1945_max']['nT_per_yr']} nT a year ({out['since_1945_max']['from']}-{out['since_1945_max']['to']})")
    print(f"Pace 2020-25 vs 1980-85: {out['ratio_2020_25_to_1980_85']}; forecast 2025-30 vs 1980-85: {out['ratio_sv2025_to_1980_85']}")
    print(f"At 16 nT a year, zero in {out['years_to_zero_at_16_from_2025']} years from 2025 (AD {out['year_zero_at_16_from_2025']})")
    print("Year each projection halves the dipole (50%) and reaches Laschamp's low (5% of the 2025 strength):")
    print(f"  straight line, 24 nT a year:            AD {sl['fastest_1970s']['year_50pct']} and AD {sl['fastest_1970s']['year_5pct']}")
    print(f"  straight line, 16 nT a year:            AD {sl['mean_1840_2010']['year_50pct']} and AD {sl['mean_1840_2010']['year_5pct']}")
    print(f"  straight line, {sv} nT a year:          AD {sl['sv_2025_30']['year_50pct']} and AD {sl['sv_2025_30']['year_5pct']}")
    print(f"  exponential decay, 5% a century:        AD {out['exp_decay']['year_50pct']} and AD {out['exp_decay']['year_5pct']} ({out['exp_decay']['pct_at_4500']}% left in AD 4500)")
    print(f"  speed-up, pace doubling every 70 years: AD {out['speedup']['year_50pct']} and AD {out['speedup']['year_5pct']}")
    n45 = need['2045']
    print(f"To halve the dipole by 2045: {n45['rate_to_50pct_nT_per_yr']} nT a year, {n45['x_sv2025_50pct']} times the 2025-30 forecast pace"
          f" and {n45['x_fastest24_50pct']} times the fastest five-year pace since 1945")
    print(f"To reach 5% by 2045: {n45['rate_to_5pct_nT_per_yr']} nT a year, {n45['x_sv2025']} times the 2025-30 forecast pace"
          f" and {n45['x_fastest24']} times the fastest five-year pace since 1945")
    print(f"Wrote {out_path}")


if __name__ == '__main__':
    args = sys.argv[1:]
    main(args[0] if len(args) > 0 else 'igrf14coeffs.txt', args[1] if len(args) > 1 else 'decline.json')