"""Hypothetical photon, gray-surface and solar-receiver accounting."""
import csv,json,math
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
P=Path(__file__).resolve().parent
H=6.62607015e-34;C=299792458.;K=1.380649e-23;E=1.602176634e-19
SIGMA=2*math.pi**5*K**4/(15*H**3*C**2)
def planck(um,t):
 if not math.isfinite(um) or not math.isfinite(t) or um<=0 or t<=0:raise ValueError('Positive finite wavelength and absolute temperature required')
 lam=um*1e-6;x=H*C/(lam*K*t)
 return 0. if x>700 else 2*H*C*C/(lam**5*math.expm1(x))*1e-6
photons=[dict(wavelength_nm=n,frequency_Hz=C/(n*1e-9),photon_energy_J=H*C/(n*1e-9),photon_energy_eV=H*C/(n*1e-9)/E,optical_power_W=1,photons_per_second=1/(H*C/(n*1e-9))) for n in range(380,751)]
rad=[]
for eps in [0,.2,.8,1]:
 for i in range(171):
  um=3+i*.1;b=planck(um,330);bg=planck(um,293)
  rad.append(dict(wavelength_um=um,surface_temperature_K=330,background_temperature_K=293,emissivity=eps,emitted_radiance_W_m2_sr_um=eps*b,reflected_background_radiance_W_m2_sr_um=(1-eps)*bg,observed_radiance_W_m2_sr_um=eps*b+(1-eps)*bg))
solar=[]
for area in [1,2,4]:
 for t in range(293,802,2):
  absorbed=800*area*.7*.9;radiative=.8*SIGMA*.2*(t**4-293**4);convective=5*.2*(t-293);balance=absorbed-radiative-convective
  solar.append(dict(projected_aperture_m2=area,receiver_temperature_K=t,absorbed_power_W=absorbed,radiative_loss_W=radiative,convective_loss_W=convective,residual_heat_W=balance,positive_extraction_available=balance>0))
for name,rows in [('photons.csv',photons),('radiance.csv',rad),('solar-balance.csv',solar)]:
 with (P/name).open('w',newline='') as f:w=csv.DictWriter(f,fieldnames=rows[0]);w.writeheader();w.writerows(rows)
assert all(math.isclose(r['photon_energy_J']*r['photons_per_second'],1,abs_tol=1e-14) for r in photons)
assert all(math.isclose(r['photon_energy_J'],H*r['frequency_Hz'],rel_tol=1e-14) for r in photons)
assert all(math.isclose(r['absorbed_power_W'],r['residual_heat_W']+r['radiative_loss_W']+r['convective_loss_W'],abs_tol=1e-10) for r in solar)
assert all(r['radiative_loss_W']==0 and r['convective_loss_W']==0 for r in solar if r['receiver_temperature_K']==293)
for eps in [0,.2,.8,1]:
 for um in [3,10,20]:assert math.isclose(eps*planck(um,293)+(1-eps)*planck(um,293),planck(um,293),rel_tol=1e-14)
# Broadband integration in wavelength, independent Stefan-Boltzmann identity.
wl=np.geomspace(.01,10000,50001);bb=np.array([planck(x,330) for x in wl]);integral=math.pi*np.trapezoid(bb,wl);relative=abs(integral/(SIGMA*330**4)-1)
assert relative<1e-6
peak_grid=np.linspace(7,11,40001);peak=float(peak_grid[np.argmax([planck(x,330) for x in peak_grid])]);assert abs(peak*330-2897.771955)<.04
negative_example=800*1*.7*.9-.8*SIGMA*.2*(800**4-293**4)-5*.2*(800-293);assert negative_example<0
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False,'svg.fonttype':'none','pdf.fonttype':42})
fig,ax=plt.subplots(1,2,figsize=(10.5,4.4))
ax[0].plot([r['wavelength_nm'] for r in photons],[r['photon_energy_eV'] for r in photons],color='#75469c',lw=2);ax[0].set(xlabel='Wavelength (nm)',ylabel='Energy per photon (eV)',title='A  |  Longer wavelength, less energy per photon')
ax[1].plot([r['wavelength_nm'] for r in photons],[r['photons_per_second']/1e18 for r in photons],color='#087e90',lw=2);ax[1].set(xlabel='Wavelength (nm)',ylabel='Photon rate (10^18 per second)',title='B  |  More photons at the same 1 W optical power')
for a in ax:a.grid(alpha=.2)
fig.suptitle('A photon energy is not a beam power',fontsize=15);fig.text(.5,.015,'Ideal monochromatic light. No material absorption, brightness, electrical efficiency or cutting performance is predicted.',ha='center',fontsize=9);fig.tight_layout(rect=[0,.05,1,.92])
for ext in ['png','svg','pdf']:fig.savefig(P/f'photon-ledger.{ext}',dpi=190)
fig,ax=plt.subplots(1,2,figsize=(11,4.8))
for eps,color in [(.2,'#087e90'),(.8,'#994994')]:
 rr=[r for r in rad if r['emissivity']==eps];x=[r['wavelength_um'] for r in rr]
 ax[0].plot(x,[r['emitted_radiance_W_m2_sr_um'] for r in rr],color=color,label=f'emitted, emissivity {eps}')
 ax[0].plot(x,[r['observed_radiance_W_m2_sr_um'] for r in rr],color=color,ls='--',label=f'with reflected background, {eps}')
ax[0].set(xlabel='Wavelength (micrometres)',ylabel='Spectral radiance (W m^-2 sr^-1 micrometre^-1)',title='A  |  Same surface temperature: 330 K');ax[0].legend(fontsize=8)
for area,color in [(1,'#087e90'),(2,'#994994'),(4,'#555')]:
 rr=[r for r in solar if r['projected_aperture_m2']==area];ax[1].plot([r['receiver_temperature_K'] for r in rr],[r['residual_heat_W'] for r in rr],color=color,label=f'{area} m² projected aperture')
ax[1].axhline(0,color='#888',lw=1);ax[1].set(xlabel='Prescribed receiver temperature (K)',ylabel='Residual heat balance (W)',title='B  |  Input minus assigned receiver losses');ax[1].legend(fontsize=8)
for a in ax:a.grid(alpha=.2)
fig.suptitle('Surfaces emit, reflect and lose heat',fontsize=15)
fig.text(.5,.052,'Hypothetical gray opaque diffuse surfaces; black surroundings at 293 K. Solar input 800 W/m², optical throughput 0.7, solar absorptivity 0.9.',ha='center',fontsize=8)
fig.text(.5,.014,'Receiver area 0.2 m², infrared emissivity 0.8, convection coefficient 5 W/(m² K). Residual is not a plant-output or equilibrium prediction.',ha='center',fontsize=8)
fig.tight_layout(rect=[0,.10,1,.92])
for ext in ['png','svg','pdf']:fig.savefig(P/f'surface-and-solar.{ext}',dpi=190)
meta={'title':'Photon accounting, gray-surface radiance and solar heat balance','status':'Analytical teaching study with hypothetical surface and collector parameters; no device performance measurements','rows':{'photons':len(photons),'radiance':len(rad),'solar':len(solar)},'constants':{'h_J_s':H,'c_m_s':C,'k_J_K':K,'elementary_charge_C':E,'derived_Stefan_Boltzmann_W_m2_K4':SIGMA},'assumptions':['Monochromatic optical power is fixed at 1 W; electrical input, efficacy and material interaction are not modeled.','Planck output is per wavelength in micrometres per area per solid angle, not integrated hemispheric flux.','Surface is opaque, diffuse and gray with wavelength-independent assigned emissivity; surroundings are black, isotropic and isothermal at 293 K; no atmospheric attenuation or instrument response.','Observed radiance includes reflected background and is not an inferred camera temperature.','Solar aperture is projected/effective collecting area, separate from receiver emitting area and land area. Direct irradiance is assigned 800 W/m².','Solar optical throughput is 0.7 and absorptivity is 0.9; infrared emissivity is separately 0.8. Solar and thermal spectral averages are not assumed equal.','Receiver temperature is prescribed; losses use view-factor-one large surroundings and assigned convection coefficient. Storage, conduction, piping and cycle conversion are excluded.','Negative residual means external heat would be needed to maintain that temperature with no extraction; it is not silently clamped. No annual yield, construction guidance or device performance follows.'],'checks':{'all_pass':True,'photon_power_and_frequency_identity':True,'solar_balance_all_rows':True,'equal_temperature_reflection_identity':True,'Planck_integrated_to_Stefan_Boltzmann_relative_error':relative,'blackbody_peak_um_at330K':peak,'negative_balance_example_W_at800K_1m2':negative_example},'formulas':{'photon':'E=h*c/lambda; rate=P/E','Planck_per_um':'2*h*c²/(lambda_m^5*expm1(h*c/(lambda_m*k*T)))*1e-6','observed':'epsilon*B(T)+(1-epsilon)*B(Tbackground)','solar':'800*Aprojected*0.7*0.9 - 0.8*sigma*0.2*(T^4-293^4) - 5*0.2*(T-293)'},'sources':[{'title':'NIST SI defining constants','url':'https://www.nist.gov/pml/special-publication-330/sp-330-section-2'},{'title':'NBS optical radiation measurements, Planck law','url':'https://nvlpubs.nist.gov/nistpubs/Legacy/TN/nbstechnicalnote910-8.pdf'},{'title':'NIST spectral emittance metrology','url':'https://www.nist.gov/laboratories/tools-instruments/system-infrared-spectral-emittance-materials'},{'title':'NASA thermal control','url':'https://www.nasa.gov/smallsat-institute/sst-soa/thermal-control/'},{'title':'DOE linear solar concentrators','url':'https://www.energy.gov/cmei/systems/linear-concentrator-system-concentrating-solar-thermal-power-basics'}]}
(P/'study.json').write_text(json.dumps(meta,indent=2)+'\n');print(json.dumps(meta['checks']))