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;NA=6.02214076e23;SIGMA=1.10e-38;RHO=2.7

def radar(t,v):
 if not math.isfinite(t)or not 0<=t<=200 or not math.isfinite(v)or not .05<=v<=.15:raise ValueError('Invalid radar inputs')
 return t*v/2

def column(length):
 if not math.isfinite(length)or not .1<=length<=100:raise ValueError('Invalid path')
 col=RHO*length*100*NA;tau=col*SIGMA
 return dict(path_m=length,nucleon_column_cm2=col,optical_depth=tau,interaction_probability=-math.expm1(-tau))

def counts(n,p,eta):
 if not all(math.isfinite(x)for x in [n,p,eta])or n<0 or not 0<=p<=1 or not 0<=eta<=1:raise ValueError('Invalid count inputs')
 mean=n*p*eta
 return mean,None if mean==0 else 1/math.sqrt(mean)
rr=[dict(two_way_time_ns=t,velocity_m_ns=v,depth_m=radar(t,v))for v in [.05,.1,.15]for t in range(201)]
cols=[column(10**(-1+3*k/100))for k in range(101)]
budgets=[]
for r in cols:
 for logn in range(9,17):
  for eta in [.1,.5,1]:
   mean,rel=counts(10**logn,r['interaction_probability'],eta);budgets.append(dict(path_m=r['path_m'],incident_particles=10**logn,efficiency=eta,expected_count=mean,relative_Poisson_std=rel))
maxTaylor=0
for r in cols:
 tau=r['optical_depth'];ref=tau-tau*tau/2+tau**3/6;maxTaylor=max(maxTaylor,abs(r['interaction_probability']/ref-1))
 # SI column check: kg/m3 * m * nucleons/kg * cm2/m2.
 independent=2700*r['path_m']*(NA*1000)/10000
 assert math.isclose(independent,r['nucleon_column_cm2'],rel_tol=1e-14)
assert maxTaylor<1e-14 and radar(100,.1)==5 and counts(0,.1,1)==(0,None) and counts(1,.1,0)==(0,None)
assert counts(4,.1,1)[1]==counts(1,.1,1)[1]/2
bad=0
for f,a in [(radar,(-1,.1)),(radar,(201,.1)),(radar,(100,.01)),(column,(0,)),(column,(101,)),(counts,(-1,.1,1)),(counts,(1,1.1,1)),(counts,(1,.1,1.1))]:
 try:f(*a)
 except ValueError:bad+=1
 else:raise AssertionError(a)
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,r in [('radar.csv',rr),('particle-columns.csv',cols),('count-budgets.csv',budgets)]:save(n,r)
profiles=[dict(lower_concentration_mg_m3=k/10,upper_concentration_mg_m3=2-k/10,slab_thickness_m=1000,column_mg_m2=1000*(k/10+2-k/10)) for k in range(21)]
assert all(abs(r['column_mg_m2']-2000)<1e-10 for r in profiles)
save('profiles.csv',profiles)
plt.rcParams.update({'font.size':10,'axes.spines.top':False,'axes.spines.right':False,'savefig.dpi':180})
fig,ax=plt.subplots(figsize=(8,4.8),layout='constrained')
for v in [.05,.1,.15]:
 a=[r for r in rr if r['velocity_m_ns']==v];ax.plot([r['two_way_time_ns']for r in a],[r['depth_m']for r in a],label=f'Assumed speed {v:.2f} m/ns')
ax.axvline(100,ls='--',color='gray');ax.set(xlabel='Vertical two-way travel time (ns)',ylabel='Inferred depth (m)',title='One radar travel time gives different depths\nUniform assumed medium; no composition or utility-map claim');ax.legend();ax.grid(alpha=.2)
for e in ['png','svg','pdf']:fig.savefig(P/f'radar-velocity.{e}')
plt.close(fig)
fig,ax=plt.subplots(1,2,figsize=(11,4.8),layout='constrained');ax[0].loglog([r['path_m']for r in cols],[r['interaction_probability']for r in cols]);ax[0].set(xlabel='Hypothetical path through material (m)',ylabel='Interaction probability per incident particle',title='A  Assigned density and cross section')
for length in [.1,1,100]:
 ns=np.logspace(9,16,101);p=column(length)['interaction_probability'];ax[1].loglog(ns,[counts(n,p,.5)[1]for n in ns],label=f'{length:g} m path')
ax[1].set(xlabel='Incident particle count',ylabel='Ideal relative Poisson standard deviation',title='B  Assigned detection efficiency 0.5');ax[1].legend();fig.suptitle('Conditional rare-interaction model: no detector or spatial-resolution prediction')
for e in ['png','svg','pdf']:fig.savefig(P/f'rare-interactions.{e}')
plt.close(fig)
(P/'validation.json').write_text(json.dumps(dict(all_pass=True,rows=dict(radar=603,columns=101,budgets=2424,profiles=21),max_taylor_relative_error=maxTaylor,independent_dimensions=True,zero_count_undefined=True,poisson_scaling=True,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(),density_g_cm3=RHO,nucleons_per_gram=NA,cross_section_cm2_per_nucleon=SIGMA,cross_section_reported_uncertainty_not_propagated=.15e-38,scope='Nominal iron charged-current cross section at about1.1GeV transferred approximately to generic rock; not material-specific prediction',assumptions=['Constant cross section','Independent ideal Poisson counts','No background dead time or correlation','No heating or operational source model']),indent=2)+'\n');print('PASS3149rows',maxTaylor)