import math,json,csv,sys,hashlib
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
P=Path(__file__).resolve().parent;G=6.67430e-11

def bulk(mass,radius_km):
    if not all(math.isfinite(v)for v in [mass,radius_km])or not .1e24<=mass<=10e24 or not 2000<=radius_km<=12000:raise ValueError('Outside teaching range')
    r=radius_km*1000;g=G*mass/r**2;rho=3*mass/(4*math.pi*r**3)
    return dict(mass_kg=mass,radius_km=radius_km,gravity_m_s2=g,mean_density_kg_m3=rho,gravity_std=g*math.hypot(.05,2*.02),density_std=rho*math.hypot(.05,3*.02),gravity_min=G*(mass*.95)/(r*1.02)**2,gravity_max=G*(mass*1.05)/(r*.98)**2,density_min=3*mass*.95/(4*math.pi*(r*1.02)**3),density_max=3*mass*1.05/(4*math.pi*(r*.98)**3))

def interior(x,mantle=3500):
    if not math.isfinite(x)or not .2<=x<=.8 or not math.isfinite(mantle)or not 1000<=mantle<=5500:raise ValueError('Outside teaching range')
    core=(5500-mantle*(1-x**3))/x**3;inertia=.4*(core*x**5+mantle*(1-x**5))/5500
    return dict(core_radius_fraction=x,mantle_density_kg_m3=mantle,core_density_kg_m3=core,mean_density_kg_m3=5500,normalized_inertia=inertia)
grid=[bulk((.1+9.9*i/20)*1e24,2000+10000*j/20)for i in range(21)for j in range(21)]
sweep=[bulk(6e24,2000+100*k)for k in range(101)]
ins=[interior(k/100)for k in range(20,81)]
def save(name,rows):
    with(P/name).open('w',newline='',encoding='utf8')as f:w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
for n,rr in [('mass-radius-grid.csv',grid),('radius-uncertainty.csv',sweep),('interior.csv',ins)]:save(n,rr)
maxmass=maxinertia=0
for row in grid:
    vol=4/3*math.pi*(row['radius_km']*1000)**3;maxmass=max(maxmass,abs(row['mean_density_kg_m3']*vol/row['mass_kg']-1))
for x in [k/100 for k in range(20,81)]:
    for mantle in range(1000,5501,100):
        a=interior(x,mantle);c=a['core_density_kg_m3'];mass=4*math.pi/3*(c*x**3+mantle*(1-x**3));moment=8*math.pi/15*(c*x**5+mantle*(1-x**5));maxinertia=max(maxinertia,abs(moment/mass-a['normalized_inertia']))
        assert abs(mass/(4*math.pi/3)-5500)<1e-9
    assert abs(interior(x,5500)['normalized_inertia']-.4)<1e-14
assert maxmass<1e-14 and maxinertia<1e-14
# Central finite differences against analytic partial derivatives at neutral default.
m=6e24;r=6e6;dm=m*1e-5;dr=r*1e-5
functions=[(lambda mm,rr:G*mm/rr**2,G/r**2,-2*G*m/r**3),(lambda mm,rr:3*mm/(4*math.pi*rr**3),3/(4*math.pi*r**3),-9*m/(4*math.pi*r**4))]
maxderiv=0
for f,am,ar in functions:
    maxderiv=max(maxderiv,abs(((f(m+dm,r)-f(m-dm,r))/(2*dm))/am-1),abs(((f(m,r+dr)-f(m,r-dr))/(2*dr))/ar-1))
assert maxderiv<1e-8
bad=0
for fn,args in [(bulk,(0,6000)),(bulk,(6e24,0)),(bulk,(math.nan,6000)),(interior,(0,3500)),(interior,(.5,0)),(interior,(.9,3500))]:
    try:fn(*args)
    except ValueError:bad+=1
    else:raise AssertionError(args)
plt.rcParams.update({'font.size':10,'axes.spines.top':False,'axes.spines.right':False,'savefig.dpi':180})
fig,ax=plt.subplots(1,2,figsize=(11,4.8),layout='constrained')
for a,field,unit in [(ax[0],'gravity','Surface gravity (m/s²)'),(ax[1],'density','Mean density (kg/m³)')]:
    actual='gravity_m_s2'if field=='gravity'else'mean_density_kg_m3';xs=[r['radius_km']for r in sweep];ys=[r[actual]for r in sweep]
    a.plot(xs,ys,label='Nominal relation',color='#7957a8');a.fill_between(xs,[r[field+'_min']for r in sweep],[r[field+'_max']for r in sweep],alpha=.2,label='Assigned bounded-input extrema');a.set(xlabel='Radius (km)',ylabel=unit);a.legend(fontsize=8);a.grid(alpha=.2)
fig.suptitle('Hypothetical fixed mass 6 × 10²⁴ kg; mass ±5%, radius ±2%\nShading is a bounded-input envelope, not a confidence interval')
for ext in ['png','svg','pdf']:fig.savefig(P/f'bulk-and-bounds.{ext}')
plt.close(fig)
fig,ax=plt.subplots(1,2,figsize=(11,4.8),layout='constrained');xs=[r['core_radius_fraction']for r in ins]
ax[0].semilogy(xs,[r['core_density_kg_m3']for r in ins],color='#7957a8');ax[0].set(xlabel='Core radius / outer radius',ylabel='Required core density (kg/m³)',title='A  Small cores require very large densities')
ax[1].plot(xs,[r['normalized_inertia']for r in ins],color='#278c89');ax[1].axhline(.4,ls='--',color='gray',label='Uniform-density sphere');ax[1].set(xlabel='Core radius / outer radius',ylabel='Moment of inertia / (mass × radius²)',title='B  Equal mean density, different mass distribution');ax[1].legend(fontsize=8)
fig.suptitle('Two constant-density layers: mean 5,500; mantle 3,500 kg/m³\nMathematical alternatives, not physically validated interiors')
for ext in ['png','svg','pdf']:fig.savefig(P/f'interior-degeneracy.{ext}')
plt.close(fig)
(P/'validation.json').write_text(json.dumps(dict(all_pass=True,rows=dict(grid=441,radius=101,interior=61),interior_checks=2806,max_mass_relative_error=maxmass,max_inertia_error=maxinertia,max_derivative_relative_error=maxderiv,invalid_guards=bad),indent=2)+'\n')
(P/'study.json').write_text(json.dumps(dict(runtime=dict(python=sys.version,numpy=np.__version__,matplotlib=matplotlib.__version__),generator_sha256=hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),G_SI=G,default_mass_kg=6e24,default_radius_km=6000,assigned_independent_standard_uncertainties=dict(mass_relative=.05,radius_relative=.02),separate_bounded_input_errors=dict(mass_fraction=.05,radius_fraction=.02),interior_mean_density=5500,interior_mantle_default=3500,limits=['No equation of state or physical plausibility filter','No rotation flattening atmosphere or altitude correction','G treated exact for this exercise','Standard uncertainty and bounded extrema are distinct assumptions']),indent=2)+'\n')
print('PASS603rows',maxderiv)