""" bhukampa_common.py — shared definitions for the "Varāhamihira's Earthquakes" essay (AyurAstro.com). Nakṣatra schemes, the four earthquake circles of Bṛhat Saṁhitā 32, sidereal Moon positions (Swiss Ephemeris, Lahiri ayanāṁśa), Gardner–Knopoff declustering, chi-square tests and power. Requirements: numpy, pandas, scipy, pyswisseph. """ from __future__ import annotations import numpy as np import pandas as pd from scipy import stats try: import swisseph as swe swe.set_sid_mode(swe.SIDM_LAHIRI) _FLAGS = swe.FLG_MOSEPH | swe.FLG_SIDEREAL except ImportError: # pragma: no cover swe = None # ---------------------------------------------------------------- nakṣatras NAK27 = ["Aśvinī", "Bharaṇī", "Kṛttikā", "Rohiṇī", "Mṛgaśīrṣa", "Ārdrā", "Punarvasu", "Puṣya", "Āśleṣā", "Maghā", "Pūrvaphalgunī", "Uttaraphalgunī", "Hasta", "Citrā", "Svātī", "Viśākhā", "Anurādhā", "Jyeṣṭhā", "Mūla", "Pūrvāṣāḍhā", "Uttarāṣāḍhā", "Śravaṇa", "Dhaniṣṭhā", "Śatabhiṣaj", "Pūrvabhādrapadā", "Uttarabhādrapadā", "Revatī"] # 28-scheme: Abhijit occupies the last quarter of Uttarāṣāḍhā (276°40′–280°) # plus the first 1/15 of Śravaṇa (280°–280°53′20″), as in the Sarvatobhadra # and Koṭa Cakra tools on the site. ABHIJIT_START = 276 + 40 / 60 # 276.6667° ABHIJIT_END = 280 + 53 / 60 + 20 / 3600 # 280.8889° NAK28 = NAK27[:21] + ["Abhijit"] + NAK27[21:] # 28-scheme bin edges in degrees (29 edges) EDGES28 = [i * 360 / 27 for i in range(21)] + [ABHIJIT_START, ABHIJIT_END] \ + [i * 360 / 27 for i in range(22, 28)] WIDTH28 = np.diff(EDGES28) # Uttarāṣāḍhā 10°, Abhijit 4°13′20″, Śravaṇa 12°26′40″ WIDTH27 = np.full(27, 360 / 27) # ------------------------------------------------- Bṛhat Saṁhitā 32 circles # 32.8 Vāyu; 32.12 Agni; 32.16 Indra; 32.20 Varuṇa (each seven asterisms; # the Indra circle contains Abhijit, hence 28 in all). MANDALA = { "Vāyu": ["Uttaraphalgunī", "Hasta", "Citrā", "Svātī", "Punarvasu", "Mṛgaśīrṣa", "Aśvinī"], "Agni": ["Puṣya", "Kṛttikā", "Viśākhā", "Bharaṇī", "Maghā", "Pūrvabhādrapadā", "Pūrvaphalgunī"], "Indra": ["Abhijit", "Śravaṇa", "Dhaniṣṭhā", "Rohiṇī", "Jyeṣṭhā", "Uttarāṣāḍhā", "Anurādhā"], "Varuṇa": ["Revatī", "Pūrvāṣāḍhā", "Ārdrā", "Āśleṣā", "Mūla", "Uttarabhādrapadā", "Śatabhiṣaj"], } MANDALA_OF = {n: m for m, ns in MANDALA.items() for n in ns} assert sorted(MANDALA_OF) == sorted(NAK28) def nak27_index(lon): """0..26 for a sidereal longitude in degrees.""" return (np.asarray(lon) % 360 // (360 / 27)).astype(int) def nak28_index(lon): """0..27 (21 = Abhijit) for a sidereal longitude in degrees.""" lon = np.asarray(lon) % 360 return np.clip(np.searchsorted(EDGES28, lon, side="right") - 1, 0, 27) # ---------------------------------------------------------------- the Moon def julian_days(times: pd.Series) -> np.ndarray: """UTC timestamps -> Julian Day (UT).""" t = pd.to_datetime(times, utc=True) days = (t - pd.Timestamp("1970-01-01", tz="UTC")) / pd.Timedelta(days=1) return 2440587.5 + np.asarray(days, float) def moon_sidereal_longitude(jd: np.ndarray) -> np.ndarray: """Geocentric sidereal (Lahiri) longitude of the Moon, degrees. Swiss Ephemeris in Moshier mode — no ephemeris files needed, better than 1″ for the Moon in 1980–2030.""" out = np.empty(len(jd)) for i, j in enumerate(jd): out[i] = swe.calc_ut(float(j), swe.MOON, _FLAGS)[0][0] return out % 360 def planet_sidereal_longitudes(jd: float) -> dict: """All nine grahas (true node) at one JD, sidereal degrees.""" bodies = {"Sun": swe.SUN, "Moon": swe.MOON, "Mars": swe.MARS, "Mercury": swe.MERCURY, "Jupiter": swe.JUPITER, "Venus": swe.VENUS, "Saturn": swe.SATURN, "Rāhu": swe.TRUE_NODE} pos = {k: swe.calc_ut(jd, v, _FLAGS)[0][0] % 360 for k, v in bodies.items()} pos["Ketu"] = (pos["Rāhu"] + 180) % 360 return pos GRAHAS = ["Sun", "Moon", "Mars", "Mercury", "Jupiter", "Venus", "Saturn", "Rāhu", "Ketu"] GRAHA_DN = {"Sun": "सूर्य", "Moon": "चन्द्र", "Mars": "मङ्गल", "Mercury": "बुध", "Jupiter": "गुरु", "Venus": "शुक्र", "Saturn": "शनि", "Rāhu": "राहु", "Ketu": "केतु"} GRAHA_SKT = {"Sun": "Sūrya", "Moon": "Candra", "Mars": "Maṅgala", "Mercury": "Budha", "Jupiter": "Guru", "Venus": "Śukra", "Saturn": "Śani", "Rāhu": "Rāhu", "Ketu": "Ketu"} _BODY = {"Sun": swe.SUN if swe else None, "Moon": swe.MOON if swe else None, "Mars": swe.MARS if swe else None, "Mercury": swe.MERCURY if swe else None, "Jupiter": swe.JUPITER if swe else None, "Venus": swe.VENUS if swe else None, "Saturn": swe.SATURN if swe else None, "Rāhu": swe.TRUE_NODE if swe else None} def graha_sidereal_longitudes(jd: np.ndarray) -> dict: """Sidereal (Lahiri) longitudes of all nine grahas at each JD (true node; Ketu = Rāhu + 180°). Returns {graha: array of degrees}.""" out = {g: np.empty(len(jd)) for g in GRAHAS if g != "Ketu"} for i, j in enumerate(jd): j = float(j) for g, body in _BODY.items(): out[g][i] = swe.calc_ut(j, body, _FLAGS)[0][0] out = {g: v % 360 for g, v in out.items()} out["Ketu"] = (out["Rāhu"] + 180) % 360 return out def ascendant_sidereal(jd: float, lat: float, lon: float) -> float: """Sidereal (Lahiri) ascendant in degrees.""" ayan = swe.get_ayanamsa_ut(jd) cusps, ascmc = swe.houses(jd, lat, lon, b"W") return (ascmc[0] - ayan) % 360 # ------------------------------------------------------------ expectations def time_weighted_expectation(start="1980-01-01", end="2026-01-01", step_hours=6): """Fraction of TIME the Moon spends in each bin over the study window (sampled every `step_hours`). This is the honest null: the Moon moves faster near perigee, so the arcs are not visited in exactly equal time. Returns (p27, p28).""" t = pd.date_range(start, end, freq=f"{step_hours}h", tz="UTC", inclusive="left") lon = moon_sidereal_longitude(julian_days(pd.Series(t))) c27 = np.bincount(nak27_index(lon), minlength=27) c28 = np.bincount(nak28_index(lon), minlength=28) return c27 / c27.sum(), c28 / c28.sum() # ---------------------------------------------------------------- tests def chisq_gof(counts, expected_p): """Chi-square goodness of fit against expected proportions. Returns dict(n, chi2, df, p, w) where w is Cohen's effect size.""" counts = np.asarray(counts, float) n = counts.sum() exp = n * np.asarray(expected_p, float) chi2 = float(((counts - exp) ** 2 / exp).sum()) df = len(counts) - 1 p = float(stats.chi2.sf(chi2, df)) w = float(np.sqrt(chi2 / n)) if n else float("nan") return {"n": int(n), "chi2": chi2, "df": df, "p": p, "w": w} def chisq_power(n, df, w=0.3, alpha=0.05): """Power of the chi-square GOF test to detect effect size w (Cohen's w; 0.1 small, 0.3 medium, 0.5 large) with n observations.""" crit = stats.chi2.ppf(1 - alpha, df) return float(stats.ncx2.sf(crit, df, n * w * w)) def n_for_power(df, w=0.3, alpha=0.05, power=0.8): """Smallest n giving at least `power`.""" lo, hi = 1, 100000 while lo < hi: mid = (lo + hi) // 2 if chisq_power(mid, df, w, alpha) >= power: hi = mid else: lo = mid + 1 return lo def binom_power(n, p0, h=0.5, alpha=0.05): """Exact power of the one-sided binomial test to detect a medium departure (Cohen's h = 0.5) above p0 with n trials.""" p1 = float(np.sin(np.arcsin(np.sqrt(p0)) + h / 2) ** 2) k = int(stats.binom.isf(alpha, n, p0)) + 1 while k > 0 and stats.binom.sf(k - 2, n, p0) <= alpha: k -= 1 while stats.binom.sf(k - 1, n, p0) > alpha: k += 1 return float(stats.binom.sf(k - 1, n, p1)) def n_for_binom_power(p0, h=0.5, alpha=0.05, power=0.8, nmax=2000): """Smallest n from which the exact one-sided binomial test has at least `power` for every larger n (the exact power is not monotone in n).""" ok = [n for n in range(5, nmax) if binom_power(n, p0, h, alpha) >= power] gaps = [n for n in range(5, nmax) if n not in set(ok)] return (max(gaps) + 1) if gaps else 5 def circle_test(nak28_idx, expected_p28, circle="Varuṇa"): """One-sided binomial test: are more quakes in the circle's seven asterisms than the time-weighted null predicts?""" members = [NAK28.index(n) for n in MANDALA[circle]] idx = np.asarray(nak28_idx) k = int(np.isin(idx, members).sum()) n = int(len(idx)) p0 = float(np.asarray(expected_p28)[members].sum()) res = stats.binomtest(k, n, p0, alternative="greater") ci = res.proportion_ci(confidence_level=0.95, method="wilson") return {"circle": circle, "k": k, "n": n, "p0": p0, "observed": k / n if n else float("nan"), "ratio": (k / n) / p0 if n else float("nan"), "p_greater": float(res.pvalue), "p_two_sided": float(stats.binomtest(k, n, p0).pvalue), "ci_low": float(ci.low), "ci_high": float(ci.high)} # ------------------------------------------------------------ declustering def gardner_knopoff_windows(mag): """Gardner & Knopoff (1974) space (km) and time (days) windows.""" mag = np.asarray(mag, float) dist = 10 ** (0.1238 * mag + 0.983) time = np.where(mag >= 6.5, 10 ** (0.032 * mag + 2.7389), 10 ** (0.5409 * mag - 0.547)) return dist, time def decluster_gk(df: pd.DataFrame, time_col="time", lat_col="latitude", lon_col="longitude", mag_col="mag") -> np.ndarray: """Boolean mask: True = mainshock (kept), False = fore/aftershock removed. Classic window method: events are visited from largest to smallest; every smaller event inside the mainshock's space window and within its time window AFTER the mainshock (aftershocks) or within 1/10 of the time window BEFORE it (foreshocks) is removed.""" from scipy.spatial import cKDTree t = np.asarray((pd.to_datetime(df[time_col], utc=True) - pd.Timestamp("1970-01-01", tz="UTC")) / pd.Timedelta(days=1), float) lat = np.radians(df[lat_col].to_numpy(float)) lon = np.radians(df[lon_col].to_numpy(float)) mag = df[mag_col].to_numpy(float) n = len(df) R = 6371.0 xyz = np.c_[np.cos(lat) * np.cos(lon), np.cos(lat) * np.sin(lon), np.sin(lat)] * R tree = cKDTree(xyz) dist_km, time_days = gardner_knopoff_windows(mag) keep = np.ones(n, bool) removed = np.zeros(n, bool) for i in np.argsort(-mag, kind="stable"): if removed[i]: continue chord = 2 * R * np.sin(min(dist_km[i], np.pi * R) / (2 * R)) idx = np.asarray(tree.query_ball_point(xyz[i], chord), int) dt = t[idx] - t[i] sel = (idx != i) & (mag[idx] <= mag[i]) & (dt <= time_days[i]) \ & (dt >= -0.1 * time_days[i]) victims = idx[sel] keep[victims] = False removed[victims] = True return keep