#!/usr/bin/env python3 """Kp, Dst, Florida's magnetic latitude and the auroral oval, behind "S0 News, November 1, 2026". Fact-Checked Science ("Our analysis"). Standard library only. Run it python3 -I aurora_numbers.py [folder] [--until YYYY-MM-DD] [--all-months] --until (default 2026-10-31, the last full day of the Kp record the page used) ignores anything later, so later storms do not change the page's counts. --all-months reads every Kyoto month from January 1957 (about 840 pages), not only the months with Kp 8- or more, those of the table's storms and March 1981. That is the test of the shorter list, which is kind to Kyoto's server: it found one storm at or below -200 nT outside those months (-215 nT on March 5, 1981, so March 1981 was added) and none at or below -250 nT, and both lists now give the same outputs. The folder (default: the script's own) must hold, or will receive, about 7 MB in all (about 11 MB with --all-months): Kp_ap_Ap_SN_F107_since_1932.txt https://kp.gfz.de/app/files/Kp_ap_Ap_SN_F107_since_1932.txt Kp from Geomagnetic Observatory Niemegk, GFZ Helmholtz Centre for Geosciences, CC BY 4.0: Matzka, Stolle, Yamazaki, Bronkalla and Morschhauser (2021), Space Weather, https://doi.org/10.1029/2020SW002641; data https://doi.org/10.5880/Kp.0001. The file grows daily; its newest values are nowcasts. Its sunspot numbers (CC BY-NC 4.0) are not used. dst_YYYYMM.html https://wdc.kugi.kyoto-u.ac.jp/dst_final/YYYYMM/index.html (then dst_provisional, then dst_realtime) Hourly Dst from the World Data Center for Geomagnetism, Kyoto, https://doi.org/10.17593/14515-74000, under the WDC's rules (https://wdc.kugi.kyoto-u.ac.jp/wdc/Sec3.html): scientific use, no commercial use, with acknowledgement of the WDC and the observatories that supply the data. Every output below that holds Dst carries that condition. igrf14coeffs.txt https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt (IAGA V-MOD; 1900-2025, SV to 2030) COV-OBS.x2-int https://www.spacecenter.dk/files/magnetic-models/COV-OBSx2/COV-OBS.x2-int (Huder et al. 2020, 1840-2020; used only for 1859, before IGRF begins) Missing files are downloaded; the two field models are checked against the SHA-256 of the copies we used. What it writes, beside what it prints aurora-numbers.json every number of ours on the page, with Kyoto's condition in its "terms" key storms-by-cycle.csv storms per solar cycle by Kp and by Dst (Kyoto's condition applies) low-latitude-aurora.csv the table of low-latitude auroras: Kp and Dst from the files, reports typed from the page's sources (Kyoto's condition applies) florida-magnetic-latitude.csv Florida's centred-dipole and apex latitude, 1900-2025 every five years, and at named storms Rules, fixed before running - A storm is a run of UT days whose highest three-hour Kp reaches the threshold, gaps of at most one day merged, as on /s0-news-2026-10-28/. Cycle starts as there (WDC-SILSO minima). Kp thirds: 8.667 = 9-, 9.000 = 9o. - Reading A counts only 9o as Kp 9 (NOAA's G5, and the episode's chart); Reading B counts 9- as well. - Dst storms: hours at or below -250 nT (the usual superstorm line) or -200 nT, joined into one storm when no more than two calendar days apart; each storm's lowest hour is listed. Kyoto's Dst begins in January 1957, so cycles 17 and 18 have none. - In the aurora table, Kp and Ap are searched over the listed UT days; Dst over those days and the next. - Centred-dipole latitude: from the epoch's dipole axis (g10, g11, h11). Apex latitude: the field line of the full model traced to its apex, lat = acos(sqrt(a / r_apex)). IGRF coefficients are interpolated linearly. Geodetic latitudes are converted to geocentric; a = 6371.2 km. - Oval scalings: cos(Lambda) ~ M^(alpha - 1/2), alpha = 1/3 (Siscoe and Chen), 0 (Siebert), 1/2 (Vogt and Glassmeier), from Glassmeier et al. (2004), p. 3669 and Table 1, written for the polar cap's boundary and applied to the oval's edge. - Storm strength: Love et al. (2025), ln L = 3.9825 - 0.5478 ln D (their ln a and b, with the signs their worked examples imply), L = (1 + h/r) / cos^2(lat), h = 250 km, r = 6371 km; the lowest latitude of discrete aurora seen overhead. Schulz (1997): the oval's inside diameter grows 4 +/- 1 degrees per -100 nT of Dst, so about 2 degrees of edge shift. - Line of sight: aurora at height h clears the horizon out to acos(a / (a + h)) of arc, with no refraction. - One standard for the Florida records. For each storm with Florida named: the highest Kp after dark in Florida, counting a three-hour interval on the storm's UT days and the next as dark if the Sun is at least 6 degrees below the horizon at Orlando at its start, middle or end (low-precision solar formula); and the arc from each named place, or from Jacksonville where no place is named, to Love et al.'s edge at the storm's lowest Dst. - Cycle 23 against cycle 25: storms at Kp 9- or more and at Dst -250 nT or below, those with Florida named, over the whole of cycle 23 and over its first stretch as long as cycle 25 has run. - Ebihara and Tanaka's about 2 degrees for a dipole 33% weaker (as Love et al. cite it) is scaled to other changes by the logarithm of the moment, as a power law would scale. """ import csv import datetime as dt import hashlib import json import math import os import re import sys import time import urllib.error import urllib.request A_KM = 6371.2 MU0 = 4e-7 * math.pi UA = "factchecked.science aurora_numbers.py" KP_URL = "https://kp.gfz.de/app/files/Kp_ap_Ap_SN_F107_since_1932.txt" MODELS = { "igrf14coeffs.txt": ("https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf14coeffs.txt", "8f8d88403028fc4ee92c4f38d97b46e0a87e2cfc496045b43c9e26c1d6b0903c"), "COV-OBS.x2-int": ("https://www.spacecenter.dk/files/magnetic-models/COV-OBSx2/COV-OBS.x2-int", "6c588becf4e3896b07be80c3f24697d299432120a7f4ab1139a914395f9d2a17"), } CYCLES = [(17, dt.date(1933, 9, 1)), (18, dt.date(1944, 2, 1)), (19, dt.date(1954, 4, 1)), (20, dt.date(1964, 10, 1)), (21, dt.date(1976, 3, 1)), (22, dt.date(1986, 9, 1)), (23, dt.date(1996, 8, 1)), (24, dt.date(2008, 12, 1)), (25, dt.date(2019, 12, 1))] DST_START = dt.date(1957, 1, 1) EXTRA_DST_MONTHS = ["198103"] # found by --all-months: -215 nT on 1981-03-05, in a month with no Kp 8- # The storms named in the record. days; Dst before 1957 from the literature; Florida named? (yes, no, lead, inferred); reports; sources AURORAE = [ (["1938-01-21", "1938-01-22"], "Dxt -340; Dcx -328", "no", "Overhead at Los Angeles (40.9 deg) with a glow overhead at Escondido (40.1); reported overhead at Noto, Sakhalin (39.2 MLAT), nearly so on Love et al.'s reading; seen from Morioka (29.5 MLAT) up to about 20 deg elevation", "Love 2025 Table 1; Hayakawa 2021 sec. 4"), (["1938-01-25", "1938-01-26"], "Dxt -350; Dcx -336", "no", "Gibraltar, Sicily, Greece; up to 45 deg elevation at Patras (37.3 MLAT); vessel off South Africa (31.8 S, -29.9 MLAT)", "Hayakawa 2021 sec. 3-4"), (["1938-04-16"], "Dxt -255", "no", "Overhead at Long Beach, California (40.5 deg)", "Love 2025 Table 1"), (["1940-03-24", "1940-03-25"], "Dst* -389", "no", "Overhead at Case Observatory, Ohio (52.5 deg)", "Love 2025 Table 1"), (["1941-03-01"], "Dst* <= -464", "no", "Overhead at Kalgoorlie, Australia (30.7 S, 41.3 deg); oval's edge about 38.5 ILAT", "Love 2025 Table 1; Hayakawa 2023 sec. 7"), (["1941-07-04", "1941-07-05"], "Dxt -451", "no", "Overhead at Gaviota, California (40.9 deg)", "Love 2025 Table 1"), (["1941-09-18", "1941-09-19"], "Dxt -415", "no", "Overhead at Riverina, Australia (43.9 deg); New Mexico", "Love 2025 Table 1; Love and Coisson 2016"), (["1946-03-28"], "Dst* <= -512", "no", "Overhead at Broken Hill, Australia (32.0 S, 41.5 deg)", "Love 2025 Table 1"), (["1957-03-01", "1957-03-02"], "", "no", "Rays to Mochita (32.3 MLAT); red glow to Kuji (30.1 MLAT), part perhaps a SAR arc; oval's edge 40.6 ILAT (35.9 if all of Kuji's glow was aurora)", "Hayakawa 2023 sec. 3"), (["1957-09-13"], "", "no", "Seen from Mori (31.9 MLAT); overhead at Fresno (43.3 deg); oval's edge 38.3 ILAT", "Hayakawa 2023 sec. 4; Love 2025 Table 1"), (["1958-02-10", "1958-02-11"], "", "no", "Mexico City (29.3 MLAT); Ogori (23.3 MLAT); rays at Niigata (27.7 MLAT); oval's edge 33.3 ILAT", "Hayakawa 2023 sec. 5"), (["1959-07-15"], "", "no", "Seen at 27.4 MLAT (China); oval's edge 35.4 ILAT; overhead in Heilungkiang (36.3 deg)", "Hayakawa 2024 Abstract; Love 2025 Table 1"), (["1960-11-13"], "", "no", "Overhead at Columbus, Indiana (50.0 deg)", "Love 2025 Table 1"), (["1967-05-25", "1967-05-26"], "", "no", "Overhead at Greenville (43.8 deg). Love's table says Alabama; its latitude and cited paper (Delta Democrat-Times) fit Greenville, Mississippi (33.4 N)", "Love 2025 Table 1; our note"), (["1972-08-04", "1972-08-05"], "", "no", "Satellite edge 53.1 deg: Kp 9o with Dst only -125 nT", "Love 2025 Table 2"), (["1981-04-12", "1981-04-13"], "", "lead", "Every 1981 storm tops out at Kp 8+. A 2025 Florida newspaper article says aurora reached South Florida in 1981, with no date or primary source", "USA TODAY Network-Florida 2025"), (["1982-07-13", "1982-07-14"], "", "no", "No report found by us", "-"), (["1986-02-08", "1986-02-09"], "", "no", "Satellite edge 40-44 deg", "Love 2025 Table 2"), (["1989-03-13", "1989-03-14"], "", "yes", "Overhead at Kosciusko, Mississippi (43.1 deg); visible to 29 MLAT. Florida: First Coast (Times-Union, as retold 2025); Bill Williams recalls seeing it in Florida in 2003 and 1989", "Love 2025 Table 1; Hayakawa 2019 Table 1; USA TODAY Network-Florida 2025; Spaceweather 27 Mar 2023"), (["1991-11-08", "1991-11-09"], "", "no", "Overhead at Mali Losinj, Croatia (44.2 deg)", "Love 2025 Table 1"), (["2000-04-06", "2000-04-07"], "", "yes", "One Florida photograph found, over Rosemary Hill Observatory, central Florida. NASA's release also says Florida but shares its author with the site; NASA's earlier APOD says as far south as North Carolina, which doesn't conflict", "Spaceweather galleries index; NASA Marshall 2000; APOD 2000-04-10"), (["2001-03-31"], "", "no", "Reached Texas, southern California and Mexico; overhead at Payson, Arizona (41.8 deg)", "Spaceweather 1 Apr 2001; Love 2025 Table 1"), (["2001-10-21", "2001-10-22"], "", "no", "Sensitive all-sky record in northern Japan, at about Florida's magnetic latitude", "Shiokawa et al., JpGU 2005 abstract"), (["2001-11-06"], "", "yes", "Reached Florida, Texas and California, the site says; no Florida photograph among its gallery's 70 or so entries", "Spaceweather 6-8 Nov 2001"), (["2001-11-24"], "", "no", "Reached Arkansas and Texas", "Spaceweather 25-26 Nov 2001"), (["2003-10-29"], "", "yes", "Reached Florida, New Mexico, Arizona and Texas, the site says; NOAA: sightings from California to Houston to Florida (29-30 Oct); no Florida photograph in the site's gallery, one from near Houston (29.8 N)", "Spaceweather 30 Oct 2003; NOAA 2004 p. 26"), (["2003-10-30"], "", "yes", "As above; Kyoto's Dst recovered to -97 nT between the two dips (-353 nT, 01 UT; -383 nT, 23 UT)", "Spaceweather 30 Oct 2003; NOAA 2004 p. 26"), (["2003-11-20"], "", "yes", "Reached Florida and Greece, the site's next-day sentence; no Florida photograph among its gallery's two dozen or so entries; the 9- came in Florida's daylight; overhead at Strathalbyn, Australia (43.7 deg)", "Spaceweather 21-22 Nov 2003; Love 2025 Table 1"), (["2004-11-07", "2004-11-08", "2004-11-09", "2004-11-10"], "", "inferred", "Every US state except Hawaii (index); Alabama and California (9 Nov). Florida only inferred", "Spaceweather galleries index; 8-9 Nov 2004"), (["2005-05-15"], "", "no", "California and Arizona; overhead at Mount Arapiles, Australia (44.9 deg)", "Spaceweather 16 May 2005; Love 2025 Table 1"), (["2017-09-07", "2017-09-08"], "", "no", "Reached Arkansas", "Spaceweather 9 Sep 2017"), (["2023-02-26", "2023-02-27"], "", "no", "All-sky imager at Rikubetsu (34.7 MLAT)", "Manu et al. 2025, preprint"), (["2023-03-23", "2023-03-24"], "", "yes", "One 26-minute exposure, Chiefland (29.4 N): a red glow the site judged most likely aurora. New Mexico (32.8 N) on video", "Spaceweather 26-27 Mar 2023; The Watchers 2023"), (["2023-04-23", "2023-04-24"], "", "no", "Terlingua, Texas (29.5 N), camera; southern California; Spain's coast (37 N); Hanle, India, all-sky camera, though the authors find the aurora itself stayed far to the north", "Spaceweather 24-25 Apr 2023; Vichare 2024"), (["2023-11-05"], "", "no", "Texas (33.6 N), Arizona; red SAR arcs", "Spaceweather 6-7 Nov 2023"), (["2024-05-10", "2024-05-11"], "", "yes", "Big Pine Key (24.7 N), by smartphone; Puerto Rico (18.1 N); Mazatlan (23.2 N); naked eye to -23.2 MLAT, oval's edge 29.8-35.5 ILAT", "Spaceweather 11 May 2024; Hayakawa 2025 pp. 34-35"), (["2024-08-12"], "", "no", "Reached Texas, Mississippi and Arizona", "Spaceweather 13 Aug 2024"), (["2024-10-10", "2024-10-11"], "", "yes", "Key West (24.7 N), photo; Lake Mary; Mazatlan; Toa Baja, Puerto Rico (18.4 N); ship south of Cuba (about 21 N)", "Spaceweather 11-12 Oct 2024; ClickOrlando 2024"), (["2025-01-01"], "", "no", "Mexico (Mexicali), Colorado and Arizona; the site's next-day page files the Mexicali photo under SAR images", "Spaceweather 2-3 Jan 2025"), (["2025-06-01"], "", "no", "Tres Piedras, New Mexico (36.6 N); a doubtful photo from Louisiana", "Spaceweather 2 Jun 2025"), (["2025-11-11", "2025-11-12"], "", "yes", "New Port Richey (28 N), Tucson (32 N), Usulutan, El Salvador (13 N), all by naked eye per the site; Shired Island, Florida, camera", "Spaceweather 13 Nov 2025; APOD 2025-11-14; Shaw 2026"), (["2026-01-19", "2026-01-20", "2026-01-21"], "", "yes", "Alabama, Georgia, New Mexico, Arizona and California; press: first night did not reach Florida; near Key West, one photographer's iPhone 15 photos, 10 s by the photographer's note, night of 20 Jan, a red glow", "Spaceweather 20-21 Jan 2026 and gallery item 230497; USA TODAY Network-Florida 20 Jan 2026"), (["2026-03-20", "2026-03-21", "2026-03-22"], "", "no", "Arizona", "Spaceweather 23 Mar 2026"), ] FIVE = ["2023-03-23", "2024-05-10", "2024-10-10", "2025-11-11", "2026-01-19"] # Storms with Florida named, as one row each: first day, last day, epoch, the places the records name (latitude, longitude) JAX, KW = ("Jacksonville", 30.33, -81.66), ("Key West", 24.56, -81.78) FLORIDA_RECORDS = [ ("1989-03-13", "1989-03-14", 1989.2, [JAX]), ("2000-04-06", "2000-04-07", 2000.26, [("Rosemary Hill Observatory", 29.40, -82.59)]), ("2001-11-06", "2001-11-06", 2001.85, [JAX]), ("2003-10-29", "2003-10-30", 2003.83, [JAX]), ("2003-11-20", "2003-11-20", 2003.89, [JAX]), ("2023-03-23", "2023-03-24", 2023.22, [("Chiefland", 29.47, -82.86)]), ("2024-05-10", "2024-05-11", 2024.36, [("Tallahassee", 30.44, -84.28), ("Big Pine Key", 24.67, -81.35)]), ("2024-10-10", "2024-10-11", 2024.78, [("Lake Mary", 28.76, -81.32), KW]), ("2025-11-11", "2025-11-12", 2025.86, [("Shired Island", 29.40, -83.20), ("New Port Richey", 28.24, -82.72)]), ("2026-01-19", "2026-01-21", 2026.05, [KW]), ] FLORIDA = [("Key West", 24.56, -81.78), ("Miami", 25.76, -80.19), ("New Port Richey", 28.24, -82.72), ("Orlando", 28.54, -81.38), ("Chiefland", 29.47, -82.86), ("Jacksonville", 30.33, -81.66), ("Tallahassee", 30.44, -84.28), ("Pensacola", 30.42, -87.22)] NAMED_EPOCHS = [1859.67, 1921.37, 1960.0, 1989.2, 2003.83, 2024.36, 2025.0] OTHERS = [ # the place, geodetic latitude, longitude, epoch, Love et al.'s geomagnetic latitude where Table 1 gives one ("Terlingua, TX (Apr 2023)", 29.33, -103.62, 2023.31, None), ("Sunspot, NM (Mar 2023)", 32.79, -105.82, 2023.22, None), ("Tucson, AZ (Nov 2025)", 32.22, -110.97, 2025.86, None), ("Usulutan, El Salvador (Nov 2025)", 13.35, -88.44, 2025.86, None), ("Hanle, India (Apr 2023)", 32.78, 78.96, 2023.31, None), ("Rikubetsu, Japan (Apr 2023)", 43.46, 143.77, 2023.31, None), ("Rikubetsu, Japan, 1960", 43.46, 143.77, 1960.0, None), ("Moshiri, Japan (2002)", 44.37, 142.27, 2002.0, None), ("Los Angeles (Jan 1938)", 34.05, -118.24, 1938.06, 40.86), ("Escondido (Jan 1938)", 33.12, -117.09, 1938.06, 40.09), ("Fresno (Feb 1958)", 36.74, -119.77, 1958.11, 43.32), ("Kosciusko, MS (Mar 1989)", 33.06, -89.59, 1989.2, 43.10), ("Payson, AZ (Mar 2001)", 34.23, -111.33, 2001.25, 41.79), ("San Juan, PR (May 1921)", 18.47, -66.11, 1921.37, 29.73), ("Greenville, MS (May 1967)", 33.41, -91.06, 1967.4, 43.79), ("Greenville, AL (May 1967)", 31.83, -86.63, 1967.4, 43.79), ("Houston, TX (Oct 2003)", 29.76, -95.37, 2003.83, None), ] def say(*a): print(*a) def cycle(d): c = 16 for n, s in CYCLES: if d >= s: c = n return c def kpname(k): i = int(round(k * 3)) base = (i + 1) // 3 third = i - 3 * base return f"{base}{'-' if third == -1 else '+' if third == 1 else 'o'}" def get(url): req = urllib.request.Request(url, headers={"User-Agent": UA}) with urllib.request.urlopen(req, timeout=300) as r: return r.read() def fetch_file(folder, name, url, want=None): path = os.path.join(folder, name) if not os.path.exists(path): say("downloading", url) with open(path, "wb") as f: f.write(get(url)) if want: with open(path, "rb") as f: if hashlib.sha256(f.read()).hexdigest() != want: say(f"Note: {name} differs from the copy we used; results may differ slightly.") return path def fetch_dst(folder, ym): path = os.path.join(folder, f"dst_{ym}.html") if os.path.exists(path): return path for kind in ("final", "provisional", "realtime"): try: data = get(f"https://wdc.kugi.kyoto-u.ac.jp/dst_{kind}/{ym}/index.html") except urllib.error.HTTPError: continue with open(path, "wb") as f: f.write(data) time.sleep(0.25) return path return None def read_kp(path, until): rows = {} for line in open(path, encoding="utf-8", errors="replace"): if line.startswith("#") or not line.strip(): continue f = line.split() day = dt.date(int(f[0]), int(f[1]), int(f[2])) kps = [float(x) for x in f[7:15]] if day <= until and max(kps) >= 0: rows[day] = {"kps": kps, "max": max(kps), "ap": int(f[23]), "d": int(f[-1])} return rows def read_dst(path): txt = re.sub(r"<[^>]+>", "", open(path, encoding="utf-8", errors="replace").read()) m = re.search(r"\n\s+([A-Z]+)\s+(\d{4})\s*\n", txt) month = dt.datetime.strptime(m.group(1).title() + " " + m.group(2), "%B %Y").date() days = {} for line in txt.splitlines(): mm = re.match(r"^\s*(\d{1,2})\s+(-?\d.*)$", line) if mm: vals = [int(v) for v in re.findall(r"-?\d+", mm.group(2).replace("-", " -"))] if len(vals) >= 24: days[month.replace(day=int(mm.group(1)))] = [v for v in vals[:24]] return days def storms(kp, thr): out = [] for d in sorted(kp): if kp[d]["max"] >= thr - 1e-3: if out and (d - out[-1][-1]).days <= 1: out[-1].append(d) else: out.append([d]) return out def sun_altitude(lat, lon, when): """The Sun's altitude in degrees, by the usual low-precision formula.""" n = (when - dt.datetime(2000, 1, 1, 12)).total_seconds() / 86400.0 g = math.radians((357.529 + 0.98560028 * n) % 360) q = (280.459 + 0.98564736 * n) % 360 ecl = math.radians(q + 1.915 * math.sin(g) + 0.020 * math.sin(2 * g)) obl = math.radians(23.439 - 0.00000036 * n) dec = math.asin(math.sin(obl) * math.sin(ecl)) ra = math.degrees(math.atan2(math.cos(obl) * math.sin(ecl), math.cos(ecl))) ha = math.radians((280.46061837 + 360.98564736629 * n) + lon - ra) la = math.radians(lat) return math.degrees(math.asin(math.sin(la) * math.sin(dec) + math.cos(la) * math.cos(dec) * math.cos(ha))) def kp_after_dark(kp, first, last): """The highest Kp in three-hour intervals, on the storm's UT days and the next, with the Sun 6 degrees or more down at Orlando.""" best, when = None, None d = first while d <= last + dt.timedelta(days=1): if d in kp: for i, v in enumerate(kp[d]["kps"]): start = dt.datetime(d.year, d.month, d.day) + dt.timedelta(hours=3 * i) dark = any(sun_altitude(28.54, -81.38, start + dt.timedelta(minutes=m)) <= -6 for m in (0, 90, 180)) if dark and (best is None or v > best + 1e-6): best, when = v, start d += dt.timedelta(days=1) return best, when # ---------------------------------------------------------------- field models def read_igrf(path): header, rows = None, [] for line in open(path, encoding="utf-8", errors="replace"): p = line.split() if not p or line.startswith("#"): continue if p[0] == "g/h": header = p elif p[0] in ("g", "h") and header is not None: rows.append(p) cols = header[3:] epochs = {float(c): {} for c in cols if "-" not in c} sv = {} for p in rows: key = (p[0], int(p[1]), int(p[2])) for c, v in zip(cols, p[3:]): if "-" in c: sv[key] = float(v) else: epochs[float(c)][key] = float(v) last = max(epochs) epochs[last + 5.0] = {k: v + 5.0 * sv.get(k, 0.0) for k, v in epochs[last].items()} # IGRF's secular variation, 2025-2030 return epochs def igrf_at(epochs, t): ks = sorted(epochs) for a, b in zip(ks, ks[1:]): if a <= t <= b: w = (t - a) / (b - a) keys = set(epochs[a]) | set(epochs[b]) return {k: (1 - w) * epochs[a].get(k, 0.0) + w * epochs[b].get(k, 0.0) for k in keys} raise ValueError(t) class CovObs: def __init__(self, path): lines = open(path, encoding="ascii").read().split("\n") head = lines[1].split() self.nmax, self.nspl, order = int(head[0]), int(head[1]), int(head[2]) self.knots = [float(x) for x in head[3:]] vals = [float(x) for line in lines[2:] for x in line.split()] n = self.nmax * (self.nmax + 2) self.coefs = [vals[j * n:(j + 1) * n] for j in range(self.nspl)] self.p = order - 1 def at(self, t): U, p = self.knots, self.p i = p while not (U[i] <= t < U[i + 1]): i += 1 N = [1.0] + [0.0] * p left, right = [0.0] * (p + 1), [0.0] * (p + 1) for j in range(1, p + 1): left[j], right[j] = t - U[i + 1 - j], U[i + j] - t saved = 0.0 for r in range(j): temp = N[r] / (right[r + 1] + left[j - r]) N[r] = saved + right[r + 1] * temp saved = left[j - r] * temp N[j] = saved out = [0.0] * len(self.coefs[0]) for j in range(p + 1): row = self.coefs[i - p + j] for q in range(len(out)): out[q] += N[j] * row[q] c = {} for n in range(1, self.nmax + 1): base = (n - 1) * (n + 1) c[("g", n, 0)] = out[base] for m in range(1, n + 1): c[("g", n, m)] = out[base + 2 * m - 1] c[("h", n, m)] = out[base + 2 * m] return c def dipole(c): g10, g11, h11 = c[("g", 1, 0)], c[("g", 1, 1)], c[("h", 1, 1)] b0 = math.sqrt(g10 ** 2 + g11 ** 2 + h11 ** 2) colat = math.degrees(math.acos(-g10 / b0)) return {"M": 4 * math.pi * (A_KM * 1e3) ** 3 / MU0 * b0 * 1e-9 / 1e22, "plat": 90.0 - colat, "plon": math.degrees(math.atan2(-h11, -g11))} def cd_latitude(lat, lon, plat, plon): la, lo, pa, po = map(math.radians, (lat, lon, plat, plon)) return math.degrees(math.asin(max(-1.0, min(1.0, math.sin(la) * math.sin(pa) + math.cos(la) * math.cos(pa) * math.cos(lo - po))))) def legendre(nmax, theta): ct, st = math.cos(theta), math.sin(theta) P = [[0.0] * (nmax + 2) for _ in range(nmax + 2)] dP = [[0.0] * (nmax + 2) for _ in range(nmax + 2)] P[0][0] = 1.0 for n in range(1, nmax + 1): if n == 1: P[1][1], dP[1][1] = st, ct else: k = math.sqrt((2 * n - 1) / (2 * n)) P[n][n] = k * st * P[n - 1][n - 1] dP[n][n] = k * (ct * P[n - 1][n - 1] + st * dP[n - 1][n - 1]) for m in range(0, n): a, b, d = 2 * n - 1, math.sqrt(max(0.0, (n - 1) ** 2 - m ** 2)), math.sqrt(n ** 2 - m ** 2) p2 = P[n - 2][m] if n >= 2 else 0.0 dp2 = dP[n - 2][m] if n >= 2 else 0.0 P[n][m] = (a * ct * P[n - 1][m] - b * p2) / d dP[n][m] = (a * (ct * dP[n - 1][m] - st * P[n - 1][m]) - b * dp2) / d return P, dP def field(c, r_km, colat_deg, lon_deg, nmax): theta, phi = math.radians(colat_deg), math.radians(lon_deg) P, dP = legendre(nmax, theta) st = math.sin(theta) or 1e-10 br = bt = bp = 0.0 ar = A_KM / r_km for n in range(1, nmax + 1): f = ar ** (n + 2) for m in range(0, n + 1): g = c.get(("g", n, m), 0.0) h = c.get(("h", n, m), 0.0) if m > 0 else 0.0 cm, sm = math.cos(m * phi), math.sin(m * phi) gh = g * cm + h * sm br += (n + 1) * f * gh * P[n][m] bt += -f * gh * dP[n][m] bp += f * m * (g * sm - h * cm) * P[n][m] / st return br, bt, bp def bdir(c, x, y, z, nmax): r = math.sqrt(x * x + y * y + z * z) th = math.acos(max(-1.0, min(1.0, z / r))) ph = math.atan2(y, x) br, bt, bp = field(c, r, math.degrees(th), math.degrees(ph), nmax) st, ct, sp, cp = math.sin(th), math.cos(th), math.sin(ph), math.cos(ph) bx = br * st * cp + bt * ct * cp - bp * sp by = br * st * sp + bt * ct * sp + bp * cp bz = br * ct - bt * st b = math.sqrt(bx * bx + by * by + bz * bz) return bx / b, by / b, bz / b, br def apex_latitude(c, lat, lon, step_km=20.0, max_steps=200000): nmax = max(n for (_, n, _) in c) th, ph = math.radians(90.0 - lat), math.radians(lon) x, y, z = A_KM * math.sin(th) * math.cos(ph), A_KM * math.sin(th) * math.sin(ph), A_KM * math.cos(th) sgn = 1.0 if bdir(c, x, y, z, nmax)[3] > 0 else -1.0 rmax = prev_r = A_KM h = step_km for _ in range(max_steps): k1 = bdir(c, x, y, z, nmax)[:3] k2 = bdir(c, x + 0.5 * h * sgn * k1[0], y + 0.5 * h * sgn * k1[1], z + 0.5 * h * sgn * k1[2], nmax)[:3] k3 = bdir(c, x + 0.5 * h * sgn * k2[0], y + 0.5 * h * sgn * k2[1], z + 0.5 * h * sgn * k2[2], nmax)[:3] k4 = bdir(c, x + h * sgn * k3[0], y + h * sgn * k3[1], z + h * sgn * k3[2], nmax)[:3] dx = [h * sgn * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6.0 for i in range(3)] x, y, z = x + dx[0], y + dx[1], z + dx[2] r = math.sqrt(x * x + y * y + z * z) rmax = max(rmax, r) if r < prev_r: if h > 0.5: x, y, z = x - dx[0], y - dx[1], z - dx[2] h /= 4.0 prev_r = math.sqrt(x * x + y * y + z * z) continue break prev_r = r h = max(step_km, 0.01 * r) if h >= step_km else h if r > 60 * A_KM: break lam = math.degrees(math.acos(math.sqrt(min(1.0, A_KM / rmax)))) return lam if lat >= 0 else -lam def geocentric(lat_deg): a, f = 6378.137, 1 / 298.257223563 la = math.radians(lat_deg) n = a / math.sqrt(1 - (2 * f - f * f) * math.sin(la) ** 2) return math.degrees(math.atan2(n * (1 - f) ** 2 * math.sin(la), n * math.cos(la))) def love_lat(D, shift=0.0): lz = math.exp(3.9825 - 0.5478 * math.log(D) + shift) v = (1 + 250.0 / 6371.0) / lz return math.degrees(math.acos(math.sqrt(min(1.0, v)))) def main(): sys.stdout.reconfigure(encoding="utf-8") args = sys.argv[1:] until = dt.date(2026, 10, 31) if "--until" in args: i = args.index("--until") until = dt.date.fromisoformat(args[i + 1]) del args[i:i + 2] all_months = "--all-months" in args args = [a for a in args if a != "--all-months"] folder = args[0] if args else os.path.dirname(os.path.abspath(__file__)) os.makedirs(folder, exist_ok=True) kp = read_kp(fetch_file(folder, "Kp_ap_Ap_SN_F107_since_1932.txt", KP_URL), until) igrf = read_igrf(fetch_file(folder, "igrf14coeffs.txt", *MODELS["igrf14coeffs.txt"])) cov = CovObs(fetch_file(folder, "COV-OBS.x2-int", *MODELS["COV-OBS.x2-int"])) out = {"terms": "Kp and Ap from GFZ, CC BY 4.0. Dst from WDC Kyoto, doi:10.17593/14515-74000: scientific use, no commercial use; we acknowledge the WDC and the observatories that supply the data. Field values from IGRF-14 and COV-OBS.x2.", "until": str(until), "last_definitive_kp_day": str(max(d for d, r in kp.items() if r["d"] >= 1))} say(f"Kp record to {max(kp)}; definitive to {out['last_definitive_kp_day']}.") # ---------------------------------------------- A. storms per cycle by Kp by = {c: {"cycle": c} for c, _ in CYCLES} for i, (c, s) in enumerate(CYCLES): end = CYCLES[i + 1][1] if i + 1 < len(CYCLES) else until + dt.timedelta(days=1) by[c]["years"] = round((end - s).days / 365.25, 1) by[c]["start"], by[c]["dst_years"] = str(s), round((end - max(s, DST_START)).days / 365.25, 1) for name, thr in (("kp_8m", 7.667), ("kp_8p", 8.333), ("kp_9m", 8.667), ("kp_9o", 9.0)): for c in by: by[c][name] = 0 for st in storms(kp, thr): if cycle(st[0]) in by: by[cycle(st[0])][name] += 1 out["kp_9m_storms"] = [{"from": str(st[0]), "to": str(st[-1]), "cycle": cycle(st[0]), "max": kpname(max(kp[d]["max"] for d in st)), "ap": max(kp[d]["ap"] for d in st)} for st in storms(kp, 8.667)] out["kp_9o_storms_total"] = len(storms(kp, 9.0)) out["kp_9m_storms_total"] = len(storms(kp, 8.667)) say("Cycle, years, storms at Kp 8-, 8+, 9-, 9o:") for c in by: say(f" {c}: {by[c]['years']}, {by[c]['kp_8m']}, {by[c]['kp_8p']}, {by[c]['kp_9m']}, {by[c]['kp_9o']}") # ---------------------------------------------- B. Dst months = {d.strftime("%Y%m") for d, r in kp.items() if r["max"] >= 7.667 and d >= DST_START} | set(EXTRA_DST_MONTHS) for days, *_ in AURORAE: for d in days: if dt.date.fromisoformat(d) >= DST_START: months.add(d[:4] + d[5:7]) nxt = dt.date.fromisoformat(d) + dt.timedelta(days=1) months.add(nxt.strftime("%Y%m")) if all_months: y, m = DST_START.year, DST_START.month while dt.date(y, m, 1) <= until: months.add(f"{y}{m:02d}") y, m = (y + 1, 1) if m == 12 else (y, m + 1) dst = {} missing = [] for ym in sorted(months): p = fetch_dst(folder, ym) if p: dst.update(read_dst(p)) else: missing.append(ym) say(f"Dst: {len(months)} months read; none found for {missing}" if missing else f"Dst: {len(months)} months read.") series = [(d, h + 1, v) for d in sorted(dst) if d <= until for h, v in enumerate(dst[d]) if v != 9999] for thr, key in ((-250, "dst_250"), (-200, "dst_200")): groups = [] for d, h, v in series: if v <= thr: if groups and (d - groups[-1][-1][0]).days <= 2: groups[-1].append((d, h, v)) else: groups.append([(d, h, v)]) for c in by: by[c][key] = 0 by[c][key + "_list"] = [] for g in groups: m = min(g, key=lambda x: x[2]) c = cycle(m[0]) by[c][key] += 1 by[c][key + "_list"].append(f"{m[0]} {m[2]}") say(f"Dst storms at or below {thr} nT, by cycle (from January 1957):") for c in by: if c >= 19: say(f" {c} ({by[c]['dst_years']} yr): {by[c][key]}: {', '.join(by[c][key + '_list'])}") for c in by: if c < 19: # before Kyoto's Dst begins by[c].update(dst_years=None, dst_250=None, dst_250_list=None, dst_200=None, dst_200_list=None) out["cycles"] = [by[c] for c in by] with open(os.path.join(folder, "storms-by-cycle.csv"), "w", encoding="utf-8", newline="\n") as f: f.write("# Kp from GFZ, CC BY 4.0. Dst from WDC Kyoto, doi:10.17593/14515-74000: scientific use, no commercial use; we acknowledge the WDC and the observatories that supply the data.\n") w = csv.writer(f, lineterminator="\n") w.writerow(["cycle", "starts", "years_in_kp", "kp_8m_or_more", "kp_8p_or_more", "kp_9m_or_more", "kp_9o", "years_in_dst", "dst_250_or_deeper", "dst_200_or_deeper"]) for c in by: b = by[c] w.writerow([c, b["start"], b["years"], b["kp_8m"], b["kp_8p"], b["kp_9m"], b["kp_9o"], b["dst_years"] if c >= 19 else "", b["dst_250"] if c >= 19 else "", b["dst_200"] if c >= 19 else ""]) # ---------------------------------------------- C and D. the table of auroras and the two readings rows = [] tally = {"A": [0, 0], "B": [0, 0]} five = [] for days, lit, fl, rep, src in AURORAE: ds = [dt.date.fromisoformat(d) for d in days] mk = max(kp[d]["max"] for d in ds) ap = max(kp[d]["ap"] for d in ds) dm = None for d in ds + [ds[-1] + dt.timedelta(days=1)]: if d in dst: vals = [x for x in dst[d] if x != 9999] dm = min(vals) if dm is None else min(dm, min(vals)) a, b = mk >= 9.0 - 1e-3, mk >= 8.667 - 1e-3 label = days[0] if len(days) == 1 else f"{days[0]} to {days[-1]}" rows.append([label, kpname(mk), f"{mk:.3f}", ap, dm if dm is not None else lit, "yes" if a else "no", "yes" if b else "no", fl, rep, src]) if fl == "yes": tally["A"][0 if a else 1] += 1 tally["B"][0 if b else 1] += 1 if days[0] in FIVE: five.append({"storm": label, "kp": kpname(mk), "ap": ap, "dst": dm, "reading_A": a, "reading_B": b}) with open(os.path.join(folder, "low-latitude-aurora.csv"), "w", encoding="utf-8", newline="\n") as f: f.write("# Kp and Ap from GFZ, CC BY 4.0. Dst from WDC Kyoto, doi:10.17593/14515-74000: scientific use, no commercial use; we acknowledge the WDC and the observatories that supply the data. Reports are summarized from the page's sources. Sightings are not a census.\n") w = csv.writer(f, lineterminator="\n") w.writerow(["storm_ut_days", "highest_kp", "highest_kp_value", "highest_daily_ap", "lowest_dst_nT", "kp9_reading_A_9o_only", "kp9_reading_B_9m_or_9o", "florida_named", "reports", "sources"]) w.writerows(rows) out["five"] = five out["five_without_kp9"] = {"reading_A": sum(not x["reading_A"] for x in five), "reading_B": sum(not x["reading_B"] for x in five)} out["florida_rows"] = {"reading_A_reach": tally["A"][0], "reading_A_short": tally["A"][1], "reading_B_reach": tally["B"][0], "reading_B_short": tally["B"][1]} say("The five:", five) say("Without Kp 9:", out["five_without_kp9"], "; rows with Florida named:", out["florida_rows"]) out["jan_2026_kp"] = {str(d): kp[d]["kps"] for d in (dt.date(2026, 1, 19), dt.date(2026, 1, 20), dt.date(2026, 1, 21)) if d in kp} out["kp_1981_max"] = max(r["max"] for d, r in kp.items() if d.year == 1981) oct03 = {str(d): dst[d] for d in (dt.date(2003, 10, 29), dt.date(2003, 10, 30)) if d in dst} out["dst_oct_2003"] = oct03 # ---------------------------------------------- E. Florida def coeffs(t): return cov.at(t) if t < 1900 else igrf_at(igrf, t) dip = {} for t in [1859.67, 1921.37, 1960.0, 1989.2, 2003.83, 2015.0, 2024.36, 2025.0]: d = dipole(coeffs(t)) dip[t] = d out["dipole_1e22"] = {str(t): {"M": round(d["M"], 3), "pole_lat": round(d["plat"], 2), "pole_lon": round(d["plon"], 2)} for t, d in dip.items()} change = lambda a, b: 100 * (dip[b]["M"] / dip[a]["M"] - 1) out["dipole_change_pct"] = {f"{a}-{b}": round(change(a, b), 2) for a, b in [(1859.67, 2025.0), (1960.0, 2025.0), (1989.2, 2025.0), (2003.83, 2024.36), (2015.0, 2025.0)]} out["dipole_rate_pct_per_century"] = {"2003.83-2024.36": round(change(2003.83, 2024.36) / (2024.36 - 2003.83) * 100, 2), "2015-2025": round(change(2015.0, 2025.0) / 10 * 100, 2)} say("Dipole:", out["dipole_change_pct"], out["dipole_rate_pct_per_century"]) lat = {} series5 = [float(y) for y in range(1900, 2030, 5)] for name, la, lo in FLORIDA: gl = geocentric(la) lat[name] = {} for t in sorted(set(series5 + NAMED_EPOCHS)): c = coeffs(t) d = dipole(c) lat[name][t] = (round(cd_latitude(gl, lo, d["plat"], d["plon"]), 2), round(apex_latitude(c, gl, lo), 2)) out["florida"] = {n: {str(t): v for t, v in lat[n].items()} for n in lat} with open(os.path.join(folder, "florida-magnetic-latitude.csv"), "w", encoding="utf-8", newline="\n") as f: f.write("# Our calculation from IGRF-14 (1900 on) and COV-OBS.x2 (1859.67, degree 14). cd = centred-dipole latitude; apex = apex latitude of the traced field line. Degrees.\n") w = csv.writer(f, lineterminator="\n") w.writerow(["epoch"] + [f"{n} {k}" for n in lat for k in ("cd", "apex")]) for t in sorted(set(series5 + NAMED_EPOCHS)): w.writerow([t] + [x for n in lat for x in lat[n][t]]) say("Florida, epoch: Key West, Orlando, Pensacola (cd/apex)") for t in NAMED_EPOCHS: say(f" {t}: " + "; ".join(f"{n} {lat[n][t][0]}/{lat[n][t][1]}" for n in ("Key West", "Orlando", "Pensacola"))) out["others"] = [] for name, la, lo, t, love in OTHERS: c = coeffs(t) d = dipole(c) gl = geocentric(la) rec = {"place": name, "epoch": t, "cd": round(cd_latitude(gl, lo, d["plat"], d["plon"]), 2), "apex": round(apex_latitude(c, gl, lo), 2), "love_table_1": love} out["others"].append(rec) say(" ", rec) band = {} d25 = dipole(igrf_at(igrf, 2025.0)) for lon in (-80, 0, 20, 140, 150): cells = [] for tgt in (36.0, 39.5, -36.0, -39.5): lo_, hi_ = -89.0, 89.0 for _ in range(60): mid = 0.5 * (lo_ + hi_) if cd_latitude(mid, lon, d25["plat"], d25["plon"]) < tgt: lo_ = mid else: hi_ = mid cells.append(round(0.5 * (lo_ + hi_), 1)) band[str(lon)] = cells out["florida_band_2025_cd36_39p5_north_south"] = band say("Band, geocentric latitude of CD 36 and 39.5, north and south, by longitude:", band) # ---------------------------------------------- F. the oval edge = {} for e in (30.0, 36.0, 45.0, 50.0, 55.0): edge[str(e)] = {} for dm in (change(2003.83, 2024.36), -3.0, -5.0, change(1859.67, 2025.0)): r = 1 + dm / 100 vals = [] for a in (1 / 3, 0.0, 0.5): vals.append(round(e - math.degrees(math.acos(math.cos(math.radians(e)) * r ** (a - 0.5))), 2)) edge[str(e)][f"{dm:.2f}"] = {"siscoe_chen": vals[0], "siebert": vals[1], "vogt_glassmeier": vals[2]} out["edge_shift_deg"] = edge say("Edge shift, degrees toward the equator:", edge) love = {str(D): {"lat": round(love_lat(D), 2), "hi": round(love_lat(D, 0.1322), 2), "lo": round(love_lat(D, -0.1322), 2)} for D in (163, 213, 217, 236, 292, 333, 383, 406, 422, 589, 964)} out["love_model"] = love out["love_shift_163_406"] = round(love_lat(163) - love_lat(406), 2) out["schulz_edge_shift"] = {str(D): round(2 * D / 100, 1) for D in (163, 217, 236, 333, 406)} out["schulz_shift_163_406"] = round(2 * (406 - 163) / 100, 1) say("Love's model:", love, "; -163 to -406 nT:", out["love_shift_163_406"], "; Schulz:", out["schulz_shift_163_406"]) need = {} for n in ("Key West", "Orlando", "Pensacola"): lz = (1 + 250.0 / 6371.0) / math.cos(math.radians(lat[n][2025.0][0])) ** 2 need[n] = round(math.exp((3.9825 - math.log(lz)) / 0.5478)) out["love_dst_for_overhead_2025"] = need say("Dst for discrete aurora overhead, Love's model:", need) out["horizon_arc_deg"] = {str(h): round(math.degrees(math.acos(A_KM / (A_KM + h))), 1) for h in (250, 300, 400, 500, 650, 830, 950)} ch = dipole(igrf_at(igrf, 2023.22)) chief_cd = cd_latitude(geocentric(29.47), -82.86, ch["plat"], ch["plon"]) out["chiefland_2023_cd"] = round(chief_cd, 2) out["chiefland_arc_to_love_edge_deg"] = {"mid": round(love_lat(163) - chief_cd, 1), "low": round(love_lat(163, -0.1322) - chief_cd, 1), "high": round(love_lat(163, 0.1322) - chief_cd, 1)} say("Horizon arcs:", out["horizon_arc_deg"], "; Chiefland to Love's edge at -163 nT:", out["chiefland_arc_to_love_edge_deg"]) moment = change(2003.83, 2024.36) / 100 out["ebihara_tanaka_scaled_deg"] = round(2.0 * math.log(1 + moment) / math.log(1 - 0.33), 2) # ---------------------------------------------- G. one standard for the Florida records recs = [] named_days = set() for first, last, epoch, places in FLORIDA_RECORDS: d0, d1 = dt.date.fromisoformat(first), dt.date.fromisoformat(last) days = [d0 + dt.timedelta(days=i) for i in range((d1 - d0).days + 1)] named_days.update(days) vals = [x for d in days + [d1 + dt.timedelta(days=1)] if d in dst for x in dst[d] if x != 9999] low = min(vals) c = coeffs(epoch) di = dipole(c) arcs = {n: round(love_lat(abs(low)) - cd_latitude(geocentric(la), lo, di["plat"], di["plon"]), 1) for n, la, lo in places} best, when = kp_after_dark(kp, d0, d1) recs.append({"first": first, "last": last, "highest_kp": kpname(max(kp[d]["max"] for d in days)), "lowest_dst": low, "love_edge": round(love_lat(abs(low)), 1), "arc_to_edge": arcs, "kp_after_dark": kpname(best), "kp_after_dark_from_utc": when.strftime("%Y-%m-%d %H:%M")}) out["florida_records"] = recs say("Florida records:", recs) # ---------------------------------------------- H. cycle 23 against cycle 25, on the same measures start = dict(CYCLES) span = until + dt.timedelta(days=1) - start[25] cut = start[23] + span def dst_day(entry): return dt.date.fromisoformat(entry.split()[0]) versus = {"same_point_ends": str(cut), "span_years": round(span.days / 365.25, 2)} kp9 = storms(kp, 8.667) for c in (23, 25): st = [x for x in kp9 if cycle(x[0]) == c] dl = [dst_day(e) for e in by[c]["dst_250_list"]] hit = lambda d: d in named_days or d - dt.timedelta(days=1) in named_days fl = [r for r in recs if cycle(dt.date.fromisoformat(r["first"])) == c] versus[str(c)] = {"florida_named": len(fl), "kp_9m": len(st), "kp_9m_named": sum(any(d in named_days for d in x) for x in st), "dst_250": len(dl), "dst_250_named": sum(hit(d) for d in dl)} if c == 23: versus["23_same_point"] = {"florida_named": sum(dt.date.fromisoformat(r["first"]) < cut for r in fl), "kp_9m": sum(x[0] < cut for x in st), "kp_9m_named": sum(x[0] < cut and any(d in named_days for d in x) for x in st), "dst_250": sum(d < cut for d in dl), "dst_250_named": sum(d < cut and hit(d) for d in dl)} out["cycle_23_against_25"] = versus out["chance_all_4_named"] = round(math.comb(8, 4) / math.comb(14, 4), 3) out["florida_named_of_kp_storms"] = {"9o": [sum(any(d in named_days for d in x) for x in storms(kp, 9.0)), out["kp_9o_storms_total"]], "9m_or_more": [sum(any(d in named_days for d in x) for x in kp9), out["kp_9m_storms_total"]]} out["kp_9o_storms_before_cycle_17"] = sum(x[0] < start[17] for x in storms(kp, 9.0)) say("Cycle 23 against 25:", versus, "; chance all 4 named:", out["chance_all_4_named"], "; named of Kp storms:", out["florida_named_of_kp_storms"], "; 9o before cycle 17:", out["kp_9o_storms_before_cycle_17"]) with open(os.path.join(folder, "aurora-numbers.json"), "w", encoding="utf-8", newline="\n") as f: json.dump(out, f, ensure_ascii=False, indent=1) say("Wrote aurora-numbers.json, storms-by-cycle.csv, low-latitude-aurora.csv, florida-magnetic-latitude.csv") if __name__ == "__main__": main()