#!/usr/bin/env python3 """Re-running the numbers of Yuan et al. (2026), "The Qianlong XXXV event", Front. Astron. Space Sci. 13:1857812, https://doi.org/10.3389/fspas.2026.1857812 (CC BY 4.0), and testing how they depend on its choices. Fact-Checked Science, https://factchecked.science/qianlong-1770-storm/ (27 October 2026). Standard library only; it reads no file and needs no network. Inputs are numbers printed in the paper (locators in brackets) and in the papers it is compared with: Eq. 1 (viewing geometry, §2.2): (a + h) cos(lam - lam0) = a + (a + h) sin(lam - lam0) tan(beta), solved exactly as lam - lam0 = arccos(a cos(beta) / (a + h)) - beta. a = 6371 km. Dongtinghu: lam0 = 19.2 deg, beta = 60/75/90 deg, h = 240 km (Fig. 2a) -> paper: 20.4 / 19.8 / 19.2 (§2.2). The paper's one 1859 point, "the US Frigate Sabine off Panama" (§3.2): lam0 = 20.45 deg, beta = 20 deg, an "assumed elevation angle" (Fig. 2b caption), h = 400 km -> paper: 28.30 (§3.2). What the ships' logs say (Hayakawa et al. 2018, ApJ 869:57, §3.2, §3.3 and Appendix 2): the Panama sighting was the USS Saranac's, at 04:00 on 29 August 1859 (the storm before Carrington's), 20.2 deg MLAT, with no elevation given; the USS Sabine was "At anchor off Greytown" (Nicaragua) on 2 September, at 23.1 deg MLAT (N11 32', W083 49': Hayakawa et al. 2019, Space Weather 17, 1553), with light "upward to an altitude of about 35°". We run the Sabine's position with the paper's assumed 20 deg and with the logged 35 deg. Endeavour, 1770: lam0 = 20.39 deg S -> paper: 28.23 deg S (§3.2). Eq. 2 (§2.3, Fig. 3a): MLAT = a2 (-Dst)^(2/3) + c, printed as a2 = -0.43, c = 71.53, R^2 = 0.563. Paper's results: 20.40 deg -> -1319 nT; 28.30 deg -> -1025 nT (§3.4). We solve a2 from those pairs. Eq. 3 (§2.3): Gumbel, p = 1 - exp(-exp((mu - x)/sigma)) per solar cycle, T = 11/p years. The paper gives no mu or sigma; we solve them from its central return periods, 10^4.15 years (x = 1319) and 10^3.16 years (x = 1025), and check its ranges 10^(4.15 +1.73 -0.85) and 10^(3.16 +1.15 -0.58) years (§4; exponents from the JATS XML). Dipole latitude of Dongtinghu (29.14 N, 113.12 E, §3.1) for the 1770 pole 80.03 N, 71.54 W (§3.2). Published equatorward boundaries of the auroral oval, in invariant latitude (ILAT): 1859 Sep: 25.1 +/- 0.5 (Hayakawa et al. 2020, EPS 72:122, abstract); 30.8 from the Sabine report (Hayakawa et al. 2018, ApJ 869:57, §3.3; Hayakawa et al. 2019, Space Weather 17, Table 1) 1872 Feb: 24.2 (Hayakawa et al. 2018, ApJ 862:15, abstract; Hayakawa et al. 2023, ApJ 959:23, abstract) 1909 Sep: 31.6 (Hayakawa et al. 2019, Space Weather 17, Table 1) 1921 May: 27.1 (same table; Hayakawa et al. 2023, abstract) Modern checks, boundaries as MLAT: 2024 May 40.99 deg with Dst -412 nT (the paper, §3.3; this is the real-time index, and Kyoto's provisional value is -406 nT, which changes nothing here); 1989 March 40.1 deg (auroral particle precipitation) and 35 deg (auroral electric field), Rich and Denig 1992 as cited by Hayakawa et al. 2019 (Space Weather 17, 1553) and 2023, with Dst -589 nT (WDC Kyoto, final). with storm sizes from magnetograms: 1859 -949 +/- 31 nT (Hayakawa et al. 2022, as cited by Hayakawa et al. 2023); 1872 <= -834 nT at Colaba, whose record has a gap at the peak, and about -1402 +/- 62 nT at Tiflis, a single mid-latitude station under overhead aurora (Hayakawa et al. 2023, §3 and §6); 1909 -595 nT (Love et al. 2019, Space Weather 17, 37); 1921 -907 +/- 132 nT (Love et al. 2019, Space Weather 17, 1281). Conversions between ILAT and the magnetic latitude of a point at height h on the same dipole field line: cos^2(ILAT) = cos^2(lam_h) / (1 + h/a). Run: python3 -I reconstruct.py [out.json] (prints every number; also writes them as JSON if a file name is given) """ import json import math import sys A = 6371.0 def eq1(lam0, beta, h): return lam0 + math.degrees(math.acos(A * math.cos(math.radians(beta)) / (A + h))) - beta def ilat_from(lam_h, h): return math.degrees(math.acos(math.sqrt(math.cos(math.radians(lam_h)) ** 2 / (1 + h / A)))) def lam_from_ilat(ilat, h): return math.degrees(math.acos(math.sqrt(math.cos(math.radians(ilat)) ** 2 * (1 + h / A)))) def dipole_lat(lat, lon, plat, plon): s = (math.sin(math.radians(lat)) * math.sin(math.radians(plat)) + math.cos(math.radians(lat)) * math.cos(math.radians(plat)) * math.cos(math.radians(lon - plon))) return math.degrees(math.asin(s)) def main(): sys.stdout.reconfigure(encoding="utf-8") out = {} # 1. viewing geometry geo = { "Dongtinghu beta60 h240": eq1(19.2, 60, 240), "Dongtinghu beta75 h240": eq1(19.2, 75, 240), "Dongtinghu beta90 h240": eq1(19.2, 90, 240), "Dongtinghu beta60 h100 (green line near 100 km, §3.1)": eq1(19.2, 60, 100), "1859, the paper's point: lam0 20.45 beta20 h400": eq1(20.45, 20, 400), "Endeavour beta20 h400": eq1(20.39, 20, 400), "1859, Sabine's position, assumed 20: lam0 23.1 beta20 h400": eq1(23.1, 20, 400), "1859, Sabine's log (Hayakawa 2018): lam0 23.1 beta35 h400": eq1(23.1, 35, 400), } out["eq1"] = geo print("Eq. 1, boundary latitude (deg):") for k, v in geo.items(): print(f" {k:60s} {v:6.2f}") out["dongtinghu_dipole_lat_1770"] = dipole_lat(29.14, 113.12, 80.03, -71.54) print(f"Dipole latitude of Dongtinghu with the 1770 pole: {out['dongtinghu_dipole_lat_1770']:.2f} (paper: 19.2)") # 2. the regression c = 71.53 a_q = (c - 20.40) / 1319 ** (2 / 3) a_c = (c - 28.30) / 1025 ** (2 / 3) a2 = (a_q + a_c) / 2 def dst(mlat, a=a2): return ((c - mlat) / a) ** 1.5 out["eq2"] = {"a_from_1770": a_q, "a_from_1859": a_c, "a_used": a2, "with_printed_a_0.43": {"1770": dst(20.40, 0.43), "1859": dst(28.30, 0.43)}, "with_solved_a": {"1770": dst(20.40), "1859": dst(28.30)}} print(f"Eq. 2: a solved from 1770 pair {a_q:.5f}, from 1859 pair {a_c:.5f}; printed -0.43") print(f" with printed a = 0.43: 1770 {dst(20.40, 0.43):.0f} nT, 1859 {dst(28.30, 0.43):.0f} nT; ratio {dst(20.40, 0.43) / dst(28.30, 0.43):.3f}") print(f" with solved a: 1770 {dst(20.40):.0f} nT, 1859 {dst(28.30):.0f} nT; ratio {dst(20.40) / dst(28.30):.3f}") print(f" paper's ratio 1319/1025 = {1319 / 1025:.3f}") # 3. Gumbel L = math.log10(11.0) k = (4.15 - 3.16) / (1319 - 1025) # log10 years per nT, tail approximation s_ln10 = 1 / k sigma = s_ln10 / math.log(10) mu = 1025 - (3.16 - L) * s_ln10 def T(x): p = 1 - math.exp(-math.exp((mu - x) / sigma)) return 11.0 / p gum = {"mu": mu, "sigma": sigma, "log10T": {str(x): math.log10(T(x)) for x in (1832, 1319, 1066, 1367, 1025, 853, 900, 850, 589)}} out["gumbel"] = gum print(f"Gumbel solved from the paper's central values: mu = {mu:.0f} nT, sigma = {sigma:.0f} nT per cycle") for x, pub in ((1832, "4.15+1.73 = 5.88"), (1066, "4.15-0.85 = 3.30"), (1367, "3.16+1.15 = 4.31"), (853, "3.16-0.58 = 2.58")): print(f" -Dst {x}: log10 T = {math.log10(T(x)):.2f} (paper: {pub})") for x in (1025, 900, 850, 589): t = T(x) print(f" -Dst {x}: T = {t:,.0f} years; chance in a decade {100 * (1 - math.exp(-10 / t)):.2f}%") # the paper's return periods and ranges in years (its exponents, §4) yrs = {k: 10 ** v for k, v in (("1770 central", 4.15), ("1770 short end", 4.15 - 0.85), ("1770 long end", 4.15 + 1.73), ("1859 central", 3.16), ("1859 short end", 3.16 - 0.58), ("1859 long end", 3.16 + 1.15))} gum["paper_return_periods_years"] = yrs print(" the paper's return periods in years: " + "; ".join(f"{k} {v:,.0f}" for k, v in yrs.items())) gum["chance_per_decade_pct"] = {str(x): 100 * (1 - math.exp(-10 / T(x))) for x in (1025, 900, 850, 589)} # How likely is the paper's own record under its own model? Cycles 2-25 (1766-2025): 24 cycle maxima. n = 24 p_q = 1 - math.exp(-math.exp((mu - 1319) / sigma)) p_c = 1 - math.exp(-math.exp((mu - 1025) / sigma)) p_both = (1 - (1 - p_q) ** n) * (1 - (1 - p_c) ** (n - 1)) # at the short ends of the paper's ranges pq_lo = 11 / 10 ** 3.30 pc_lo = 11 / 10 ** 2.58 p_both_lo = (1 - (1 - pq_lo) ** n) * (1 - (1 - pc_lo) ** (n - 1)) out["own_record_probability"] = {"cycles": n, "p_cycle_ge_1319": p_q, "p_cycle_ge_1025": p_c, "p_both_central": p_both, "p_both_short_ends": p_both_lo} print(f"Chance of a record like 1770+1859 (one cycle max >= 1319 nT and another >= 1025 nT in {n} cycles):" f" {100 * p_both:.2f}% with the central model; {100 * p_both_lo:.1f}% with the short ends of the paper's ranges") # 4. frames: the paper's own boundaries in ILAT fr = {"1770 Yuan (20.40 at 240 km)": ilat_from(20.40, 240), "1770 zenith (19.2) at 240 km": ilat_from(19.2, 240), "1770 zenith (19.2) at 100 km": ilat_from(19.2, 100), "1859 Yuan (28.30 at 400 km)": ilat_from(28.30, 400)} out["ilat_of_paper_boundaries"] = fr print("The paper's boundaries as invariant latitude:") for kk, v in fr.items(): print(f" {kk:32s} {v:5.2f} deg ILAT") # 5. sensitivity: published boundaries through the paper's own regression, in two frames pubs = {"1859 Hayakawa 2020": (25.1, -949), "1859 lower (24.6)": (24.6, -949), "1859 upper (25.6)": (25.6, -949), "1859 Sabine 30.8": (30.8, -949), "1872": (24.2, -834), "1909": (31.6, -595), "1921": (27.1, -907)} sens = {} q240 = dst(20.40) q400 = dst(lam_from_ilat(ilat_from(20.40, 240), 400)) print("Published boundaries through the paper's regression (MLAT at 240 km | at 400 km), against magnetogram values:") print(f" 1770 (Yuan's boundary): {q240:.0f} nT | {q400:.0f} nT") for kk, (il, mag) in pubs.items(): l240, l400 = lam_from_ilat(il, 240), lam_from_ilat(il, 400) d240, d400 = dst(l240), dst(l400) sens[kk] = {"ilat": il, "mlat240": l240, "mlat400": l400, "dst240": d240, "dst400": d400, "magnetogram": mag, "ratio_1770_to_this_240": q240 / d240, "ratio_1770_to_this_400": q400 / d400} print(f" {kk:20s} ILAT {il:4.1f}: MLAT {l240:5.2f}|{l400:5.2f} -> {d240:5.0f}|{d400:5.0f} nT;" f" magnetogram {mag} nT; 1770 is {100 * (q240 / d240 - 1):+.0f}% | {100 * (q400 / d400 - 1):+.0f}%") t_lo, t_hi = 1402 - 62, 1402 + 62 s72 = sorted((sens["1872"]["dst240"], sens["1872"]["dst400"])) print(f" 1872 at Tiflis (one mid-latitude station, under overhead aurora): about -1402 +/- 62 nT, so {t_lo}-{t_hi};" f" the fit's {s72[0]:.0f}-{s72[1]:.0f} lies between Colaba's lower limit and Tiflis") sens["1872"]["tiflis"] = [-1402, 62] out["sensitivity"] = sens out["dst_1770_frames"] = {"240km": q240, "400km": q400} # 5b. the 1859 input: the paper's point and the two readings of the Sabine's log, each in both frames print("The 1859 input through the paper's regression (edges at 400 km; 1770 in the same frame, or as the paper" " compares them, 1770 at 240 km against 1859 at 400 km):") in59 = {"the paper's point (Panama's 20.45, assumed 20 deg)": eq1(20.45, 20, 400), "the Sabine's position (23.1) with the paper's assumed 20 deg": eq1(23.1, 20, 400), "the Sabine's log (23.1, light to about 35 deg)": eq1(23.1, 35, 400)} out["input_1859"] = {} for kk, e400 in in59.items(): il = ilat_from(e400, 400) e240 = lam_from_ilat(il, 240) d240, d400 = dst(e240), dst(e400) m240, m400, mixed = 100 * (q240 / d240 - 1), 100 * (q400 / d400 - 1), 100 * (q240 / d400 - 1) out["input_1859"][kk] = {"mlat400": e400, "ilat": il, "mlat240": e240, "dst240": d240, "dst400": d400, "margin_pct_240": m240, "margin_pct_400": m400, "margin_pct_paper_style": mixed} print(f" {kk:60s} edge {e400:5.2f} at 400 km = {il:5.2f} ILAT -> {d240:5.0f}|{d400:5.0f} nT;" f" 1770 stronger by {m240:.0f}% | {m400:.0f}% (paper's style: {mixed:.0f}%)") # 5c. how high the 1770 light stood (the paper allows 60-90 deg, §2.2) against the 2020 edge for 1859 print("1770's margin over 1859 at 25.1 +/- 0.5 deg ILAT, by the elevation of the Dongtinghu light (both frames):") out["margin_by_1770_elevation"] = {} for beta, e in ((60, 20.40), (75, 19.8), (90, 19.2)): # the paper's edges for 60/75/90 deg (§2.2) il = ilat_from(e, 240) qa, qb = dst(e), dst(lam_from_ilat(il, 400)) ms = [100 * (q / dst(lam_from_ilat(i59, h)) - 1) for i59 in (24.6, 25.1, 25.6) for h, q in ((240, qa), (400, qb))] out["margin_by_1770_elevation"][str(beta)] = {"edge240": e, "ilat": il, "dst240": qa, "dst400": qb, "margin_pct_min": min(ms), "margin_pct_max": max(ms)} print(f" light {beta} deg up: edge {e:5.2f} ({il:5.2f} ILAT), {qa:.0f}|{qb:.0f} nT; 1770 stronger by" f" {min(ms):.0f}-{max(ms):.0f}%") # 6. the paper's own modern check, and 1989, through the same regression (boundaries given as MLAT) modern = {"2024 May, the paper's boundary 40.99 deg (§3.3) vs Dst -412 nT (§3.3)": (40.99, -412), "1989 March, precipitation boundary 40.1 deg MLAT (Rich and Denig 1992, via Hayakawa et al. 2019, 2023) vs -589 nT": (40.1, -589), "1989 March, auroral electric field 35 deg MLAT (same sources) vs -589 nT": (35.0, -589)} out["modern_checks"] = {} print("Modern storms through the paper's regression:") for kk, (lat, mag) in modern.items(): d = dst(lat) fit_lat = c - a2 * abs(mag) ** (2 / 3) out["modern_checks"][kk] = {"mlat": lat, "dst_from_fit": d, "measured": mag, "fit_lat_at_measured_dst": fit_lat} print(f" {kk}: fit gives {d:.0f} nT ({100 * (d / abs(mag) - 1):+.0f}%); at {mag} nT the fit puts the boundary at {fit_lat:.1f} deg") if len(sys.argv) > 1: with open(sys.argv[1], "w", encoding="utf-8", newline="\n") as f: json.dump(out, f, indent=1) print("wrote", sys.argv[1]) if __name__ == "__main__": main()