Review index · Download raw source
reference/design/calculate_design_scenarios.py
1"""Reproducible receiver/optics design algebra, not measured device performance."""2from pathlib import Path3from datetime import datetime,timezone4import csv, hashlib,json,math,sys,time5import numpy as np6sys.path.insert(0,r'D:/Photometer_research_20261022/research/cloud_transport/private_readers')7from scipy.integrate import quad89REV=Path(__file__).resolve().parents[1];OUT=REV/'analysis/design';OUT.mkdir(exist_ok=True)10started=datetime.now(timezone.utc);tick=time.monotonic()11h=6.62607015e-34;c=299792458.;qe=1.602176634e-1912def table(name,rows):13 with (OUT/(name+'.csv')).open('w',newline='',encoding='utf-8') as f:14 w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)15def integrate_ramp(t,duration):16 a=np.maximum(t,0)17 return np.where(a<duration,a*a/(2*duration),a-duration/2)18def sample_rise(tr,rate,phase):19 dt=1/rate;D=tr*1e-6/.8;t0=phase*dt20 edges=np.arange(-2*dt,t0+D+4*dt,dt)21 y=np.diff(integrate_ramp(edges-t0,D))/dt;x=(edges[:-1]+edges[1:])/222 # Average samples are assigned their actual integration midpoints.23 t10=float(np.interp(.1,y,x));t90=float(np.interp(.9,y,x))24 return (t90-t10)*1e62526etendue=[]27for d in [1.,5.,25.]:28 for full in [30.,60.,80.,90.]:29 A=math.pi*(d*1e-3)**2/4;sa=math.sin(math.radians(full/2))30 for beta in [5.,10.,20.]:31 etendue.append(dict(aperture_diameter_mm=d,full_circular_field_deg=full,32 internal_filter_half_angle_deg=beta,projected_etendue_m2_sr=A*math.pi*sa**2,33 minimum_detector_area_mm2=A*1e6*sa**2,minimum_detector_diameter_mm=d*sa,34 minimum_filter_diameter_mm=d*sa/math.sin(math.radians(beta)),35 assumptions='Air, lossless passive optics, circular filled cones, necessary bounds only'))3637filters=[];qa=[]38for sw in [.5,1.,2.]:39 ss=sw/(2*math.sqrt(2*math.log(2)))40 for width in [1.,2.,4.,8.]:41 sf=width/(2*math.sqrt(2*math.log(2)));eq=sf*math.sqrt(2*math.pi)42 ref=(4/math.sqrt(sw*sw+16))/math.sqrt(4*math.sqrt(math.pi/(4*math.log(2))))43 for angle in range(0,31):44 neff=1.8;center=337.1*math.sqrt(1-math.sin(math.radians(angle))**2/neff**2)45 delta=center-337.146 cap=sf/math.sqrt(ss*ss+sf*sf)*math.exp(-delta*delta/(2*(ss*ss+sf*sf)))47 filters.append(dict(assumed_source_Gaussian_FWHM_nm=sw,filter_Gaussian_FWHM_nm=width,48 incidence_deg=angle,assumed_effective_index=neff,filter_center_nm=center,49 shift_nm=delta,equivalent_rectangular_band_nm=eq,unit_peak_capture_fraction=cap,50 conditional_background_dominated_SNR_relative_to_centered4nm=cap/math.sqrt(eq)/ref,51 source_model='Hypothetical normalized Gaussian, NOT a fitted N2 band'))52 if angle in [0,5,10,20]:53 num=quad(lambda x:math.exp(-x*x/(2*ss*ss))/(ss*math.sqrt(2*math.pi))*54 math.exp(-(x-delta)**2/(2*sf*sf)),-20,20,epsabs=1e-12)[0]55 assert abs(num-cap)<1e-956 qa.append(abs(num-cap))5758sampling=[]59for tr in [30.,80.,200.]:60 for rate in [10000,50000,100000]:61 for phase in np.linspace(0,1,101)[:-1]:62 sampling.append(dict(true_10_90_us=tr,sample_rate_s=rate,phase_fraction=float(phase),63 noiseless_reconstructed_10_90_us=sample_rise(tr,rate,float(phase)),64 sample_definition='Exact average over nonoverlapping bins, linearly interpolated at bin midpoints'))6566droop=[]67for duration in [.003,.005,.01]:68 for f in [.01,.05,.1]:69 tau=-duration/math.log1p(-f)70 droop.append(dict(plateau_duration_s=duration,maximum_fractional_droop=f,71 minimum_highpass_timeconstant_s=tau,maximum_highpass_corner_Hz=1/(2*math.pi*tau),72 model='Single-pole AC coupling, electrical step plateau'))7374blocking=[]75for ratio in np.logspace(2,9,29):76 for allowed in [.01,.1]:77 blocking.append(dict(assumed_unblocked_detector_weighted_out_to_in_background_ratio=float(ratio),78 allowed_leakage_fraction_of_inband_background=allowed,required_uniform_blocking_OD=math.log10(ratio/allowed),79 model='Allocation algebra; actual source spectra, filter curve and detector response required'))8081nep=[]82for wave in [337.1,777.4]:83 E=h*c/(wave*1e-9)84 for eff in [.05,.2,.5]:85 for dark in [1.,50.,1000.,1e6]:86 for F in [1.,1.5,2.]:87 nep.append(dict(wavelength_nm=wave,assumed_end_to_end_efficiency=eff,88 assumed_dark_primary_cps=dark,assumed_excess_noise_factor=F,89 dark_only_entrance_referred_NEP_W_sqrtHz=E*math.sqrt(2*F*dark)/eff,90 model='Ideal primary-count shot noise; one-sided PSD, no sky, amplifier, afterpulse or radiation forecast'))9192count_tests=[]93for wave in [337.1,777.4]:94 E=h*c/(wave*1e-9)95 for d in [1.,5.,25.]:96 A=math.pi*(d*.001)**2/497 for full in [60.,90.]:98 G=A*math.pi*math.sin(math.radians(full/2))**299 for L in [0.,1.77e-9,3e-9,3e-7]:100 for dt in [30e-6,80e-6,200e-6,.003]:101 sky=.2*G*4*L/E102 count_tests.append(dict(wavelength_nm=wave,aperture_diameter_mm=d,full_circular_FOV_deg=full,103 efficiency=.2,assumed_continuum_W_m2_sr_nm=L,equivalent_band_nm=4,window_s=dt,104 sky_primary_cps=sky,dark_primary_cps=50,expected_background_count=(sky+50)*dt,105 input_aperture_fluence_ph_m2=1e8,registered_signal_count_for_input_fluence=.2*A*1e8))106107for name,rows in [('etendue',etendue),('parametric_filter',filters),('sampling_phase',sampling),108 ('ac_coupling',droop),('blocking_allocation',blocking),('dark_only_nep',nep),109 ('browser_budget_reference',count_tests)]:table(name,rows)110result={'schema':'DESIGN_CONDITIONAL_ALGEBRA_R4','created_utc':datetime.now(timezone.utc).isoformat(),111 'started_utc':started.isoformat(),'wall_seconds':time.monotonic()-tick,112 'evidence_state':'EXECUTED_CONDITIONAL_CALCULATIONS_NOT_DEVICE_MEASUREMENTS',113 'constants':{'h_J_s':h,'c_m_s':c,'electron_charge_C':qe,114 'photon_energy_337p1_J':h*c/(337.1e-9),'photon_energy_777p4_J':h*c/(777.4e-9),115 'one_Rayleigh_ph_m2_sr_s':1e10/(4*math.pi)},116 'etendue':etendue,'parametric_filter':filters,'sampling_phase':sampling,'ac_coupling':droop,117 'blocking_allocation':blocking,'dark_only_nep':nep,'browser_budget_reference':count_tests,118 'QA':{'Gaussian_overlap_independent_quadratures':len(qa),'max_abs_capture_error':max(qa),119 'ramp_bin_integration_exact':True,'rows':{name:len(rows) for name,rows in [120 ('etendue',etendue),('parametric_filter',filters),('sampling_phase',sampling),121 ('ac_coupling',droop),('blocking_allocation',blocking),('dark_only_nep',nep),('browser_budget_reference',count_tests)]}},122 'limits':[123 'Gaussian source widths are hypothetical design inputs. They do not repair unvalidated nitrogen line data or establish an optimal flight filter.',124 'Effective-index angle formula is a coating approximation; real angular, temperature and polarization curves must be measured.',125 'Aperture fluence already includes atmospheric delivery. Do not apply cloud transmission a second time to observed receiver fluence.',126 'A noise-equivalent-power value here is an ideal dark-only model, not a detector-class specification or measured end-of-life receiver NEP.',127 'Sampling curves are noise-free averaging/interpolation examples, not a population depth classifier.',128 'Étendue and filter bounds are necessary conditions, not a ray-traced1U packaging solution.',129 'BlockingOD is required over the actual detector-weighted leak range; no universal OD6 guarantee.'130 ],131 'sources':[132 {'id':'MMIA2019','url':'https://link.springer.com/article/10.1007/s11214-019-0593-y','scope':'Optical relay, aperture/FOV and count architecture'},133 {'id':'ETOUTPUT','url':'https://et-enterprises.com/images/technical_papers/rp065_design_of_PMT_output_circuits.pdf','scope':'DC/AC PMT readout design'},134 {'id':'FEMTO','url':'https://www.femto.de/images/pdf-dokumente/de-dlpca-200.pdf','scope':'Actual amplifier gain/bandwidth/noise need their stated conditions'}135 ]}136(OUT/'design_results.json').write_text(json.dumps(result,indent=2,allow_nan=False),encoding='utf-8')137(OUT/'execution_receipt.json').write_text(json.dumps({'status':'PASS','created_utc':datetime.now(timezone.utc).isoformat(),138 'command':[sys.executable,'-B',str(Path(__file__).resolve())],'script_sha256':hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),139 'rows':result['QA']['rows'],'independent_Gaussian_checks':len(qa),'max_error':max(qa),'wall_seconds':result['wall_seconds']},indent=2),encoding='utf-8')140print(json.dumps({'status':'PASS','rows':result['QA']['rows'],'max_error':max(qa)},indent=2))141