"""Reproducible receiver/optics design algebra, not measured device performance."""
from pathlib import Path
from datetime import datetime,timezone
import csv, hashlib,json,math,sys,time
import numpy as np
sys.path.insert(0,r'D:/Photometer_research_20261022/research/cloud_transport/private_readers')
from scipy.integrate import quad

REV=Path(__file__).resolve().parents[1];OUT=REV/'analysis/design';OUT.mkdir(exist_ok=True)
started=datetime.now(timezone.utc);tick=time.monotonic()
h=6.62607015e-34;c=299792458.;qe=1.602176634e-19
def table(name,rows):
    with (OUT/(name+'.csv')).open('w',newline='',encoding='utf-8') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
def integrate_ramp(t,duration):
    a=np.maximum(t,0)
    return np.where(a<duration,a*a/(2*duration),a-duration/2)
def sample_rise(tr,rate,phase):
    dt=1/rate;D=tr*1e-6/.8;t0=phase*dt
    edges=np.arange(-2*dt,t0+D+4*dt,dt)
    y=np.diff(integrate_ramp(edges-t0,D))/dt;x=(edges[:-1]+edges[1:])/2
    # Average samples are assigned their actual integration midpoints.
    t10=float(np.interp(.1,y,x));t90=float(np.interp(.9,y,x))
    return (t90-t10)*1e6

etendue=[]
for d in [1.,5.,25.]:
 for full in [30.,60.,80.,90.]:
  A=math.pi*(d*1e-3)**2/4;sa=math.sin(math.radians(full/2))
  for beta in [5.,10.,20.]:
   etendue.append(dict(aperture_diameter_mm=d,full_circular_field_deg=full,
    internal_filter_half_angle_deg=beta,projected_etendue_m2_sr=A*math.pi*sa**2,
    minimum_detector_area_mm2=A*1e6*sa**2,minimum_detector_diameter_mm=d*sa,
    minimum_filter_diameter_mm=d*sa/math.sin(math.radians(beta)),
    assumptions='Air, lossless passive optics, circular filled cones, necessary bounds only'))

filters=[];qa=[]
for sw in [.5,1.,2.]:
 ss=sw/(2*math.sqrt(2*math.log(2)))
 for width in [1.,2.,4.,8.]:
  sf=width/(2*math.sqrt(2*math.log(2)));eq=sf*math.sqrt(2*math.pi)
  ref=(4/math.sqrt(sw*sw+16))/math.sqrt(4*math.sqrt(math.pi/(4*math.log(2))))
  for angle in range(0,31):
   neff=1.8;center=337.1*math.sqrt(1-math.sin(math.radians(angle))**2/neff**2)
   delta=center-337.1
   cap=sf/math.sqrt(ss*ss+sf*sf)*math.exp(-delta*delta/(2*(ss*ss+sf*sf)))
   filters.append(dict(assumed_source_Gaussian_FWHM_nm=sw,filter_Gaussian_FWHM_nm=width,
    incidence_deg=angle,assumed_effective_index=neff,filter_center_nm=center,
    shift_nm=delta,equivalent_rectangular_band_nm=eq,unit_peak_capture_fraction=cap,
    conditional_background_dominated_SNR_relative_to_centered4nm=cap/math.sqrt(eq)/ref,
    source_model='Hypothetical normalized Gaussian, NOT a fitted N2 band'))
   if angle in [0,5,10,20]:
    num=quad(lambda x:math.exp(-x*x/(2*ss*ss))/(ss*math.sqrt(2*math.pi))*
             math.exp(-(x-delta)**2/(2*sf*sf)),-20,20,epsabs=1e-12)[0]
    assert abs(num-cap)<1e-9
    qa.append(abs(num-cap))

sampling=[]
for tr in [30.,80.,200.]:
 for rate in [10000,50000,100000]:
  for phase in np.linspace(0,1,101)[:-1]:
   sampling.append(dict(true_10_90_us=tr,sample_rate_s=rate,phase_fraction=float(phase),
    noiseless_reconstructed_10_90_us=sample_rise(tr,rate,float(phase)),
    sample_definition='Exact average over nonoverlapping bins, linearly interpolated at bin midpoints'))

droop=[]
for duration in [.003,.005,.01]:
 for f in [.01,.05,.1]:
  tau=-duration/math.log1p(-f)
  droop.append(dict(plateau_duration_s=duration,maximum_fractional_droop=f,
    minimum_highpass_timeconstant_s=tau,maximum_highpass_corner_Hz=1/(2*math.pi*tau),
    model='Single-pole AC coupling, electrical step plateau'))

blocking=[]
for ratio in np.logspace(2,9,29):
 for allowed in [.01,.1]:
  blocking.append(dict(assumed_unblocked_detector_weighted_out_to_in_background_ratio=float(ratio),
    allowed_leakage_fraction_of_inband_background=allowed,required_uniform_blocking_OD=math.log10(ratio/allowed),
    model='Allocation algebra; actual source spectra, filter curve and detector response required'))

nep=[]
for wave in [337.1,777.4]:
 E=h*c/(wave*1e-9)
 for eff in [.05,.2,.5]:
  for dark in [1.,50.,1000.,1e6]:
   for F in [1.,1.5,2.]:
    nep.append(dict(wavelength_nm=wave,assumed_end_to_end_efficiency=eff,
     assumed_dark_primary_cps=dark,assumed_excess_noise_factor=F,
     dark_only_entrance_referred_NEP_W_sqrtHz=E*math.sqrt(2*F*dark)/eff,
     model='Ideal primary-count shot noise; one-sided PSD, no sky, amplifier, afterpulse or radiation forecast'))

count_tests=[]
for wave in [337.1,777.4]:
 E=h*c/(wave*1e-9)
 for d in [1.,5.,25.]:
  A=math.pi*(d*.001)**2/4
  for full in [60.,90.]:
   G=A*math.pi*math.sin(math.radians(full/2))**2
   for L in [0.,1.77e-9,3e-9,3e-7]:
    for dt in [30e-6,80e-6,200e-6,.003]:
     sky=.2*G*4*L/E
     count_tests.append(dict(wavelength_nm=wave,aperture_diameter_mm=d,full_circular_FOV_deg=full,
      efficiency=.2,assumed_continuum_W_m2_sr_nm=L,equivalent_band_nm=4,window_s=dt,
      sky_primary_cps=sky,dark_primary_cps=50,expected_background_count=(sky+50)*dt,
      input_aperture_fluence_ph_m2=1e8,registered_signal_count_for_input_fluence=.2*A*1e8))

for name,rows in [('etendue',etendue),('parametric_filter',filters),('sampling_phase',sampling),
                  ('ac_coupling',droop),('blocking_allocation',blocking),('dark_only_nep',nep),
                  ('browser_budget_reference',count_tests)]:table(name,rows)
result={'schema':'DESIGN_CONDITIONAL_ALGEBRA_R4','created_utc':datetime.now(timezone.utc).isoformat(),
 'started_utc':started.isoformat(),'wall_seconds':time.monotonic()-tick,
 'evidence_state':'EXECUTED_CONDITIONAL_CALCULATIONS_NOT_DEVICE_MEASUREMENTS',
 'constants':{'h_J_s':h,'c_m_s':c,'electron_charge_C':qe,
  'photon_energy_337p1_J':h*c/(337.1e-9),'photon_energy_777p4_J':h*c/(777.4e-9),
  'one_Rayleigh_ph_m2_sr_s':1e10/(4*math.pi)},
 'etendue':etendue,'parametric_filter':filters,'sampling_phase':sampling,'ac_coupling':droop,
 'blocking_allocation':blocking,'dark_only_nep':nep,'browser_budget_reference':count_tests,
 'QA':{'Gaussian_overlap_independent_quadratures':len(qa),'max_abs_capture_error':max(qa),
       'ramp_bin_integration_exact':True,'rows':{name:len(rows) for name,rows in [
       ('etendue',etendue),('parametric_filter',filters),('sampling_phase',sampling),
       ('ac_coupling',droop),('blocking_allocation',blocking),('dark_only_nep',nep),('browser_budget_reference',count_tests)]}},
 'limits':[
  'Gaussian source widths are hypothetical design inputs. They do not repair unvalidated nitrogen line data or establish an optimal flight filter.',
  'Effective-index angle formula is a coating approximation; real angular, temperature and polarization curves must be measured.',
  'Aperture fluence already includes atmospheric delivery. Do not apply cloud transmission a second time to observed receiver fluence.',
  'A noise-equivalent-power value here is an ideal dark-only model, not a detector-class specification or measured end-of-life receiver NEP.',
  'Sampling curves are noise-free averaging/interpolation examples, not a population depth classifier.',
  'Étendue and filter bounds are necessary conditions, not a ray-traced1U packaging solution.',
  'BlockingOD is required over the actual detector-weighted leak range; no universal OD6 guarantee.'
 ],
 'sources':[
  {'id':'MMIA2019','url':'https://link.springer.com/article/10.1007/s11214-019-0593-y','scope':'Optical relay, aperture/FOV and count architecture'},
  {'id':'ETOUTPUT','url':'https://et-enterprises.com/images/technical_papers/rp065_design_of_PMT_output_circuits.pdf','scope':'DC/AC PMT readout design'},
  {'id':'FEMTO','url':'https://www.femto.de/images/pdf-dokumente/de-dlpca-200.pdf','scope':'Actual amplifier gain/bandwidth/noise need their stated conditions'}
 ]}
(OUT/'design_results.json').write_text(json.dumps(result,indent=2,allow_nan=False),encoding='utf-8')
(OUT/'execution_receipt.json').write_text(json.dumps({'status':'PASS','created_utc':datetime.now(timezone.utc).isoformat(),
 'command':[sys.executable,'-B',str(Path(__file__).resolve())],'script_sha256':hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
 'rows':result['QA']['rows'],'independent_Gaussian_checks':len(qa),'max_error':max(qa),'wall_seconds':result['wall_seconds']},indent=2),encoding='utf-8')
print(json.dumps({'status':'PASS','rows':result['QA']['rows'],'max_error':max(qa)},indent=2))