import urllib.request, csv, io, json, math
from pathlib import Path
from collections import Counter
H=Path(__file__).resolve().parent
q=("select+pl_name,pl_rade,pl_orbper,pl_bmasse,st_teff,st_rad,discoverymethod "
   "from+ps+where+pl_rade+is+not+null+and+pl_orbper+is+not+null")
q=q.replace(" ","+")
u=f"https://exoplanetarchive.ipac.caltech.edu/TAP/sync?query={q}&format=csv"
r=urllib.request.Request(u,headers={"User-Agent":"incoprea-research/1.0"})
body=urllib.request.urlopen(r,timeout=120).read().decode()
rows=list(csv.DictReader(io.StringIO(body)))
print("planets:",len(rows))
data=[]
for x in rows:
    try:
        rad=float(x["pl_rade"]); per=float(x["pl_orbper"])
    except: continue
    if 0.3<=rad<=6 and 0.2<=per<=1000: data.append((rad,per))
print("usable:",len(data))
def valley(sub):
    if len(sub)<80: return None
    bins=[(0.5+0.1*i,0.6+0.1*i) for i in range(45)]
    h=[sum(1 for r,_ in sub if lo<=r<hi) for lo,hi in bins]
    i=min(range(1,len(h)-1), key=lambda k:h[k] if h[k-1]>0 and h[k+1]>0 else 10**9)
    return {"n":len(sub),"valley_Rearth":round((bins[i][0]+bins[i][1])/2,2),"count_at_min":h[i]}
allv=valley(data)
short=[d for d in data if d[1]<10]
long_=[d for d in data if d[1]>=10]
vs=valley(short); vl=valley(long_)
# radius distribution table
bins=[(0.5+0.1*i,0.6+0.1*i) for i in range(35)]
dist=[{"radius_bin":f"{(lo+hi)/2:.1f}","count":sum(1 for r,_ in data if lo<=r<hi)} for lo,hi in bins]
(H/"raw_counts.json").write_text(json.dumps({"total":len(rows),"usable":len(data),
 "all":allv,"period_lt10":vs,"period_ge10":vl},indent=2))
with (H/"radius_distribution.csv").open("w",newline="") as f:
    w=csv.DictWriter(f,fieldnames=["radius_bin","count"]); w.writeheader(); w.writerows(dist)
L=["# The exoplanet radius valley","",
f"NASA Exoplanet Archive, confirmed planets with radius and period: **{len(data):,}** usable of {len(rows):,}.","",
"Radius histogram in 0.1 R⊕ bins from 0.5 to 6 R⊕; the valley is the local minimum between the",
"super-Earth peak and the sub-Neptune peak.","",
"| Sample | N | Valley position (R⊕) |","|---|---:|---:|"]
for lab,v in [("All",allv),("P < 10 d",vs),("P ≥ 10 d",vl)]:
    L.append(f"| {lab} | {v['n']:,} | {v['valley_Rearth']:.2f} |" if v else f"| {lab} | — | insufficient |")
L+=["","## What this is","",
"The radius valley is a real, measured feature of the exoplanet population — a deficit of planets near",
"1.8–2 R⊕ between the rocky super-Earths and the volatile-rich sub-Neptunes. It is thought to arise from",
"photoevaporation and core-powered mass loss stripping envelopes from small planets.","",
"If the valley position shifts with orbital period, that is a signature of the mechanism: closer planets",
"receive more irradiation and lose envelopes more efficiently.","",
"## How to read the result","",
"Compare the two period bins above. A measurable shift in valley position between them is consistent with",
"irradiation-driven mass loss; no shift would be evidence against the simple picture. This run measures the",
"shift; it does not model the physics.","",
"## Limits","",
"- Detection bias: small and long-period planets are under-detected, which distorts the histogram.",
"- Radius precision varies by method; transit radii are precise, others are not.",
"- The archive is a heterogeneous compilation, not a homogeneous survey.",
"- No stellar-radius or insolation correction is applied here; period is a crude proxy.",
"- Bin width and the minimum-finding rule affect the measured valley position.",""]
(H/"FINDINGS.md").write_text("\n".join(L))
print(json.dumps({"all":allv,"lt10":vs,"ge10":vl},indent=1))