"""The sealed rule for "earthquakes come at the full and new moon". Sealed at the 1f916 registry before the catalogue for this study was fetched and before a single lunar phase was computed. What is already on this server is `~/data/comcat_m65.csv`, an M>=6.5 event list downloaded on 2026-09-11 for an unrelated study about how often big earthquakes happen; its times have never been turned into a moon phase by me or by anything I have run. This study fetches its own M>=6.0 list. THE CLAIM --------- People say big earthquakes come at the full moon, or at the full and new moon together -- the "supermoon earthquake". The serious version of the same idea is tidal triggering: the sun and moon in line raise the largest solid-earth and ocean tides, the extra stress is small but real, and it might push a fault already close to failure over the edge. The popular version says you could notice it. This study tests the popular version on the global catalogue. POPULATION ---------- USGS ComCat, FDSN event service, `eventtype=earthquake`, `minmagnitude=6.0`, origin time from 1973-01-01T00:00:00Z to 2025-12-31T23:59:59Z, worldwide, all depths. Fetched year by year as GeoJSON/CSV and concatenated. No other filter. 1973 is where the previous study found ComCat's global M>=6 reporting to be stable; 2025 is the last complete year. THE PHASE OF THE MOON --------------------- For each event, theta = the geocentric elongation of the moon from the sun at the origin instant, in degrees, 0 at new moon and 180 at full, computed with `ephem` (`ephem.Moon(...).elong`, wrapped into [0, 360)). The tidal claim is twice a month, because the tide-raising force does not care which of the two bodies is in front. So the quantity is the SYZYGY DISTANCE d = min(theta, 360 - theta, |theta - 180|) in [0, 90] which is 0 at both the new and the full moon and 90 at both quarters. "Near syzygy" means d <= 45: exactly half the circle of phase. THE STATISTIC ------------- f_obs = the fraction of events with d <= 45. f_exp is NOT 0.5. The moon's elongation does not advance at a constant rate -- the orbit is elliptical -- so the catalogue window does not contain equal amounts of syzygy time and quadrature time. f_exp is the fraction of the 1973-2025 window spent at d <= 45, computed by evaluating d on a grid of one point per hour across the whole window and taking the mean. Using 0.5 here would put the orbit's own asymmetry into the answer. E = f_obs / f_exp - 1 the SYZYGY EXCESS, reported as a per cent. E = 0 means earthquakes fall on lunar phase exactly as time does. E = +0.10 means ten per cent more of them happen near syzygy than time alone predicts. AFTERSHOCKS ARE THE WHOLE DIFFICULTY ------------------------------------ Earthquakes are not independent. One M7 can be followed by hundreds of events inside a week and a hundred kilometres, all at nearly the same lunar phase as their mainshock. Counting them as separate evidence would shrink every interval by a factor that has nothing to do with the moon. Two answers, both pre-registered: 1. DECLUSTER, and make the mainshock the unit. Gardner & Knopoff (1974) windows, in the standard form: L(M) = 10 ** (0.1238 * M + 0.983) km, epicentral T(M) = 10 ** (0.032 * M + 2.7389) days, for M >= 6.5 T(M) = 10 ** (0.5409 * M - 0.547) days, for M < 6.5 Algorithm, fixed here so it is reproducible: take the events in order of decreasing magnitude (ties broken by earlier origin time). Each event not already assigned becomes a mainshock and claims, into its own cluster, every unassigned event whose origin time lies in [t, t + T(M)] and whose epicentre lies within L(M) great-circle kilometres. The set of mainshocks is the declustered catalogue. Windows run FORWARD only, as Gardner and Knopoff wrote them; a symmetric variant (|dt| <= T) is a pre-registered robustness check, not the primary. 2. BOOTSTRAP THE CLUSTER, NOT THE EVENT. The 95 % interval on E comes from 10,000 resamples. The primary estimate resamples mainshocks with replacement (declustering has already made the unit a cluster). The robustness estimate on the FULL undeclustered catalogue resamples whole Gardner-Knopoff clusters with replacement, keeping every event inside a drawn cluster. An interval whose unit is the event is not computed at all. A circular time shift is NOT used as a null here, and the reason is worth writing down before anyone suggests it: shifting every event by the same lag rotates every lunar phase by nearly the same angle, so any statistic that measures how concentrated the phases are is almost invariant under it. It would produce a null distribution with no width and a p-value of nothing. SCHUSTER -------- Reported beside E, on the declustered catalogue, because it is the classical test for this exact question and a reader will look for it. On the doubled angle (2 * theta, so that new and full moon fall together): R = |sum exp(i * 2 * theta_k)| , p = exp(-R^2 / n) Its p-value assumes independent events, which is why declustering comes first and why the verdict is carried by E and its cluster bootstrap, not by p. THE BAR, AND THE KILL RULE -------------------------- The claim under test is that you could notice it. The bar is +5 %: one event in twenty more than time alone would put near syzygy. KILLED if the 95 % interval on E lies wholly below +5 %. SURVIVED if the 95 % interval on E lies wholly above +5 %. INCONCLUSIVE otherwise. If the interval lies wholly below +5 % and also wholly above 0, the verdict is still KILLED and the entry must say in the same breath that a small real effect was measured, because that is the literature's claim and not the folklore's. THE FLOOR, COMPUTED BEFORE THE VERDICT -------------------------------------- Positive control. Simulate catalogues of the same size as the declustered one, with phases drawn from the true time measure, then inject a syzygy excess of 2, 5, 10 and 20 per cent by resampling. For each injected value, run the whole primary estimator on 200 simulated catalogues. The FLOOR is the smallest injected excess recovered with a 95 % interval excluding 0 in at least 90 % of them. If the floor is above +5 %, then the SURVIVED branch of the kill rule cannot fire on any data, and the study is recorded `void: unfailable` in that direction, said plainly and before the real number is looked at. NULL CONTROLS, both run before the real statistic ------------------------------------------------- N1 Times replaced by uniform draws from the catalogue window, N kept. E must come back near 0 with its interval covering 0, and the false-positive rate of "interval excludes 0" over 200 draws must be near 5 %. N2 Each event's time jittered by a uniform draw from +/- 365 days. This keeps the catalogue's secular and seasonal shape and destroys only the phase. E must come back near 0. A jitter control that failed while N1 passed would mean the statistic was reading something about the calendar. TIME-MEASURE CONTROL -------------------- f_exp computed on the hourly grid must agree to within 0.001 with f_exp computed on a grid of one point per 7 minutes, and both must differ from 0.500 by less than 0.02. A large departure would mean the elongation code is wrong, not that the orbit is strange. SECONDARIES, pre-registered, scored, never the verdict ------------------------------------------------------ S1 Shallow events only, depth <= 70 km. This is where the tidal stress change is largest and where the literature looks. S2 M >= 7.0, all depths -- the "great earthquake" the folklore actually means. Its own floor is computed the same way; I expect it to be above the bar, in which case S2 is void: unfailable and is published as such. S3 Full moon alone: |theta - 180| <= 45 against its own time measure. The folklore usually says full moon, not syzygy. S4 The symmetric declustering window. BRANCH ODDS, WRITTEN BEFORE THE DATA ------------------------------------ P(KILLED) 0.80, P(INCONCLUSIVE) 0.15, P(SURVIVED) 0.05. The literature that I have read about reports tidal triggering effects of a few per cent, in particular settings, needing large catalogues to see at all; a global M>=6 catalogue at a +5 % bar should therefore kill the folkloric version while possibly leaving a small positive estimate standing. WHAT WOULD MAKE THIS STUDY VOID ------------------------------- - The floor exceeding the bar (the SURVIVED branch cannot fire): void in that direction, stated. - Either null control returning a non-zero E: the statistic is wrong and no verdict may be published from it. - Declustering leaving fewer than 1,000 mainshocks: the interval would be too wide for the bar to mean anything, and the floor computation will say so. Everything below is the estimator, written before any of this study's data existed. `analyse.py` imports it and adds nothing to it. """ import math import ephem import numpy as np SEED = 20260922 BAR = 0.05 # the +5 % bar HALF = 45.0 # "near syzygy" is d <= 45 degrees BOOTSTRAPS = 10_000 START = "1973-01-01T00:00:00Z" END = "2025-12-31T23:59:59Z" # ---------------------------------------------------------------- lunar phase def elongation(dt64): """Geocentric elongation of the moon from the sun, degrees in [0, 360).""" d = ephem.Date(_to_pyephem(dt64)) m = ephem.Moon(d) return math.degrees(m.elong) % 360.0 def _to_pyephem(dt64): """numpy datetime64 -> the (y, m, d, H, M, S) tuple pyephem accepts.""" s = np.datetime64(dt64, "s").astype("datetime64[s]").astype(object) return (s.year, s.month, s.day, s.hour, s.minute, s.second) def phases(times): return np.array([elongation(t) for t in times], dtype=float) def syzygy_distance(theta): """0 at new and full moon, 90 at the quarters.""" t = np.asarray(theta, dtype=float) % 360.0 return np.minimum(np.minimum(t, 360.0 - t), np.abs(t - 180.0)) def full_distance(theta): """Angular distance from the full moon alone, 0..180.""" t = np.asarray(theta, dtype=float) % 360.0 return np.abs(((t - 180.0 + 180.0) % 360.0) - 180.0) def time_measure(start, end, step_hours=1.0, distance=syzygy_distance, half=HALF): """Fraction of the window spent inside the phase window. The null expectation. Not 0.5: the elongation rate varies over the month. """ t0 = np.datetime64(start[:-1], "s") t1 = np.datetime64(end[:-1], "s") n = int((t1 - t0) / np.timedelta64(1, "s") // (step_hours * 3600)) + 1 grid = t0 + (np.arange(n) * int(step_hours * 3600)).astype("timedelta64[s]") return float((distance(phases(grid)) <= half).mean()) # ------------------------------------------------------------------- excess def excess(theta, f_exp, distance=syzygy_distance, half=HALF): """The syzygy excess E = f_obs / f_exp - 1.""" f_obs = float((distance(theta) <= half).mean()) return f_obs / f_exp - 1.0 def bootstrap_excess(theta, f_exp, rng, draws=BOOTSTRAPS, distance=syzygy_distance, half=HALF): """Percentile interval, resampling the rows given (mainshocks).""" inside = (distance(theta) <= half).astype(float) n = inside.size idx = rng.integers(0, n, size=(draws, n)) f = inside[idx].mean(axis=1) e = f / f_exp - 1.0 return float(np.percentile(e, 2.5)), float(np.percentile(e, 97.5)) def bootstrap_excess_by_cluster(theta, cluster_id, f_exp, rng, draws=BOOTSTRAPS, distance=syzygy_distance, half=HALF): """Robustness on the full catalogue: whole clusters are the unit.""" inside = (distance(theta) <= half).astype(float) ids = np.unique(cluster_id) members = [np.flatnonzero(cluster_id == c) for c in ids] k = len(members) out = np.empty(draws) for b in range(draws): pick = rng.integers(0, k, size=k) rows = np.concatenate([members[i] for i in pick]) out[b] = inside[rows].mean() / f_exp - 1.0 return float(np.percentile(out, 2.5)), float(np.percentile(out, 97.5)) def schuster(theta, doubled=True): """Classical Schuster test. Returns (R, p). Assumes independence.""" a = np.radians(np.asarray(theta, dtype=float) * (2.0 if doubled else 1.0)) n = a.size r = math.hypot(float(np.cos(a).sum()), float(np.sin(a).sum())) return r, math.exp(-(r ** 2) / n) # -------------------------------------------------------------- declustering def gk_windows(mag): """Gardner & Knopoff (1974): (km, days).""" L = 10.0 ** (0.1238 * mag + 0.983) T = 10.0 ** (0.032 * mag + 2.7389) if mag >= 6.5 else 10.0 ** (0.5409 * mag - 0.547) return L, T def haversine(lat1, lon1, lat2, lon2): p1, p2 = np.radians(lat1), np.radians(lat2) dp = p2 - p1 dl = np.radians(lon2 - lon1) a = np.sin(dp / 2) ** 2 + np.cos(p1) * np.cos(p2) * np.sin(dl / 2) ** 2 return 6371.0 * 2 * np.arcsin(np.sqrt(a)) def decluster(times, lat, lon, mag, symmetric=False): """Gardner-Knopoff. Returns (mainshock_index, cluster_id_per_event).""" n = len(mag) t = np.asarray(times, dtype="datetime64[s]").astype("int64") / 86400.0 order = np.lexsort((t, -np.asarray(mag, dtype=float))) cluster = np.full(n, -1, dtype=int) mains = [] for i in order: if cluster[i] != -1: continue c = len(mains) mains.append(i) cluster[i] = c L, T = gk_windows(float(mag[i])) dt = t - t[i] near_t = (np.abs(dt) <= T) if symmetric else ((dt >= 0) & (dt <= T)) free = (cluster == -1) & near_t if free.any(): j = np.flatnonzero(free) d = haversine(lat[i], lon[i], np.asarray(lat)[j], np.asarray(lon)[j]) cluster[j[d <= L]] = c return np.array(mains, dtype=int), cluster # -------------------------------------------------------- the verdict itself def verdict(lo, hi, bar=BAR): if hi < bar: return "KILLED" if lo > bar: return "SURVIVED" return "INCONCLUSIVE"