#!/usr/bin/env python3
"""Distances and time gaps behind our review of S0 News, October 23, 2026.

Great-circle distances (haversine, mean Earth radius 6,371 km) from the Sicily epicentres to Campi Flegrei
and other places, and the time gaps the page quotes. Needs only Python 3:  python3 -I distances.py
"""
import math
from datetime import datetime, timedelta

R = 6371.0  # km


def km(a, b):
    la1, lo1 = map(math.radians, a)
    la2, lo2 = map(math.radians, b)
    h = math.sin((la2 - la1) / 2) ** 2 + math.cos(la1) * math.cos(la2) * math.sin((lo2 - lo1) / 2) ** 2
    return 2 * R * math.asin(math.sqrt(h))


# Epicentres (latitude, longitude): INGV ISIDe and USGS ComCat, October 23, 2026
EPICENTRES = {
    "INGV Mw 5.5, 10:19:15 UTC": (38.7938, 14.1717),
    "USGS M 5.7, 10:19:16 UTC": (38.7199, 14.1520),
    "INGV Mw 4.8, 09:51:38 UTC": (38.7277, 14.1720),
    "USGS M 4.9, 09:51:39 UTC": (38.7219, 14.2527),
}
# Reference points (approximate coordinates of the named places)
PLACES = {
    "Alicudi": (38.54, 14.35),
    "Salina (Leni)": (38.56, 14.83),
    "Stromboli": (38.79, 15.21),
    "Palermo": (38.116, 13.361),
    "Campi Flegrei: the Solfatara, Pozzuoli": (40.827, 14.139),
    "Campi Flegrei: southern rim near Capo Miseno": (40.785, 14.085),
    "Naples, city centre": (40.852, 14.268),
    "Vesuvius": (40.821, 14.426),
    "Epicentre of 7 September 2009, Mw 4.7 (INGV)": (38.728, 14.044),
    "Epicentre of 6 September 2002, ML 5.6, offshore Palermo (INGV)": (38.381, 13.701),
}

for name, e in EPICENTRES.items():
    print(name)
    for place, p in PLACES.items():
        print("   %-64s %6.1f km" % (place, km(e, p)))

print("\nTime gaps")
print("  USGS M 4.9 to M 5.7: ", datetime(2026, 10, 23, 10, 19, 16) - datetime(2026, 10, 23, 9, 51, 39))
print("  Tibet M 7.0, 14:33:40 UTC Oct 22 = Beijing", datetime(2026, 10, 22, 14, 33, 40) + timedelta(hours=8))
upload = datetime(2026, 10, 23, 11, 8)  # the episode's upload, UTC
print("  Upload 11:08 UTC Oct 23 = Phoenix (MST, UTC-7)", upload - timedelta(hours=7))
print("  Upload minus Tibet quake:", upload - datetime(2026, 10, 22, 14, 33, 40))
print("  Upload minus Sicily M 5.7:", upload - datetime(2026, 10, 23, 10, 19, 16))
hail = datetime(2026, 10, 21, 15, 25)  # first of 18 Phoenix hail reports, MST (NWS Phoenix)
prev = datetime(2026, 10, 22, 11, 55) - timedelta(hours=7)  # the October 22 episode's upload, in MST
print("  October 22 episode in MST:", prev, "; hours after the Phoenix hail began:", (prev - hail))