Review index · Download raw source
modules/design_review.py
1"""Independent SI design algebra and ideal-model diagnostics.23This module does not import the original calculation scripts. Expected tables are4loaded only as comparison targets and as declared input-parameter grids. The5count-threshold S/n columns in the FOV and auroral tables are deliberately inputs:6their exact Poisson validation belongs to the statistics verification module.7Run from any directory: python modules/design_review.py8"""9from pathlib import Path10from datetime import datetime, timezone11from functools import lru_cache12import ast, hashlib, json, math, platform, time13import numpy as np1415ROOT = Path(__file__).resolve().parents[1]16H = 6.62607015e-3417C = 299792458.018Q_E = 1.602176634e-1919K_B = 1.380649e-2320RAYLEIGH_PHOTONS = 1e10 / (4 * math.pi)2122def load(name):23 return json.loads((ROOT / 'data' / name).read_text(encoding='utf-8'))2425@lru_cache(None)26def nodes(order):27 return np.polynomial.legendre.leggauss(order)2829def integrate(fun, lo, hi, order=192):30 """Fixed Gauss-Legendre integral, including an explicit interval Jacobian."""31 if hi <= lo:32 return 0.033 x, w = nodes(order)34 a = (hi - lo) / 235 return float(a * np.dot(w, fun((hi + lo) / 2 + a * x)))3637def bisect(fun, lo, hi, iterations=72):38 """Bracketed root; no original root-solver implementation is reused."""39 flo, fhi = fun(lo), fun(hi)40 if flo == 0: return lo41 if fhi == 0: return hi42 if flo * fhi > 0: raise ValueError(('unbracketed', lo, hi, flo, fhi))43 for _ in range(iterations):44 mid = (lo + hi) / 245 fm = fun(mid)46 if fm == 0: return mid47 if flo * fm <= 0: hi = mid48 else: lo, flo = mid, fm49 return (lo + hi) / 25051def photon_energy(wavelength_m):52 """D01: energy in J for a wavelength explicitly supplied in metres."""53 if wavelength_m <= 0: raise ValueError('wavelength must be positive')54 return H * C / wavelength_m5556def circular_angles(half_angle_rad):57 """D02: ordinary solid angle and projected angular integral, in sr."""58 return (4 * math.pi * math.sin(half_angle_rad / 2)**2,59 math.pi * math.sin(half_angle_rad)**2)6061def square_angles(diagonal_half_angle_rad):62 """D03: square field; its diagonal half-angle is not a circular half-angle."""63 u = math.tan(diagonal_half_angle_rad) / math.sqrt(2)64 omega = 4 * math.atan(u*u / math.sqrt(1 + 2*u*u))65 x, w = nodes(96); xx = u*x66 projected = float(u*u*np.einsum('i,j,ij->', w, w,67 1 / (1 + xx[:, None]**2 + xx[None, :]**2)**2))68 return omega, projected6970def aperture_area(diameter_m):71 return math.pi * diameter_m**2 / 47273def counts_from_received_fluence(fluence_ph_m2, area_m2, efficiency, incidence_rad=0.):74 """D04: receiver-plane photon fluence -> detected expected counts."""75 return fluence_ph_m2 * area_m2 * efficiency * math.cos(incidence_rad)7677def continuum_rate(L_W_m2_sr_per_nm, band_nm, area_m2, half_angle_rad,78 efficiency, wavelength_m):79 """D05: uniform energy radiance, rectangular/equivalent band, flat response."""80 return (L_W_m2_sr_per_nm * band_nm * area_m2 *81 circular_angles(half_angle_rad)[1] * efficiency / photon_energy(wavelength_m))8283def etendue_bounds(diameter_m, field_half_rad, filter_half_rad):84 """D06: necessary air-to-air phase-space bounds; no real optical prescription."""85 a = aperture_area(diameter_m)86 g = a * circular_angles(field_half_rad)[1]87 detector = g / math.pi # detector acceptance is an entire air hemisphere88 filt = g / (math.pi * math.sin(filter_half_rad)**2)89 return g, detector, math.sqrt(4*detector/math.pi), math.sqrt(4*filt/math.pi)9091def shifted_center(wavelength_m, incidence_rad, effective_index):92 """D07: scalar effective-index coating approximation, not measured coating data."""93 return wavelength_m * math.sqrt(1 - (math.sin(incidence_rad)/effective_index)**2)9495def sigma_from_fwhm(fwhm_m):96 return fwhm_m / math.sqrt(8 * math.log(2))9798def gaussian_overlap(source_sigma_m, filter_sigma_m, shift_m):99 """D08: normalized photon-source Gaussian times a unit-peak filter Gaussian."""100 variance = source_sigma_m**2 + filter_sigma_m**2101 return filter_sigma_m / math.sqrt(variance) * math.exp(-shift_m**2/(2*variance))102103def gaussian_overlap_integral(source_sigma_m, filter_sigma_m, shift_m):104 """D08 independent quadrature in standard-normal source coordinates."""105 return integrate(lambda z: np.exp(-z*z/2)/math.sqrt(2*math.pi) *106 np.exp(-(source_sigma_m*z-shift_m)**2/(2*filter_sigma_m**2)), -12., 12., 256)107108def gaussian_equivalent_band(filter_sigma_m):109 """D09: spectral equivalent rectangular width in m, distinct from electrical ENBW."""110 return math.sqrt(2 * math.pi) * filter_sigma_m111112def electrical_enbw_first_order(corner_Hz):113 """D10: one-sided ENBW of H(f)=1/(1+i*f/fc)."""114 return math.pi * corner_Hz / 2115116def boxcar_enbw(exposure_s):117 """D10: one-sided ENBW for a unit-DC-gain rectangular average."""118 return 1 / (2 * exposure_s)119120def blocking_od(unblocked_out_to_in_ratio, allowed_leak_fraction):121 """D11: uniform blocking bound on a detector-weighted photon-rate ratio."""122 return max(0., math.log10(unblocked_out_to_in_ratio/allowed_leak_fraction))123124def dark_nep(wavelength_m, efficiency, primary_dark_s, excess_noise_factor=1.):125 """D12: entrance-referred one-sided amplitude spectral density, W/sqrt(Hz)."""126 return photon_energy(wavelength_m)/efficiency * math.sqrt(2*excess_noise_factor*primary_dark_s)127128def hp_droop_timeconstant(plateau_s, maximum_droop):129 """D13: electrical step y=e^(-t/tau); return minimum tau and maximum fc."""130 tau = plateau_s / (-math.log1p(-maximum_droop))131 return tau, 1/(2*math.pi*tau)132133def ramp_bin_average(left_s, right_s, onset_s, ramp_s):134 """D14: integrate ramp by separating its linear piece and plateau rectangle."""135 a, b = left_s-onset_s, right_s-onset_s136 lower, upper = max(a, 0.), min(b, ramp_s)137 triangle = (upper-lower)*(upper+lower)/(2*ramp_s) if upper>lower else 0.138 rectangle = max(0., b-max(a, ramp_s))139 return (triangle + rectangle)/(right_s-left_s)140141def sampled_rise(rise_s, sample_rate_s, phase):142 """D14: bin-average linear ramp, then interpolate crossings at exposure midpoints."""143 dt = 1/sample_rate_s; ramp = rise_s/.8; onset = phase*dt144 edges = np.arange(-2, math.ceil((onset+ramp)/dt)+5, dtype=float)*dt145 y = np.array([ramp_bin_average(a,b,onset,ramp) for a,b in zip(edges[:-1],edges[1:])])146 x = (edges[:-1]+edges[1:])/2147 def crossing(level):148 j = int(np.nonzero(y >= level)[0][0])149 return x[j-1] + (x[j]-x[j-1])*(level-y[j-1])/(y[j]-y[j-1])150 return crossing(.9)-crossing(.1)151152class ReferenceAtmosphere:153 """P01: documented USSA1976 layer constants; no new polar atmosphere is invented."""154 R = 6356766.; g = 9.80665; gas = 8314.32; M = 28.9644155 levels = np.array([0.,11000.,20000.,32000.,47000.,51000.,71000.,84852.])156 lapse = np.array([-.0065,0.,.001,.0028,0.,-.0028,-.002])157 def __init__(self):158 self.T=[288.15]; self.p=[101325.]159 for j,L in enumerate(self.lapse):160 dh=self.levels[j+1]-self.levels[j]; t=self.T[-1]161 self.T.append(t+L*dh)162 log_ratio=(-self.g*self.M*dh/(self.gas*t) if L==0 else163 -self.g*self.M/(self.gas*L)*math.log1p(L*dh/t))164 self.p.append(self.p[-1]*math.exp(log_ratio))165 def geopotential(self, z_m): return self.R*z_m/(self.R+z_m)166 def state_H(self, H_m):167 i=min(6,max(0,int(np.searchsorted(self.levels,H_m,side='right')-1)))168 L=self.lapse[i];delta=H_m-self.levels[i]; T=self.T[i]+L*delta169 logp=(-self.g*self.M*delta/(self.gas*self.T[i]) if L==0 else170 -self.g*self.M/(self.gas*L)*math.log1p(L*delta/self.T[i]))171 return self.p[i]*math.exp(logp),T172 def state(self,z_m): return self.state_H(self.geopotential(z_m))173 def pressure_column(self,z_m,order=96):174 """Independent pressure-coordinate quadrature; original used density trapezoids."""175 H0=self.geopotential(z_m); total=0.176 for i,L in enumerate(self.lapse):177 lo=max(H0,self.levels[i]); hi=self.levels[i+1]178 if lo>=hi: continue179 pl,_=self.state_H(lo); ph,_=self.state_H(hi)180 def jac(p):181 if L==0:182 hh=self.levels[i]-self.gas*self.T[i]/(self.g*self.M)*np.log(p/self.p[i])183 else:184 hh=self.levels[i]+self.T[i]/L*np.expm1(-self.gas*L/(self.g*self.M)*np.log(p/self.p[i]))185 return 1/(1-hh/self.R)**2186 total+=integrate(jac,ph,pl,order)187 # Retain and expose the SAME small pressure-tail approximation at ~86 km.188 tail=self.p[-1]/(1-self.levels[-1]/self.R)**2189 constant=self.gas/(self.g*self.M*K_B)190 return constant*(total+tail), constant*tail191192def rayleigh_bodhaine(wavelength_m):193 """P02: Bodhaine1999 Eq29 fit; wavelength argument to rational fit is micrometres."""194 u=wavelength_m*1e6195 r=(1.0455996-341.29061/u**2-.90230850*u*u)/(1+.0027059889/u**2-85.968563*u*u)196 return r*1e-32, .0021520*r197198def rayleigh_frozen(wavelength_m):199 """P02: separate frozen CloudScat source arithmetic, not a replacement atmosphere."""200 u=wavelength_m*1e6; n0=2.6867811e25201 n=1+1e-8*(8060.51+2480990/(132.274-1/u**2)+17455.7/(39.32957-1/u**2))202 king=(780840*(1.034+3.17e-4/u**2)+209476*(1.096+1.385e-3/u**2+1.448e-4/u**4)+9340+314*1.15)/1e6203 sigma=24*math.pi**3*((n*n-1)/(n*n+2))**2*king/(wavelength_m**4*(n0*273.15/288.15)**2)204 return sigma,n0,7200.205206def direct_fraction(area_m2, distance_m, ray_tau=0., cloud_tau=0.):207 """P03: ballistic photons per emitted isotropic photon, not total cloudy light."""208 return area_m2/(4*math.pi*distance_m**2)*math.exp(-ray_tau-cloud_tau)209210class LambertPlane:211 """P04-P05: ideal instantaneous Lambert surface, expressed in log radial coordinate."""212 def __init__(self, observer_height_m, source_height_m, albedo, field_half_rad):213 self.H=observer_height_m;self.h=source_height_m;self.rho=albedo214 if not 0<=source_height_m<observer_height_m:raise ValueError('0 <= h < H required')215 self.radius=observer_height_m*math.tan(field_half_rad)216 self.qmax=math.log1p(self.radius/source_height_m) if source_height_m else 0.217 self.B=self.integral(0,self.qmax) if source_height_m else 2*albedo218 def density(self,q):219 u=np.expm1(q); a=self.h/self.H220 return 2*self.rho*(1-a)**2*u*np.exp(q)/(np.hypot(1,u)**3*(1+(a*u)**2)**2)221 def delay(self,q):222 u=np.expm1(q); x=self.h*u/self.H223 return (self.h*(np.hypot(1,u)+1)+self.H*x*x/(np.hypot(1,x)+1))/C224 def integral(self,a,b): return integrate(self.density,a,b,192)225 def cumulative_q(self,t):226 if t<=self.delay(0.):return 0.227 if t>=self.delay(self.qmax):return self.qmax228 return bisect(lambda q:float(self.delay(q))-t,0.,self.qmax)229 def received_ramp(self,t,ramp_s):230 if not self.h:return float(np.clip(t/ramp_s,0,1))231 a=self.cumulative_q(t-ramp_s);b=self.cumulative_q(t)232 reflected=self.integral(0,a)+integrate(lambda q:self.density(q)*(t-self.delay(q))/ramp_s,a,b,128)233 return (float(np.clip(t/ramp_s,0,1))+reflected)/(1+self.B)234 def rise(self,input_10_90_s):235 if not self.h or not self.rho:return input_10_90_s236 ramp=input_10_90_s/.8;hi=ramp+float(self.delay(self.qmax))237 lo=bisect(lambda t:self.received_ramp(t,ramp)-.1,0.,hi)238 up=bisect(lambda t:self.received_ramp(t,ramp)-.9,0.,hi)239 return up-lo240241def diffusion_kernel(t_s, depth_m, D_m2_s, gamma_s=0.):242 """P06: first-hitting density s^-1; an acausal diffusion approximation, not early RTE."""243 t=np.asarray(t_s,dtype=float);safe=np.maximum(t,np.finfo(float).tiny)244 with np.errstate(over='ignore',under='ignore'):245 log_density=math.log(depth_m)-.5*math.log(4*math.pi*D_m2_s)-1.5*np.log(safe)-depth_m**2/(4*D_m2_s*safe)-gamma_s*safe246 result=np.exp(log_density)247 return np.where(t>0,result,0.)248249def diffusion_peak(depth_m,D_m2_s,gamma_s=0.):250 a=depth_m**2/(4*D_m2_s)251 return 2*a/(1.5+math.sqrt(2.25+4*gamma_s*a))252253def diffusion_rise(depth_m,D_m2_s,gamma_s=0.):254 peak=diffusion_peak(depth_m,D_m2_s,gamma_s); a=depth_m**2/(4*D_m2_s)255 def log_ratio(t):return -1.5*math.log(t/peak)-a*(1/t-1/peak)-gamma_s*(t-peak)256 t10=bisect(lambda t:log_ratio(t)-math.log(.1),peak*1e-9,peak)257 t90=bisect(lambda t:log_ratio(t)-math.log(.9),peak*1e-9,peak)258 return t10,t90,peak259260def diffusion_exit_integral(depth_m,D_m2_s,gamma_s=0.,window_s=None):261 """P06 independent log-u quadrature of h(t), avoiding the original scipy.integrate route."""262 a=gamma_s*depth_m**2/(4*D_m2_s)263 lower=-32. if window_s is None else math.log(depth_m/math.sqrt(4*D_m2_s*window_s))264 return integrate(lambda v:2/math.sqrt(math.pi)*np.exp(v-np.exp(2*v)-a*np.exp(-2*v)),lower,math.log(12.),512)265266def diffusion_cdf(window_s,depth_m,D_m2_s,gamma_s=0.):267 """P06 closed finite-time exit probability, evaluated in this bounded grid."""268 a=depth_m/math.sqrt(4*D_m2_s*window_s)269 q=math.sqrt(gamma_s*window_s);b=depth_m*math.sqrt(gamma_s/D_m2_s)270 return .5*(math.exp(-b)*math.erfc(a-q)+math.exp(b)*math.erfc(a+q))271272def deadtime_mean(rate_s,deadtime_s):273 """P07: stationary nonparalyzable renewal mean, no Poisson output assertion."""274 return 1/(1/rate_s+deadtime_s) if rate_s else 0.275276def deadtime_fractions(background_s,signal_s,deadtime_s):277 """P07 inverse mean-rate route, independent of original closed shape formula."""278 mb=deadtime_mean(background_s,deadtime_s);mp=deadtime_mean(background_s+signal_s,deadtime_s)279 values=[]280 for y in [.1,.9]:281 measured=mb+y*(mp-mb);true=measured/(1-measured*deadtime_s)282 values.append((true-background_s)/signal_s)283 return values[0],values[1],(values[1]-values[0])/.8284285def spherical_field(half_rad,observer_radius_m,cloud_radius_m):286 """P08: central cap angle and cap area for the near shell intersection."""287 psi=math.asin(observer_radius_m/cloud_radius_m*math.sin(half_rad))-half_rad288 return psi,4*math.pi*cloud_radius_m**2*math.sin(psi/2)**2289290def shell_geometry(psi,observer_radius_m,cloud_radius_m):291 """P08: independent shell-central-angle coordinates, not observer-ray quadrature."""292 s2=np.sin(np.asarray(psi)/2)**2;h=observer_radius_m-cloud_radius_m293 distance=np.sqrt(h*h+4*observer_radius_m*cloud_radius_m*s2)294 cos_view=(h+2*cloud_radius_m*s2)/distance295 return distance,cos_view296297def population_weight(required_detected_count,area_m2,efficiency,half_rad,alpha,298 observer_radius_m=6921000.,cloud_radius_m=6388000.,pivot_ph=1e20):299 """P09: hypothetical bright-enough weight, not a full detection integral or rate."""300 end,_=spherical_field(half_rad,observer_radius_m,cloud_radius_m)301 def integrand(psi):302 distance,cv=shell_geometry(psi,observer_radius_m,cloud_radius_m)303 n=required_detected_count*4*math.pi*distance**2/(efficiency*area_m2*cv)304 return (n/pivot_ph)**(-alpha)*2*math.pi*cloud_radius_m**2*np.sin(psi)305 return integrate(integrand,0.,end,192)306307class Checks:308 def __init__(self):self.fields={};self.failures=[];self.comparisons=0309 def add(self,group,row,column,fresh,expected,rtol=2e-10,atol=0.):310 key=group+'.'+column311 x=self.fields.setdefault(key,{'comparisons':0,'max_absolute_error':0.,'max_relative_error':0.})312 for label,value in [('relative_tolerance',rtol),('absolute_tolerance_in_column_units',atol)]:313 interval=x.setdefault(label+'_range',[value,value]);interval[0]=min(interval[0],value);interval[1]=max(interval[1],value)314 x['comparisons']+=1;self.comparisons+=1315 if not math.isfinite(float(fresh)) or not math.isfinite(float(expected)):316 self.failures.append({'group':group,'row':row,'column':column,'fresh':repr(fresh),'expected':repr(expected),'reason':'nonfinite value rejected'})317 return318 err=abs(float(fresh)-float(expected));x['max_absolute_error']=max(x['max_absolute_error'],err)319 if float(expected)!=0:x['max_relative_error']=max(x['max_relative_error'],err/abs(float(expected)))320 else:x['zero_reference_comparisons']=x.get('zero_reference_comparisons',0)+1321 if err>atol+rtol*abs(float(expected)):self.failures.append({'group':group,'row':row,'column':column,'fresh':fresh,'expected':expected,'absolute_error':err,'rtol':rtol,'atol':atol})322 def result(self):return {'status':'PASS' if not self.failures else 'FAIL','comparisons':self.comparisons,'fields':self.fields,'failures':self.failures}323324def check_design(write=True):325 d=load('design_results.json');t=Checks();fresh={};quad=[]326 for wave,key in [(337.1,'photon_energy_337p1_J'),(777.4,'photon_energy_777p4_J')]:t.add('constants',wave,key,photon_energy(wave*1e-9),d['constants'][key])327 t.add('constants',0,'one_Rayleigh_ph_m2_sr_s',RAYLEIGH_PHOTONS,d['constants']['one_Rayleigh_ph_m2_sr_s'])328 for name,rows in [(k,d[k]) for k in ['etendue','parametric_filter','sampling_phase','ac_coupling','blocking_allocation','dark_only_nep','browser_budget_reference']]:329 fresh[name]=[]330 for i,r in enumerate(rows):331 if name=='etendue':332 g,a,dd,df=etendue_bounds(r['aperture_diameter_mm']*1e-3,math.radians(r['full_circular_field_deg']/2),math.radians(r['internal_filter_half_angle_deg']))333 v=dict(projected_etendue_m2_sr=g,minimum_detector_area_mm2=a*1e6,minimum_detector_diameter_mm=dd*1e3,minimum_filter_diameter_mm=df*1e3)334 elif name=='parametric_filter':335 ss=sigma_from_fwhm(r['assumed_source_Gaussian_FWHM_nm']*1e-9);sf=sigma_from_fwhm(r['filter_Gaussian_FWHM_nm']*1e-9)336 center=shifted_center(337.1e-9,math.radians(r['incidence_deg']),r['assumed_effective_index']);shift=center-337.1e-9337 cap=gaussian_overlap(ss,sf,shift);band=gaussian_equivalent_band(sf)338 ref_sf=sigma_from_fwhm(4e-9);ref=gaussian_overlap(ss,ref_sf,0)/math.sqrt(gaussian_equivalent_band(ref_sf))339 v=dict(filter_center_nm=center*1e9,shift_nm=shift*1e9,equivalent_rectangular_band_nm=band*1e9,unit_peak_capture_fraction=cap,conditional_background_dominated_SNR_relative_to_centered4nm=cap/math.sqrt(band)/ref)340 if r['incidence_deg'] in [0,5,10,20]:341 integral=gaussian_overlap_integral(ss,sf,shift);t.add('Gaussian_quadrature',i,'capture',integral,cap,rtol=1e-10,atol=2e-13);quad.append({'row':i,'formula':cap,'independent_integral':integral,'error':abs(integral-cap)})342 elif name=='sampling_phase':v={'noiseless_reconstructed_10_90_us':sampled_rise(r['true_10_90_us']*1e-6,r['sample_rate_s'],r['phase_fraction'])*1e6}343 elif name=='ac_coupling':344 tau,fc=hp_droop_timeconstant(r['plateau_duration_s'],r['maximum_fractional_droop']);v=dict(minimum_highpass_timeconstant_s=tau,maximum_highpass_corner_Hz=fc)345 elif name=='blocking_allocation':v={'required_uniform_blocking_OD':blocking_od(r['assumed_unblocked_detector_weighted_out_to_in_background_ratio'],r['allowed_leakage_fraction_of_inband_background'])}346 elif name=='dark_only_nep':v={'dark_only_entrance_referred_NEP_W_sqrtHz':dark_nep(r['wavelength_nm']*1e-9,r['assumed_end_to_end_efficiency'],r['assumed_dark_primary_cps'],r['assumed_excess_noise_factor'])}347 else:348 area=aperture_area(r['aperture_diameter_mm']*1e-3)349 rate=continuum_rate(r['assumed_continuum_W_m2_sr_nm'],r['equivalent_band_nm'],area,math.radians(r['full_circular_FOV_deg']/2),r['efficiency'],r['wavelength_nm']*1e-9)350 v=dict(sky_primary_cps=rate,expected_background_count=(rate+r['dark_primary_cps'])*r['window_s'],registered_signal_count_for_input_fluence=counts_from_received_fluence(r['input_aperture_fluence_ph_m2'],area,r['efficiency']))351 for k,x in v.items():352 atol=2e-7 if name=='sampling_phase' else 0.353 if name=='parametric_filter':atol=(1e-12 if k in ['filter_center_nm','shift_nm'] else 2e-14)354 t.add(name,i,k,x,r[k],atol=atol)355 fresh[name].append({'row':i,**v})356 result=t.result();result['rows_checked']={k:len(v) for k,v in fresh.items()};result['independent_Gaussian_quadratures']=len(quad);result['Gaussian_quadratures']=quad357 assert len(quad)==48358 limits=unit_and_limiting_checks();result['unit_and_limiting_checks']=limits359 if not all(limits.values()):result['status']='FAIL'360 if write:save('design',result,fresh)361 return result,fresh362363def check_diagnostics(write=True):364 d=load('analytic_results.json');t=Checks();fresh={};atm=ReferenceAtmosphere()365 sr,tr0=rayleigh_bodhaine(337.1e-9);sc,n0,scale=rayleigh_frozen(337.1e-9)366 fresh['rayleigh']=[]; rayleigh_grid_audit=[]367 for i,r in enumerate(d['rayleigh']['rows']):368 z=r['geometric_altitude_m'];p,T=atm.state(z);column,tail=atm.pressure_column(z);tau=column*sr;frozen=sc*n0*scale*math.exp(-z/scale)369 v=dict(geopotential_altitude_m=atm.geopotential(z),pressure_Pa=p,temperature_K=T,number_column_m2=column,tau_R_standard_density_integral=tau,T_R_standard_direct=math.exp(-tau),tau_R_pressure_ratio_simple=tr0*p/101325.,tau_R_frozen_CloudScat_exp7200=frozen,T_R_frozen_CloudScat_direct=math.exp(-frozen),standard_minus_frozen_direct_percentage_points=100*(math.exp(-tau)-math.exp(-frozen)))370 # The source table is a 1 m density trapezoid, not the exact hydrostatic371 # integral. A relative tolerance on a subtraction of two transmissions372 # is ill-conditioned. Use its recorded 1 m/2 m refinement envelope,373 # propagated to percentage points; this is a diagnostic tolerance,374 # not a rigorous bound on atmospheric or discretization uncertainty.375 refinement=abs(r['number_column_1m_2m_relative_difference'])376 pp_tolerance=100*max(tau,r['tau_R_standard_density_integral'])*refinement+5e-12377 col_high,_=atm.pressure_column(z,192)378 t.add('rayleigh_quadrature',i,'96_vs_192_relative',abs(col_high-column)/column,0,rtol=0,atol=2e-13)379 for k,x in v.items():380 t.add('rayleigh',i,k,x,r[k],rtol=1e-8,381 atol=(pp_tolerance if k=='standard_minus_frozen_direct_percentage_points' else 1e-9))382 rayleigh_grid_audit.append({'altitude_m':z,'original_1m_2m_column_relative_difference':refinement,'independent_96_vs_192_column_relative_difference':abs(col_high-column)/column,'propagated_comparison_tolerance_percentage_points':pp_tolerance,'actual_difference_percentage_points':v['standard_minus_frozen_direct_percentage_points']-r['standard_minus_frozen_direct_percentage_points']})383 fresh['rayleigh'].append({'row':i,**v,'number_tail_m2':tail})384 fresh['reflection']=[];cache={}385 for i,r in enumerate(d['reflection']['rows']):386 h=r['source_above_plane_m'];H0=(r['satellite_altitude_km']-r['plane_altitude_km'])*1e3;half=math.radians(r['full_circular_FOV_deg']/2);rho=r['assumed_Lambert_albedo'];key=(H0,h,half,rho)387 if key not in cache:cache[key]=LambertPlane(H0,h,rho,half)388 a=cache[key];B=a.B389 v=dict(nadir_fluence_ratio_to_vacuum=1+B,reflected_to_direct_ratio=B,independent_radial_integral=B,downward_source_fraction_intercepted_within_view=.5*(1-h/math.hypot(h,a.radius)) if h else .5,all_directions_total_upward_escape_fraction=(1+rho)/2,ideal_25mm_aperture_fraction_of_emitted_photons=direct_fraction(aperture_area(.025),H0-h)*(1+B),minimum_reflected_extra_delay_us=2*h/C*1e6,maximum_reflected_extra_delay_in_FOV_us=float(a.delay(a.qmax))*1e6 if h else 0.,observed_10_90_of_final_plateau_rise_us=a.rise(r['assumed_intrinsic_10_90_linear_rise_us']*1e-6)*1e6)390 for k,x in v.items():t.add('reflection',i,k,x,r[k],rtol=2e-8,atol=(5e-6 if k.endswith('_us') else 0.))391 fresh['reflection'].append({'row':i,**v})392 coefficient=(diffusion_rise(1.,1.)[1]-diffusion_rise(1.,1.)[0]);fresh['diffusion']=[]393 for i,r in enumerate(d['diffusion']['rows']):394 D=r['assumed_D_m2_s'];gamma=r['assumed_uniform_absorption_frequency_s'];L=math.sqrt(D*r['input_nonabsorbing_10_90_rise_us']*1e-6/coefficient)395 a,b,peak=diffusion_rise(L,D,gamma);total=math.exp(-L*math.sqrt(gamma/D));captured=diffusion_exit_integral(L,D,gamma,.003);numeric=diffusion_exit_integral(L,D,gamma)396 v=dict(inferred_depth_m=L,mode_onset_to_peak_us=peak*1e6,leading_10_90_us=(b-a)*1e6,earliest_physical_straight_path_us=L/C*1e6,diffusion_t10_us=a*1e6,equivalent_transport_mean_free_path_m=3*D/C,equivalent_transport_depth=L/(3*D/C),all_time_diffusion_exit_probability=total,within3ms_exit_probability=captured,fraction_of_eventually_exiting_weight_within3ms=captured/total)397 for k,x in v.items():t.add('diffusion',i,k,x,r[k],rtol=2e-9,atol=2e-8)398 t.add('diffusion_integral',i,'all_time_integral',numeric,total,rtol=1e-10,atol=2e-11)399 t.add('diffusion_integral',i,'finite_time_quadrature_vs_CDF',captured,diffusion_cdf(.003,L,D,gamma),rtol=1e-10,atol=2e-11)400 t.add('diffusion',i,'causal_flag',int(a<L/C),int(r['diffusion_t10_precedes_light_flight']),rtol=0,atol=0)401 fresh['diffusion'].append({'row':i,**v,'all_time_independent_integral':numeric})402 fresh['deadtime']=[]403 for i,r in enumerate(d['deadtime']['rows']):404 b=r['background_true_cps'];s=r['peak_true_signal_cps'];tau=r['assumed_nonparalyzable_deadtime_s'];m=deadtime_mean(b+s,tau);lo,hi,ratio=deadtime_fractions(b,s,tau)405 v=dict(peak_total_true_cps=b+s,peak_mean_registered_cps=m,baseline_mean_registered_cps=deadtime_mean(b,tau),peak_total_fraction_lost=1-m/(b+s),effective_signal_occupancy=s*tau/(1+b*tau),output_uncorrected_10_90_rise_us=r['input_linear_10_90_rise_us']*ratio,output_over_input_rise=ratio,source_fraction_at_output10pct=lo,source_fraction_at_output90pct=hi,inverse_mean_rate_slope_at_peak=1/(1-m*tau)**2)406 for k,x in v.items():t.add('deadtime',i,k,x,r[k],rtol=2e-10,atol=2e-10)407 fresh['deadtime'].append({'row':i,**v})408 fresh['fov_population']=[];weights={}409 for i,r in enumerate(d['fov_population']['rows']):410 half=math.radians(r['full_circular_FOV_deg']/2);area=aperture_area(r['diameter_mm']*1e-3);rs=6371000.+r['orbit_altitude_km']*1e3;rc=6371000.+r['cloud_shell_altitude_km']*1e3;eta=r['end_to_end_efficiency'];s=r['required_source_count_for_95pct_detection'];alpha=r['assumed_cumulative_source_number_slope_alpha'];psi,cap=spherical_field(half,rs,rc);dist,cv=shell_geometry(psi,rs,rc)411 axis=s*4*math.pi*(rs-rc)**2/(eta*area);edge=s*4*math.pi*float(dist)**2/(eta*area*float(cv));rate=continuum_rate(r['assumed_received_background_W_m2_sr_nm'],4,area,half,eta,337.1e-9);weight=population_weight(s,area,eta,half,alpha,rs,rc)412 v=dict(projected_etendue_m2_sr=area*circular_angles(half)[1],background_photoelectron_cps=rate,footprint_km2=cap/1e6,axis_required_isotropic_source_photons_in_window=axis,edge_required_isotropic_source_photons_in_window=edge,edge_over_axis_source_threshold=edge/axis,arbitrary_pivot_1e20ph_weighted_exposure_area_m2=weight,hypothetical_17ns_counter_background_occupancy=rate*17e-9)413 for k,x in v.items():t.add('fov_population',i,k,x,r[k],rtol=1e-8,atol=0.)414 key=(r['assumed_received_background_W_m2_sr_nm'],r['diameter_mm'],r['full_circular_FOV_deg'],alpha);weights[key]=weight;fresh['fov_population'].append({'row':i,**v})415 for i,r in enumerate(d['fov_population']['rows']):416 L=r['assumed_received_background_W_m2_sr_nm'];diam=r['diameter_mm'];alpha=r['assumed_cumulative_source_number_slope_alpha'];own=fresh['fov_population'][i]['arbitrary_pivot_1e20ph_weighted_exposure_area_m2']417 for col,denom in [('relative_population_weight_to_same_aperture60deg',weights[L,diam,60.,alpha]),('relative_population_weight_to_25mm60deg',weights[L,25.,60.,alpha])]:418 x=own/denom;t.add('fov_population',i,col,x,r[col],rtol=1e-8);fresh['fov_population'][i][col]=x419 fresh['auroral_conversion']=[]420 for i,r in enumerate(d['auroral_counts']['rows']):421 area=aperture_area(r['diameter_mm']*1e-3);g=area*circular_angles(math.radians(r['full_circular_FOV_deg']/2))[1];rate=r['assumed_received_band_equivalent_R']*RAYLEIGH_PHOTONS*g*r['eta'];B=(rate+50)*r['observation_window_s'];S=r['signal_mean_for_95pct_power'];inc=S/r['eta']422 v=dict(projected_etendue_m2_sr=g,nominal_four_nm_mean_W_m2_sr_nm=r['assumed_received_band_equivalent_R']*RAYLEIGH_PHOTONS*photon_energy(337.1e-9)/4,incident_linear_photoelectron_rate_s=rate,unmultiplied_mean_photocurrent_A=rate*Q_E,known_background_mean_count=B,Poisson_sigma_count=math.sqrt(B),aperture_photon_equivalent_rms=math.sqrt(B)/r['eta'],aperture_photons_for_95pct_power=inc,aperture_fluence_ph_m2_for_95pct_power=inc/area,rectangular_source_irradiance_W_m2_for_95pct_power=inc*photon_energy(337.1e-9)/(area*r['observation_window_s']),hypothetical_17ns_deadtime_occupancy=rate*17e-9,hypothetical_nonparalyzable_fraction_lost=rate*17e-9/(1+rate*17e-9))423 for k,x in v.items():t.add('auroral_conversion',i,k,x,r[k],rtol=2e-10)424 fresh['auroral_conversion'].append({'row':i,**v})425 result=t.result();result['rows_checked']={k:len(v) for k,v in fresh.items()};result['diffusion_rise_coefficient']=coefficient426 result['rayleigh_grid_tolerance_audit']=rayleigh_grid_audit427 result['scope_exclusions']=['Poisson threshold n and required source S are inputs in FOV/auroral conversion checks, not revalidated by this module.','Noll/Bianchi original raw spectrum/sky cube are not present in this curated bundle: original source spectral folds are not independently re-integrated.','No scattering Monte Carlo, physical atmospheric validation, actual mission event rate, source current, radiation lifetime or source spectrum optimization is asserted.']428 if write:save('diagnostics',result,fresh)429 return result,fresh430431def dimensional_checks():432 """Declared SI exponent checks (mass,length,time,current,temperature,amount).433434 This checks the equations' dimensions; it is not an automatic unit parser for435 arbitrary caller inputs. Angular units and photon/count labels remain semantic.436 """437 mass,length,seconds,current,kelvin,amount=np.eye(6)438 zero=np.zeros(6);energy=mass+2*length-2*seconds;power=energy-seconds439 hdim=energy+seconds;cdim=length-seconds440 pressure=mass-length-2*seconds;gas=energy-kelvin-amount441 gravity=length-2*seconds;molmass=mass-amount;boltzmann=energy-kelvin442 equal=lambda a,b:bool(np.allclose(a,b,atol=0,rtol=0))443 return {444 'dimension_Egamma_J':equal(hdim+cdim-length,energy),445 'dimension_radiance_integral_rate_per_s':equal((power-3*length)+2*length+length-energy,-seconds),446 'dimension_fluence_area_count':equal(-2*length+2*length,zero),447 'dimension_Rayleigh_tau':equal(2*length-2*length,zero),448 'dimension_pressure_column_m-2':equal(gas-gravity-molmass-boltzmann+pressure,-2*length),449 'dimension_NEP_W_per_sqrtHz':equal(energy-.5*seconds,power+.5*seconds),450 'dimension_NEP_squared_ENBW_power_variance':equal(2*(energy-.5*seconds)-seconds,2*power),451 'dimension_diffusion_first_exit_per_s':equal(length-.5*(2*length-seconds)-1.5*seconds,-seconds),452 'dimension_killing_exponent':equal(length+.5*(-seconds-(2*length-seconds)),zero),453 'dimension_deadtime_occupancy':equal(-seconds+seconds,zero)}454455def unit_and_limiting_checks():456 """Small independently interpretable hand limits, not mirrors of every table row."""457 o,op=circular_angles(math.pi/2);sq,sqp=square_angles(math.radians(40))458 ss=sigma_from_fwhm(1e-9)459 tau,fc=hp_droop_timeconstant(.003,.01)460 a=aperture_area(.025);nep=dark_nep(337.1e-9,.2,50.,1.)461 guard=Checks();guard.add('negative_finite_guard',0,'control',math.nan,0.)462 return {**dimensional_checks(),'nonfinite_comparison_rejected':guard.result()['status']=='FAIL','SI_photon_J':abs(photon_energy(337.1e-9)-5.892749502073357e-19)<1e-32,463 'hemisphere_solid_2pi_projected_pi':abs(o-2*math.pi)<1e-14 and abs(op-math.pi)<1e-14,464 'square80diag_smaller_than_circle80':sqp<circular_angles(math.radians(40))[1] and sq>sqp,465 'isotropic_aperture_dimensionless':0<direct_fraction(a,534000.)<1,466 'fluence_referred_count':abs(counts_from_received_fluence(1e8,a,.2)-9817.477042468105)<1e-8,467 'Gaussian_identical_centered_1oversqrt2':abs(gaussian_overlap(ss,ss,0)-1/math.sqrt(2))<1e-14,468 'Gaussian_4nm_equivalent_4p2579nm':abs(gaussian_equivalent_band(sigma_from_fwhm(4e-9))*1e9-4.257868077724905)<1e-12,469 'ENBW_onepole_10kHz':abs(electrical_enbw_first_order(1e4)-15707.963267948966)<1e-10,470 'NEP_boxcar_variance':abs((nep*math.sqrt(boxcar_enbw(.003)))**2-(photon_energy(337.1e-9)/.2)**2*50/.003)<1e-45,471 'AC_3ms_1pct_tau_0p2985s':abs(tau-.2984974874202665)<1e-12,472 'deadtime_zero_gives_unit_rise':abs(deadtime_fractions(50,1e6,0)[2]-1)<1e-12,473 'Lambert_zero_height_gain3_escape1':LambertPlane(534000,0,1,math.pi/6).B==2,474 'diffusion_degenerate_family':np.allclose(diffusion_kernel(np.array([1e-5,1e-4,1e-3]),1000.,3e9),diffusion_kernel(np.array([1e-5,1e-4,1e-3]),2000.,1.2e10),rtol=1e-14,atol=1e-20)}475476def save(prefix,result,recomputed):477 out=ROOT/'checks';out.mkdir(exist_ok=True)478 result['created_utc']=datetime.now(timezone.utc).isoformat();result['module_sha256']=hashlib.sha256(Path(__file__).read_bytes()).hexdigest()479 result['inputs']=[{'path':'data/'+name,'sha256':hashlib.sha256((ROOT/'data'/name).read_bytes()).hexdigest()} for name in ['design_results.json','analytic_results.json']]480 (out/(prefix+'_checks.json')).write_text(json.dumps(result,indent=2,allow_nan=False),encoding='utf8')481 (out/(prefix+'_recomputed.json')).write_text(json.dumps(recomputed,indent=2,allow_nan=False),encoding='utf8')482483if __name__=='__main__':484 start=time.monotonic();a,_=check_design();b,_=check_diagnostics()485 print(json.dumps({'design':a['status'],'diagnostics':b['status'],'design_comparisons':a['comparisons'],'diagnostics_comparisons':b['comparisons'],'design_failures':a['failures'][:8],'diagnostics_failures':b['failures'][:8],'wall_seconds':time.monotonic()-start},indent=2))486 if a['status']!='PASS' or b['status']!='PASS':raise SystemExit(1)487