import json, urllib.request, math, csv, time
from pathlib import Path
from collections import Counter
H=Path(__file__).resolve().parent
UA={"User-Agent":"incoprea-research/1.0"}
REGIONS={
 "Japan":(30,46,128,146,"subduction"),
 "Chile":(-56,-17,-76,-66,"subduction"),
 "California":(32,42,-125,-114,"transform"),
 "Iceland":(63,67,-25,-13,"rift"),
}
def fetch(box):
    la0,la1,lo0,lo1=box
    u=("https://earthquake.usgs.gov/fdsnws/event/1/query?format=geojson"
       f"&starttime=2015-01-01&endtime=2026-01-01&minmagnitude=4.0"
       f"&minlatitude={la0}&maxlatitude={la1}&minlongitude={lo0}&maxlongitude={lo1}&limit=20000")
    d=json.loads(urllib.request.urlopen(urllib.request.Request(u,headers=UA),timeout=120).read())
    return [f["properties"]["mag"] for f in d["features"] if f["properties"].get("mag") is not None]
def bvalue(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)<50: return None
    mean=sum(m)/len(m); b=math.log10(math.e)/(mean-(Mc-0.05))
    return {"Mc":Mc,"n":len(m),"b":round(b,3),"se":round(b/math.sqrt(len(m)),3)}
rows=[]
for name,(la0,la1,lo0,lo1,setting) in REGIONS.items():
    try:
        mags=fetch((la0,la1,lo0,lo1)); time.sleep(0.3)
        r=bvalue(mags)
        if r:
            r.update({"region":name,"setting":setting,"total_M4plus":len(mags)})
            rows.append(r); print(name,r)
        else: print(name,"too few events")
    except Exception as e:
        print(name,"FAILED",repr(e))
if rows:
    cols=["region","setting","Mc","n","b","se","total_M4plus"]
    with (H/"bvalues.csv").open("w",newline="") as f:
        w=csv.DictWriter(f,fieldnames=cols); w.writeheader()
        for r in rows: w.writerow({k:r[k] for k in cols})
    L=["# Gutenberg-Richter b-values across tectonic settings","",
    "USGS catalogue, M ≥ 4.0, 2015–2025. Magnitude of completeness from maximum curvature;",
    "b from the Aki maximum-likelihood estimator, SE = b/√N.","",
    "| Region | Setting | Mc | N ≥ Mc | b | SE |","|---|---|---:|---:|---:|---:|"]
    for r in rows:
        L.append(f"| {r['region']} | {r['setting']} | {r['Mc']} | {r['n']:,} | {r['b']:.2f} | {r['se']:.2f} |")
    L+=["","## What this is","",
    "b is the slope of the frequency–magnitude distribution — the ratio of small to large earthquakes.",
    "It varies with stress state, fault heterogeneity and depth. These values are computed here from",
    "the catalogue, not recalled. Whether the ordering matches published regional b-values is a check",
    "this note does not perform, so treat it as a measurement on real data, not a new result.","",
    "## Limits","",
    "- Catalogue completeness differs by network density; Mc differences propagate into b.",
    "- One fixed magnitude threshold for all regions.",
    "- No depth separation, though subduction b varies strongly with depth.",
    "- Aftershocks violate independence and inflate N.",""]
    (H/"FINDINGS.md").write_text("\n".join(L))
    print("wrote",len(rows),"regions")
else:
    (H/"FINDINGS.md").write_text("# Seismicity\n\nNo region returned a usable sample. No results reported.\n")