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;AU=149597870.7;RP=6371;RADARC=180/math.pi*3600

def phase(deg):
 if not math.isfinite(deg) or not 0<=deg<=180:raise ValueError('Phase domain0..180 degrees')
 if deg==180:return 0.
 a=math.radians(deg);return (math.sin(a)+(math.pi-a)*math.cos(a))/math.pi
rows=[dict(phase_deg=k,phase_function=phase(k),flux_ratio=.3*(RP/AU)**2*phase(k))for k in range(181)]
degeneracy=[]
for k in range(50,201):
 r=k/100;ag=.3/r**2
 if ag<=.6:degeneracy.append(dict(radius_reference_units=r,geometric_albedo=ag,full_phase_flux_ratio=ag*(r*RP/AU)**2,passive_lambert_energy_consistent=ag<=2/3))
sep=[dict(distance_pc=d,projected_separation_AU=1,separation_arcsec=1/d,aperture_m=D,wavelength_nm=550,rayleigh_arcsec=1.22*550e-9/D*RADARC)for d in range(1,51)for D in range(1,11)]
# Independent visible-sphere integration: z=mu, area element dmu dphi.
def sphere(deg,n):
 mu=(np.arange(n)+.5)/n;phi=(np.arange(2*n)+.5)*math.pi/n;a=math.radians(deg)
 inc=np.sqrt(1-mu[:,None]**2)*np.cos(phi)[None,:]*math.sin(a)+mu[:,None]*math.cos(a)
 integral=np.sum(mu[:,None]*np.maximum(inc,0))/n*(math.pi/n)
 return integral/(2*math.pi/3)
checks=[]
for a in [0,30,60,90,120,150,180]:
 low=sphere(a,400);hi=sphere(a,800);checks.append(dict(phase_deg=a,analytic=phase(a),quadrature_400=low,quadrature_800=hi,absolute_error=abs(hi-phase(a))))
assert max(r['absolute_error']for r in checks)<2e-5
assert rows[0]['phase_function']==1 and rows[-1]['phase_function']==0 and abs(phase(90)-1/math.pi)<1e-15
assert all(a['phase_function']>=b['phase_function']>=0 for a,b in zip(rows,rows[1:]))
assert all(abs(r['full_phase_flux_ratio']/(.3*(RP/AU)**2)-1)<1e-14 for r in degeneracy)
assert sep[0]['separation_arcsec']==1
for k in [-1,181,float('nan')]:
 try:phase(k)
 except ValueError:pass
 else:raise AssertionError(k)
def save(n,rs):
 with(P/n).open('w',newline='',encoding='utf8')as f:w=csv.DictWriter(f,fieldnames=list(rs[0]));w.writeheader();w.writerows(rs)
for n,rs in [('phase.csv',rows),('brightness-degeneracy.csv',degeneracy),('separation.csv',sep),('sphere-check.csv',checks)]:save(n,rs)
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.5),layout='constrained');ax[0].plot([r['phase_deg']for r in rows],[r['flux_ratio']for r in rows]);ax[0].set(xlabel='Star–planet–observer phase angle (degrees)',ylabel='Planet/star reflected flux ratio',title='A  Assigned geometric albedo 0.3; radius 6,371 km')
ax[1].plot([r['radius_reference_units']for r in degeneracy],[r['geometric_albedo']for r in degeneracy]);ax[1].axhline(2/3,ls='--',color='gray',label='Passive Lambert sphere energy bound');ax[1].set(xlabel='Radius / 6,371 km',ylabel='Geometric albedo giving the same full-phase flux',title='B  Brightness alone leaves an ambiguity');ax[1].legend(fontsize=8);fig.suptitle('Ideal reflected-light geometry at 1 AU; no detection or habitability claim')
for e in ['png','svg','pdf']:fig.savefig(P/f'phase-and-brightness.{e}')
plt.close(fig)
fig,ax=plt.subplots(figsize=(8,5),layout='constrained');ds=np.arange(1,51);ax.plot(ds,1/ds,label='Projected 1 AU separation',color='#7957a8')
for D,dash in zip([1,3,10],['--','-.',':']):
 ax.axhline(1.22*550e-9/D*RADARC,ls=dash,label=f'Ideal Rayleigh scale: {D} m aperture')
 ax.text(49,1.22*550e-9/D*RADARC*1.06,f'{D} m',ha='right',va='bottom',fontsize=9)
ax.set(xlabel='Distance (pc)',ylabel='Angular scale (arcseconds)',yscale='log',title='Angular separation and diffraction are separate constraints\n550 nm; projected separation is an assigned sky-plane distance');ax.legend(fontsize=9);ax.grid(alpha=.2)
for e in ['png','svg','pdf']:fig.savefig(P/f'separation-and-diffraction.{e}')
plt.close(fig)
(P/'validation.json').write_text(json.dumps(dict(all_pass=True,phase_rows=181,degeneracy_rows=len(degeneracy),separation_rows=500,quadrature_rows=7,max_sphere_error=max(r['absolute_error']for r in checks),endpoints_and_monotonicity=True,brightness_scaling=True,invalid_phase_guards=3),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(),AU_km=AU,reference_radius_km=RP,default_albedo=.3,default_projected_separation_AU=1,default_distance_pc=10,wavelength_nm=550,scope='Uniform Lambert phase function; formal geometric-albedo degeneracy capped at0.6 with additional passive-Lambert bound flag',limits=['No glint clouds atmosphere thermal emission or habitability score','No contrast floor coronagraph inner working angle photon noise or detection probability','Projected separation is not generally orbital radius']),indent=2)+'\n')
print('PASS',len(degeneracy),max(r['absolute_error']for r in checks))