import json, urllib.request, math, csv, datetime as dt
from pathlib import Path
from collections import Counter
H=Path(__file__).resolve().parent
UA={"User-Agent":"incoprea-research/1.0"}
u=("https://earthquake.usgs.gov/fdsnws/event/1/query?format=geojson"
   "&starttime=2010-01-01&endtime=2026-01-01&minmagnitude=3.5"
   "&minlatitude=30&maxlatitude=46&minlongitude=128&maxlongitude=146&limit=20000")
d=json.loads(urllib.request.urlopen(urllib.request.Request(u,headers=UA),timeout=120).read())
ev=[]
for f in d["features"]:
    p=f["properties"]; m=p.get("mag"); t=p.get("time")
    if m is None or t is None: continue
    jst=dt.datetime.fromtimestamp(t/1000, dt.timezone.utc)+dt.timedelta(hours=9)
    ev.append((m, jst.hour))
(H/"diurnal_raw.json").write_text(json.dumps({"n":len(ev),"note":"Japan box, M>=3.5, 2010-2025, local hour JST"}))
def bval(mags):
    c=Counter(round(m,1) for m in mags); Mc=max(c.items(), key=lambda kv: kv[1])[0]
    m=[x for x in mags if x>=Mc]
    if len(m)<100: return None
    mean=sum(m)/len(m); b=math.log10(math.e)/(mean-(Mc-0.05))
    return Mc, len(m), b, b/math.sqrt(len(m))
day=[m for m,h in ev if 6<=h<18]
night=[m for m,h in ev if h<6 or h>=18]
res={}
for label,s in [("day (06-18 JST)",day),("night (18-06 JST)",night)]:
    r=bval(s)
    if r: res[label]={"Mc":r[0],"n":r[1],"b":round(r[2],3),"se":round(r[3],3)}
diff=res["night (18-06 JST)"]["b"]-res["day (06-18 JST)"]["b"]
se=math.sqrt(res["night (18-06 JST)"]["se"]**2+res["day (06-18 JST)"]["se"]**2)
(H/"diurnal.json").write_text(json.dumps({"n_events":len(ev),"day":len(day),"night":len(night),
  "results":res,"b_difference":round(diff,3),"se_difference":round(se,3),"z":round(diff/se,2)},indent=2))
with (H/"diurnal.csv").open("w",newline="") as f:
    w=csv.writer(f); w.writerow(["window","Mc","n","b","se"])
    for k,v in res.items(): w.writerow([k,v["Mc"],v["n"],v["b"],v["se"]])
L=["# Does the measured b-value depend on time of day?","",
"Japan box (30–46°N, 128–146°E), M ≥ 3.5, 2010–2025. Events split by local (JST) hour of occurrence.","",
"| Window | Mc | N | b | SE |","|---|---:|---:|---:|---:|"]
for k,v in res.items(): L.append(f"| {k} | {v['Mc']} | {v['n']:,} | {v['b']:.2f} | {v['se']:.2f} |")
L+=["",f"Difference (night − day): **{diff:+.3f} ± {se:.3f}** (z = {diff/se:.1f}).","",
"## Why this matters","",
"Seismic networks detect more small events when cultural noise is low. If that is true here, the",
"magnitude of completeness falls at night, more small events enter the sample, and the fitted b rises —",
"not because the physics changed, but because the detector did.","",
"A diurnal difference is therefore a **catalog artefact**, and a warning: any study comparing b-values",
"across regions or epochs is comparing detection thresholds as well as tectonics.","",
"## What was found","",
"See the table and z-score above. A positive difference with |z| well above 2 supports the artefact",
"hypothesis; a null result means detection bias is not visible at this magnitude threshold in this region.","",
"## Limits","",
"- Single region and threshold; the effect should be strongest near Mc, which is where M ≥ 3.5 sits.",
"- Local hour is a proxy for cultural noise; no station-level noise data used.",
"- Aftershock sequences are not declustered and are unevenly distributed across hours.",
"- USGS 'time' is origin time, not detection time.",""]
(H/"DIURNAL.md").write_text("\n".join(L))
print(json.dumps({"day":res["day (06-18 JST)"],"night":res["night (18-06 JST)"],"diff":round(diff,3),"se":round(se,3),"z":round(diff/se,2)},indent=1))