"""Distances and time gaps on the page "S0 News, October 24, 2026" (factchecked.science).

Great-circle distances by the haversine formula on a sphere of radius 6,371 km, rounded to the kilometre.
Track points: the Japan Meteorological Agency's analyses of Typhoon No. 30 (Champi), as listed by Digital Typhoon
(National Institute of Informatics), and the Joint Typhoon Warning Center's final warning. Epicentres: the USGS's
for event us6000u4j2 and the warning centres' preliminary one. Place coordinates are approximate, to about 0.01°.
Times are UTC. Needs only Python 3: python3 -I distances.py
"""
from math import asin, cos, radians, sin, sqrt


def km(a, b):
    (la1, lo1), (la2, lo2) = a, b
    p1, p2 = radians(la1), radians(la2)
    h = sin((p2 - p1) / 2) ** 2 + cos(p1) * cos(p2) * sin(radians(lo2 - lo1) / 2) ** 2
    return 2 * 6371.0 * asin(sqrt(h))


PLACES = {
    "Sandiao Cape, Taiwan's easternmost point": (25.008, 122.003),
    "Su-ao, Yilan": (24.595, 121.865),
    "Keelung": (25.130, 121.740),
    "Yonaguni, Japan": (24.467, 123.004),
    "Ishigaki, Japan": (24.340, 124.156),
    "Kagoshima, Kyushu": (31.597, 130.557),
    "Nagasaki, Kyushu": (32.750, 129.877),
}
# JMA analyses (date and hour UTC, latitude N, longitude E)
JMA = [("23 Oct 06", 23.8, 123.3), ("23 Oct 09", 24.2, 123.0), ("23 Oct 12", 24.9, 122.8),
       ("23 Oct 15", 25.0, 122.4), ("23 Oct 18", 24.8, 122.3), ("23 Oct 21", 24.9, 122.7),
       ("24 Oct 00", 24.9, 123.3), ("24 Oct 06", 26.7, 123.3), ("24 Oct 12", 27.7, 123.1),
       ("25 Oct 00", 29.0, 124.5), ("25 Oct 06", 30.0, 125.4), ("25 Oct 12", 31.0, 126.0)]
# JTWC: the final warning's position and its 12-hour forecast point
JTWC = [("25 Oct 18 (JTWC)", 31.2, 127.4), ("26 Oct 06 (JTWC forecast)", 31.9, 130.5)]

print("Champi: the closest track point to each place")
for name, pos in PLACES.items():
    d, t = min((km((la, lo), pos), t) for t, la, lo in JMA + JTWC)
    print(f"  {name:42s} {d:5.0f} km  ({t})")
print("\nChampi on 25 October: distance to Kagoshima and Nagasaki")
for t, la, lo in JMA[-3:] + JTWC:
    print(f"  {t:26s} Kagoshima {km((la, lo), PLACES['Kagoshima, Kyushu']):4.0f} km   Nagasaki {km((la, lo), PLACES['Nagasaki, Kyushu']):4.0f} km")

EPI = (-14.0558, 165.5455)       # USGS, us6000u4j2
PRELIM = (-13.8, 165.7)          # PTWC and NTWC, first messages
VANUATU = {"Port-Olry, Espiritu Santo": (-15.0494, 167.0511), "Sola, Vanua Lava": (-13.876, 167.551),
           "Luganville": (-15.5333, 167.1667), "Port Vila": (-17.734, 168.322)}
print("\nVanuatu M 6.9 of 24 October: the USGS epicentre to")
for name, pos in VANUATU.items():
    print(f"  {name:28s} {km(EPI, pos):4.0f} km")
print(f"  the warning centres' preliminary epicentre: {km(EPI, PRELIM):.0f} km")


def minutes(h, m, s=0):
    return h * 60 + m + s / 60


quake = minutes(2, 31, 28.6)
print("\nTime gaps, 24 October (UTC)")
for label, t in [("PTWC threat message 1, M 7.2 (02:37)", minutes(2, 37)),
                 ("NTWC statement, M 7.2 (02:39)", minutes(2, 39)),
                 ("USGS origin, Mww 6.9 (02:45:49)", minutes(2, 45, 49)),
                 ("PTWC final message, M 6.9 (03:28)", minutes(3, 28))]:
    print(f"  {label:40s} {t - quake:4.0f} min after the quake")
print(f"  PTWC threat message 1 to its final message: {minutes(3, 28) - minutes(2, 37):.0f} min")
episode = minutes(11, 34, 42)
print(f"  the episode (11:34:42): {(episode - quake) / 60:.1f} h after the quake")
print(f"  Taiwan's land warning lifted (06:30): {(episode - minutes(6, 30)) / 60:.1f} h before the episode")
print(f"  IMD's press release (14:30 IST = 09:00): {(episode - minutes(9, 0)) / 60:.1f} h before the episode")