"""Independent SI design algebra and ideal-model diagnostics.

This module does not import the original calculation scripts. Expected tables are
loaded only as comparison targets and as declared input-parameter grids. The
count-threshold S/n columns in the FOV and auroral tables are deliberately inputs:
their exact Poisson validation belongs to the statistics verification module.
Run from any directory: python modules/design_review.py
"""
from pathlib import Path
from datetime import datetime, timezone
from functools import lru_cache
import ast, hashlib, json, math, platform, time
import numpy as np

ROOT = Path(__file__).resolve().parents[1]
H = 6.62607015e-34
C = 299792458.0
Q_E = 1.602176634e-19
K_B = 1.380649e-23
RAYLEIGH_PHOTONS = 1e10 / (4 * math.pi)

def load(name):
    return json.loads((ROOT / 'data' / name).read_text(encoding='utf-8'))

@lru_cache(None)
def nodes(order):
    return np.polynomial.legendre.leggauss(order)

def integrate(fun, lo, hi, order=192):
    """Fixed Gauss-Legendre integral, including an explicit interval Jacobian."""
    if hi <= lo:
        return 0.0
    x, w = nodes(order)
    a = (hi - lo) / 2
    return float(a * np.dot(w, fun((hi + lo) / 2 + a * x)))

def bisect(fun, lo, hi, iterations=72):
    """Bracketed root; no original root-solver implementation is reused."""
    flo, fhi = fun(lo), fun(hi)
    if flo == 0: return lo
    if fhi == 0: return hi
    if flo * fhi > 0: raise ValueError(('unbracketed', lo, hi, flo, fhi))
    for _ in range(iterations):
        mid = (lo + hi) / 2
        fm = fun(mid)
        if fm == 0: return mid
        if flo * fm <= 0: hi = mid
        else: lo, flo = mid, fm
    return (lo + hi) / 2

def photon_energy(wavelength_m):
    """D01: energy in J for a wavelength explicitly supplied in metres."""
    if wavelength_m <= 0: raise ValueError('wavelength must be positive')
    return H * C / wavelength_m

def circular_angles(half_angle_rad):
    """D02: ordinary solid angle and projected angular integral, in sr."""
    return (4 * math.pi * math.sin(half_angle_rad / 2)**2,
            math.pi * math.sin(half_angle_rad)**2)

def square_angles(diagonal_half_angle_rad):
    """D03: square field; its diagonal half-angle is not a circular half-angle."""
    u = math.tan(diagonal_half_angle_rad) / math.sqrt(2)
    omega = 4 * math.atan(u*u / math.sqrt(1 + 2*u*u))
    x, w = nodes(96); xx = u*x
    projected = float(u*u*np.einsum('i,j,ij->', w, w,
                        1 / (1 + xx[:, None]**2 + xx[None, :]**2)**2))
    return omega, projected

def aperture_area(diameter_m):
    return math.pi * diameter_m**2 / 4

def counts_from_received_fluence(fluence_ph_m2, area_m2, efficiency, incidence_rad=0.):
    """D04: receiver-plane photon fluence -> detected expected counts."""
    return fluence_ph_m2 * area_m2 * efficiency * math.cos(incidence_rad)

def continuum_rate(L_W_m2_sr_per_nm, band_nm, area_m2, half_angle_rad,
                   efficiency, wavelength_m):
    """D05: uniform energy radiance, rectangular/equivalent band, flat response."""
    return (L_W_m2_sr_per_nm * band_nm * area_m2 *
            circular_angles(half_angle_rad)[1] * efficiency / photon_energy(wavelength_m))

def etendue_bounds(diameter_m, field_half_rad, filter_half_rad):
    """D06: necessary air-to-air phase-space bounds; no real optical prescription."""
    a = aperture_area(diameter_m)
    g = a * circular_angles(field_half_rad)[1]
    detector = g / math.pi  # detector acceptance is an entire air hemisphere
    filt = g / (math.pi * math.sin(filter_half_rad)**2)
    return g, detector, math.sqrt(4*detector/math.pi), math.sqrt(4*filt/math.pi)

def shifted_center(wavelength_m, incidence_rad, effective_index):
    """D07: scalar effective-index coating approximation, not measured coating data."""
    return wavelength_m * math.sqrt(1 - (math.sin(incidence_rad)/effective_index)**2)

def sigma_from_fwhm(fwhm_m):
    return fwhm_m / math.sqrt(8 * math.log(2))

def gaussian_overlap(source_sigma_m, filter_sigma_m, shift_m):
    """D08: normalized photon-source Gaussian times a unit-peak filter Gaussian."""
    variance = source_sigma_m**2 + filter_sigma_m**2
    return filter_sigma_m / math.sqrt(variance) * math.exp(-shift_m**2/(2*variance))

def gaussian_overlap_integral(source_sigma_m, filter_sigma_m, shift_m):
    """D08 independent quadrature in standard-normal source coordinates."""
    return integrate(lambda z: np.exp(-z*z/2)/math.sqrt(2*math.pi) *
              np.exp(-(source_sigma_m*z-shift_m)**2/(2*filter_sigma_m**2)), -12., 12., 256)

def gaussian_equivalent_band(filter_sigma_m):
    """D09: spectral equivalent rectangular width in m, distinct from electrical ENBW."""
    return math.sqrt(2 * math.pi) * filter_sigma_m

def electrical_enbw_first_order(corner_Hz):
    """D10: one-sided ENBW of H(f)=1/(1+i*f/fc)."""
    return math.pi * corner_Hz / 2

def boxcar_enbw(exposure_s):
    """D10: one-sided ENBW for a unit-DC-gain rectangular average."""
    return 1 / (2 * exposure_s)

def blocking_od(unblocked_out_to_in_ratio, allowed_leak_fraction):
    """D11: uniform blocking bound on a detector-weighted photon-rate ratio."""
    return max(0., math.log10(unblocked_out_to_in_ratio/allowed_leak_fraction))

def dark_nep(wavelength_m, efficiency, primary_dark_s, excess_noise_factor=1.):
    """D12: entrance-referred one-sided amplitude spectral density, W/sqrt(Hz)."""
    return photon_energy(wavelength_m)/efficiency * math.sqrt(2*excess_noise_factor*primary_dark_s)

def hp_droop_timeconstant(plateau_s, maximum_droop):
    """D13: electrical step y=e^(-t/tau); return minimum tau and maximum fc."""
    tau = plateau_s / (-math.log1p(-maximum_droop))
    return tau, 1/(2*math.pi*tau)

def ramp_bin_average(left_s, right_s, onset_s, ramp_s):
    """D14: integrate ramp by separating its linear piece and plateau rectangle."""
    a, b = left_s-onset_s, right_s-onset_s
    lower, upper = max(a, 0.), min(b, ramp_s)
    triangle = (upper-lower)*(upper+lower)/(2*ramp_s) if upper>lower else 0.
    rectangle = max(0., b-max(a, ramp_s))
    return (triangle + rectangle)/(right_s-left_s)

def sampled_rise(rise_s, sample_rate_s, phase):
    """D14: bin-average linear ramp, then interpolate crossings at exposure midpoints."""
    dt = 1/sample_rate_s; ramp = rise_s/.8; onset = phase*dt
    edges = np.arange(-2, math.ceil((onset+ramp)/dt)+5, dtype=float)*dt
    y = np.array([ramp_bin_average(a,b,onset,ramp) for a,b in zip(edges[:-1],edges[1:])])
    x = (edges[:-1]+edges[1:])/2
    def crossing(level):
        j = int(np.nonzero(y >= level)[0][0])
        return x[j-1] + (x[j]-x[j-1])*(level-y[j-1])/(y[j]-y[j-1])
    return crossing(.9)-crossing(.1)

class ReferenceAtmosphere:
    """P01: documented USSA1976 layer constants; no new polar atmosphere is invented."""
    R = 6356766.; g = 9.80665; gas = 8314.32; M = 28.9644
    levels = np.array([0.,11000.,20000.,32000.,47000.,51000.,71000.,84852.])
    lapse = np.array([-.0065,0.,.001,.0028,0.,-.0028,-.002])
    def __init__(self):
        self.T=[288.15]; self.p=[101325.]
        for j,L in enumerate(self.lapse):
            dh=self.levels[j+1]-self.levels[j]; t=self.T[-1]
            self.T.append(t+L*dh)
            log_ratio=(-self.g*self.M*dh/(self.gas*t) if L==0 else
                       -self.g*self.M/(self.gas*L)*math.log1p(L*dh/t))
            self.p.append(self.p[-1]*math.exp(log_ratio))
    def geopotential(self, z_m): return self.R*z_m/(self.R+z_m)
    def state_H(self, H_m):
        i=min(6,max(0,int(np.searchsorted(self.levels,H_m,side='right')-1)))
        L=self.lapse[i];delta=H_m-self.levels[i]; T=self.T[i]+L*delta
        logp=(-self.g*self.M*delta/(self.gas*self.T[i]) if L==0 else
              -self.g*self.M/(self.gas*L)*math.log1p(L*delta/self.T[i]))
        return self.p[i]*math.exp(logp),T
    def state(self,z_m): return self.state_H(self.geopotential(z_m))
    def pressure_column(self,z_m,order=96):
        """Independent pressure-coordinate quadrature; original used density trapezoids."""
        H0=self.geopotential(z_m); total=0.
        for i,L in enumerate(self.lapse):
            lo=max(H0,self.levels[i]); hi=self.levels[i+1]
            if lo>=hi: continue
            pl,_=self.state_H(lo); ph,_=self.state_H(hi)
            def jac(p):
                if L==0:
                    hh=self.levels[i]-self.gas*self.T[i]/(self.g*self.M)*np.log(p/self.p[i])
                else:
                    hh=self.levels[i]+self.T[i]/L*np.expm1(-self.gas*L/(self.g*self.M)*np.log(p/self.p[i]))
                return 1/(1-hh/self.R)**2
            total+=integrate(jac,ph,pl,order)
        # Retain and expose the SAME small pressure-tail approximation at ~86 km.
        tail=self.p[-1]/(1-self.levels[-1]/self.R)**2
        constant=self.gas/(self.g*self.M*K_B)
        return constant*(total+tail), constant*tail

def rayleigh_bodhaine(wavelength_m):
    """P02: Bodhaine1999 Eq29 fit; wavelength argument to rational fit is micrometres."""
    u=wavelength_m*1e6
    r=(1.0455996-341.29061/u**2-.90230850*u*u)/(1+.0027059889/u**2-85.968563*u*u)
    return r*1e-32, .0021520*r

def rayleigh_frozen(wavelength_m):
    """P02: separate frozen CloudScat source arithmetic, not a replacement atmosphere."""
    u=wavelength_m*1e6; n0=2.6867811e25
    n=1+1e-8*(8060.51+2480990/(132.274-1/u**2)+17455.7/(39.32957-1/u**2))
    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)/1e6
    sigma=24*math.pi**3*((n*n-1)/(n*n+2))**2*king/(wavelength_m**4*(n0*273.15/288.15)**2)
    return sigma,n0,7200.

def direct_fraction(area_m2, distance_m, ray_tau=0., cloud_tau=0.):
    """P03: ballistic photons per emitted isotropic photon, not total cloudy light."""
    return area_m2/(4*math.pi*distance_m**2)*math.exp(-ray_tau-cloud_tau)

class LambertPlane:
    """P04-P05: ideal instantaneous Lambert surface, expressed in log radial coordinate."""
    def __init__(self, observer_height_m, source_height_m, albedo, field_half_rad):
        self.H=observer_height_m;self.h=source_height_m;self.rho=albedo
        if not 0<=source_height_m<observer_height_m:raise ValueError('0 <= h < H required')
        self.radius=observer_height_m*math.tan(field_half_rad)
        self.qmax=math.log1p(self.radius/source_height_m) if source_height_m else 0.
        self.B=self.integral(0,self.qmax) if source_height_m else 2*albedo
    def density(self,q):
        u=np.expm1(q); a=self.h/self.H
        return 2*self.rho*(1-a)**2*u*np.exp(q)/(np.hypot(1,u)**3*(1+(a*u)**2)**2)
    def delay(self,q):
        u=np.expm1(q); x=self.h*u/self.H
        return (self.h*(np.hypot(1,u)+1)+self.H*x*x/(np.hypot(1,x)+1))/C
    def integral(self,a,b): return integrate(self.density,a,b,192)
    def cumulative_q(self,t):
        if t<=self.delay(0.):return 0.
        if t>=self.delay(self.qmax):return self.qmax
        return bisect(lambda q:float(self.delay(q))-t,0.,self.qmax)
    def received_ramp(self,t,ramp_s):
        if not self.h:return float(np.clip(t/ramp_s,0,1))
        a=self.cumulative_q(t-ramp_s);b=self.cumulative_q(t)
        reflected=self.integral(0,a)+integrate(lambda q:self.density(q)*(t-self.delay(q))/ramp_s,a,b,128)
        return (float(np.clip(t/ramp_s,0,1))+reflected)/(1+self.B)
    def rise(self,input_10_90_s):
        if not self.h or not self.rho:return input_10_90_s
        ramp=input_10_90_s/.8;hi=ramp+float(self.delay(self.qmax))
        lo=bisect(lambda t:self.received_ramp(t,ramp)-.1,0.,hi)
        up=bisect(lambda t:self.received_ramp(t,ramp)-.9,0.,hi)
        return up-lo

def diffusion_kernel(t_s, depth_m, D_m2_s, gamma_s=0.):
    """P06: first-hitting density s^-1; an acausal diffusion approximation, not early RTE."""
    t=np.asarray(t_s,dtype=float);safe=np.maximum(t,np.finfo(float).tiny)
    with np.errstate(over='ignore',under='ignore'):
        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*safe
        result=np.exp(log_density)
    return np.where(t>0,result,0.)

def diffusion_peak(depth_m,D_m2_s,gamma_s=0.):
    a=depth_m**2/(4*D_m2_s)
    return 2*a/(1.5+math.sqrt(2.25+4*gamma_s*a))

def diffusion_rise(depth_m,D_m2_s,gamma_s=0.):
    peak=diffusion_peak(depth_m,D_m2_s,gamma_s); a=depth_m**2/(4*D_m2_s)
    def log_ratio(t):return -1.5*math.log(t/peak)-a*(1/t-1/peak)-gamma_s*(t-peak)
    t10=bisect(lambda t:log_ratio(t)-math.log(.1),peak*1e-9,peak)
    t90=bisect(lambda t:log_ratio(t)-math.log(.9),peak*1e-9,peak)
    return t10,t90,peak

def diffusion_exit_integral(depth_m,D_m2_s,gamma_s=0.,window_s=None):
    """P06 independent log-u quadrature of h(t), avoiding the original scipy.integrate route."""
    a=gamma_s*depth_m**2/(4*D_m2_s)
    lower=-32. if window_s is None else math.log(depth_m/math.sqrt(4*D_m2_s*window_s))
    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)

def diffusion_cdf(window_s,depth_m,D_m2_s,gamma_s=0.):
    """P06 closed finite-time exit probability, evaluated in this bounded grid."""
    a=depth_m/math.sqrt(4*D_m2_s*window_s)
    q=math.sqrt(gamma_s*window_s);b=depth_m*math.sqrt(gamma_s/D_m2_s)
    return .5*(math.exp(-b)*math.erfc(a-q)+math.exp(b)*math.erfc(a+q))

def deadtime_mean(rate_s,deadtime_s):
    """P07: stationary nonparalyzable renewal mean, no Poisson output assertion."""
    return 1/(1/rate_s+deadtime_s) if rate_s else 0.

def deadtime_fractions(background_s,signal_s,deadtime_s):
    """P07 inverse mean-rate route, independent of original closed shape formula."""
    mb=deadtime_mean(background_s,deadtime_s);mp=deadtime_mean(background_s+signal_s,deadtime_s)
    values=[]
    for y in [.1,.9]:
        measured=mb+y*(mp-mb);true=measured/(1-measured*deadtime_s)
        values.append((true-background_s)/signal_s)
    return values[0],values[1],(values[1]-values[0])/.8

def spherical_field(half_rad,observer_radius_m,cloud_radius_m):
    """P08: central cap angle and cap area for the near shell intersection."""
    psi=math.asin(observer_radius_m/cloud_radius_m*math.sin(half_rad))-half_rad
    return psi,4*math.pi*cloud_radius_m**2*math.sin(psi/2)**2

def shell_geometry(psi,observer_radius_m,cloud_radius_m):
    """P08: independent shell-central-angle coordinates, not observer-ray quadrature."""
    s2=np.sin(np.asarray(psi)/2)**2;h=observer_radius_m-cloud_radius_m
    distance=np.sqrt(h*h+4*observer_radius_m*cloud_radius_m*s2)
    cos_view=(h+2*cloud_radius_m*s2)/distance
    return distance,cos_view

def population_weight(required_detected_count,area_m2,efficiency,half_rad,alpha,
                      observer_radius_m=6921000.,cloud_radius_m=6388000.,pivot_ph=1e20):
    """P09: hypothetical bright-enough weight, not a full detection integral or rate."""
    end,_=spherical_field(half_rad,observer_radius_m,cloud_radius_m)
    def integrand(psi):
        distance,cv=shell_geometry(psi,observer_radius_m,cloud_radius_m)
        n=required_detected_count*4*math.pi*distance**2/(efficiency*area_m2*cv)
        return (n/pivot_ph)**(-alpha)*2*math.pi*cloud_radius_m**2*np.sin(psi)
    return integrate(integrand,0.,end,192)

class Checks:
    def __init__(self):self.fields={};self.failures=[];self.comparisons=0
    def add(self,group,row,column,fresh,expected,rtol=2e-10,atol=0.):
        key=group+'.'+column
        x=self.fields.setdefault(key,{'comparisons':0,'max_absolute_error':0.,'max_relative_error':0.})
        for label,value in [('relative_tolerance',rtol),('absolute_tolerance_in_column_units',atol)]:
            interval=x.setdefault(label+'_range',[value,value]);interval[0]=min(interval[0],value);interval[1]=max(interval[1],value)
        x['comparisons']+=1;self.comparisons+=1
        if not math.isfinite(float(fresh)) or not math.isfinite(float(expected)):
            self.failures.append({'group':group,'row':row,'column':column,'fresh':repr(fresh),'expected':repr(expected),'reason':'nonfinite value rejected'})
            return
        err=abs(float(fresh)-float(expected));x['max_absolute_error']=max(x['max_absolute_error'],err)
        if float(expected)!=0:x['max_relative_error']=max(x['max_relative_error'],err/abs(float(expected)))
        else:x['zero_reference_comparisons']=x.get('zero_reference_comparisons',0)+1
        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})
    def result(self):return {'status':'PASS' if not self.failures else 'FAIL','comparisons':self.comparisons,'fields':self.fields,'failures':self.failures}

def check_design(write=True):
    d=load('design_results.json');t=Checks();fresh={};quad=[]
    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])
    t.add('constants',0,'one_Rayleigh_ph_m2_sr_s',RAYLEIGH_PHOTONS,d['constants']['one_Rayleigh_ph_m2_sr_s'])
    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']]:
        fresh[name]=[]
        for i,r in enumerate(rows):
            if name=='etendue':
                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']))
                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)
            elif name=='parametric_filter':
                ss=sigma_from_fwhm(r['assumed_source_Gaussian_FWHM_nm']*1e-9);sf=sigma_from_fwhm(r['filter_Gaussian_FWHM_nm']*1e-9)
                center=shifted_center(337.1e-9,math.radians(r['incidence_deg']),r['assumed_effective_index']);shift=center-337.1e-9
                cap=gaussian_overlap(ss,sf,shift);band=gaussian_equivalent_band(sf)
                ref_sf=sigma_from_fwhm(4e-9);ref=gaussian_overlap(ss,ref_sf,0)/math.sqrt(gaussian_equivalent_band(ref_sf))
                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)
                if r['incidence_deg'] in [0,5,10,20]:
                    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)})
            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}
            elif name=='ac_coupling':
                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)
            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'])}
            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'])}
            else:
                area=aperture_area(r['aperture_diameter_mm']*1e-3)
                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)
                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']))
            for k,x in v.items():
                atol=2e-7 if name=='sampling_phase' else 0.
                if name=='parametric_filter':atol=(1e-12 if k in ['filter_center_nm','shift_nm'] else 2e-14)
                t.add(name,i,k,x,r[k],atol=atol)
            fresh[name].append({'row':i,**v})
    result=t.result();result['rows_checked']={k:len(v) for k,v in fresh.items()};result['independent_Gaussian_quadratures']=len(quad);result['Gaussian_quadratures']=quad
    assert len(quad)==48
    limits=unit_and_limiting_checks();result['unit_and_limiting_checks']=limits
    if not all(limits.values()):result['status']='FAIL'
    if write:save('design',result,fresh)
    return result,fresh

def check_diagnostics(write=True):
    d=load('analytic_results.json');t=Checks();fresh={};atm=ReferenceAtmosphere()
    sr,tr0=rayleigh_bodhaine(337.1e-9);sc,n0,scale=rayleigh_frozen(337.1e-9)
    fresh['rayleigh']=[]; rayleigh_grid_audit=[]
    for i,r in enumerate(d['rayleigh']['rows']):
        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)
        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)))
        # The source table is a 1 m density trapezoid, not the exact hydrostatic
        # integral. A relative tolerance on a subtraction of two transmissions
        # is ill-conditioned. Use its recorded 1 m/2 m refinement envelope,
        # propagated to percentage points; this is a diagnostic tolerance,
        # not a rigorous bound on atmospheric or discretization uncertainty.
        refinement=abs(r['number_column_1m_2m_relative_difference'])
        pp_tolerance=100*max(tau,r['tau_R_standard_density_integral'])*refinement+5e-12
        col_high,_=atm.pressure_column(z,192)
        t.add('rayleigh_quadrature',i,'96_vs_192_relative',abs(col_high-column)/column,0,rtol=0,atol=2e-13)
        for k,x in v.items():
            t.add('rayleigh',i,k,x,r[k],rtol=1e-8,
                  atol=(pp_tolerance if k=='standard_minus_frozen_direct_percentage_points' else 1e-9))
        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']})
        fresh['rayleigh'].append({'row':i,**v,'number_tail_m2':tail})
    fresh['reflection']=[];cache={}
    for i,r in enumerate(d['reflection']['rows']):
        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)
        if key not in cache:cache[key]=LambertPlane(H0,h,rho,half)
        a=cache[key];B=a.B
        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)
        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.))
        fresh['reflection'].append({'row':i,**v})
    coefficient=(diffusion_rise(1.,1.)[1]-diffusion_rise(1.,1.)[0]);fresh['diffusion']=[]
    for i,r in enumerate(d['diffusion']['rows']):
        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)
        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)
        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)
        for k,x in v.items():t.add('diffusion',i,k,x,r[k],rtol=2e-9,atol=2e-8)
        t.add('diffusion_integral',i,'all_time_integral',numeric,total,rtol=1e-10,atol=2e-11)
        t.add('diffusion_integral',i,'finite_time_quadrature_vs_CDF',captured,diffusion_cdf(.003,L,D,gamma),rtol=1e-10,atol=2e-11)
        t.add('diffusion',i,'causal_flag',int(a<L/C),int(r['diffusion_t10_precedes_light_flight']),rtol=0,atol=0)
        fresh['diffusion'].append({'row':i,**v,'all_time_independent_integral':numeric})
    fresh['deadtime']=[]
    for i,r in enumerate(d['deadtime']['rows']):
        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)
        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)
        for k,x in v.items():t.add('deadtime',i,k,x,r[k],rtol=2e-10,atol=2e-10)
        fresh['deadtime'].append({'row':i,**v})
    fresh['fov_population']=[];weights={}
    for i,r in enumerate(d['fov_population']['rows']):
        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)
        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)
        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)
        for k,x in v.items():t.add('fov_population',i,k,x,r[k],rtol=1e-8,atol=0.)
        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})
    for i,r in enumerate(d['fov_population']['rows']):
        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']
        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])]:
            x=own/denom;t.add('fov_population',i,col,x,r[col],rtol=1e-8);fresh['fov_population'][i][col]=x
    fresh['auroral_conversion']=[]
    for i,r in enumerate(d['auroral_counts']['rows']):
        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']
        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))
        for k,x in v.items():t.add('auroral_conversion',i,k,x,r[k],rtol=2e-10)
        fresh['auroral_conversion'].append({'row':i,**v})
    result=t.result();result['rows_checked']={k:len(v) for k,v in fresh.items()};result['diffusion_rise_coefficient']=coefficient
    result['rayleigh_grid_tolerance_audit']=rayleigh_grid_audit
    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.']
    if write:save('diagnostics',result,fresh)
    return result,fresh

def dimensional_checks():
    """Declared SI exponent checks (mass,length,time,current,temperature,amount).

    This checks the equations' dimensions; it is not an automatic unit parser for
    arbitrary caller inputs. Angular units and photon/count labels remain semantic.
    """
    mass,length,seconds,current,kelvin,amount=np.eye(6)
    zero=np.zeros(6);energy=mass+2*length-2*seconds;power=energy-seconds
    hdim=energy+seconds;cdim=length-seconds
    pressure=mass-length-2*seconds;gas=energy-kelvin-amount
    gravity=length-2*seconds;molmass=mass-amount;boltzmann=energy-kelvin
    equal=lambda a,b:bool(np.allclose(a,b,atol=0,rtol=0))
    return {
      'dimension_Egamma_J':equal(hdim+cdim-length,energy),
      'dimension_radiance_integral_rate_per_s':equal((power-3*length)+2*length+length-energy,-seconds),
      'dimension_fluence_area_count':equal(-2*length+2*length,zero),
      'dimension_Rayleigh_tau':equal(2*length-2*length,zero),
      'dimension_pressure_column_m-2':equal(gas-gravity-molmass-boltzmann+pressure,-2*length),
      'dimension_NEP_W_per_sqrtHz':equal(energy-.5*seconds,power+.5*seconds),
      'dimension_NEP_squared_ENBW_power_variance':equal(2*(energy-.5*seconds)-seconds,2*power),
      'dimension_diffusion_first_exit_per_s':equal(length-.5*(2*length-seconds)-1.5*seconds,-seconds),
      'dimension_killing_exponent':equal(length+.5*(-seconds-(2*length-seconds)),zero),
      'dimension_deadtime_occupancy':equal(-seconds+seconds,zero)}

def unit_and_limiting_checks():
    """Small independently interpretable hand limits, not mirrors of every table row."""
    o,op=circular_angles(math.pi/2);sq,sqp=square_angles(math.radians(40))
    ss=sigma_from_fwhm(1e-9)
    tau,fc=hp_droop_timeconstant(.003,.01)
    a=aperture_area(.025);nep=dark_nep(337.1e-9,.2,50.,1.)
    guard=Checks();guard.add('negative_finite_guard',0,'control',math.nan,0.)
    return {**dimensional_checks(),'nonfinite_comparison_rejected':guard.result()['status']=='FAIL','SI_photon_J':abs(photon_energy(337.1e-9)-5.892749502073357e-19)<1e-32,
      'hemisphere_solid_2pi_projected_pi':abs(o-2*math.pi)<1e-14 and abs(op-math.pi)<1e-14,
      'square80diag_smaller_than_circle80':sqp<circular_angles(math.radians(40))[1] and sq>sqp,
      'isotropic_aperture_dimensionless':0<direct_fraction(a,534000.)<1,
      'fluence_referred_count':abs(counts_from_received_fluence(1e8,a,.2)-9817.477042468105)<1e-8,
      'Gaussian_identical_centered_1oversqrt2':abs(gaussian_overlap(ss,ss,0)-1/math.sqrt(2))<1e-14,
      'Gaussian_4nm_equivalent_4p2579nm':abs(gaussian_equivalent_band(sigma_from_fwhm(4e-9))*1e9-4.257868077724905)<1e-12,
      'ENBW_onepole_10kHz':abs(electrical_enbw_first_order(1e4)-15707.963267948966)<1e-10,
      'NEP_boxcar_variance':abs((nep*math.sqrt(boxcar_enbw(.003)))**2-(photon_energy(337.1e-9)/.2)**2*50/.003)<1e-45,
      'AC_3ms_1pct_tau_0p2985s':abs(tau-.2984974874202665)<1e-12,
      'deadtime_zero_gives_unit_rise':abs(deadtime_fractions(50,1e6,0)[2]-1)<1e-12,
      'Lambert_zero_height_gain3_escape1':LambertPlane(534000,0,1,math.pi/6).B==2,
      '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)}

def save(prefix,result,recomputed):
    out=ROOT/'checks';out.mkdir(exist_ok=True)
    result['created_utc']=datetime.now(timezone.utc).isoformat();result['module_sha256']=hashlib.sha256(Path(__file__).read_bytes()).hexdigest()
    result['inputs']=[{'path':'data/'+name,'sha256':hashlib.sha256((ROOT/'data'/name).read_bytes()).hexdigest()} for name in ['design_results.json','analytic_results.json']]
    (out/(prefix+'_checks.json')).write_text(json.dumps(result,indent=2,allow_nan=False),encoding='utf8')
    (out/(prefix+'_recomputed.json')).write_text(json.dumps(recomputed,indent=2,allow_nan=False),encoding='utf8')

if __name__=='__main__':
    start=time.monotonic();a,_=check_design();b,_=check_diagnostics()
    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))
    if a['status']!='PASS' or b['status']!='PASS':raise SystemExit(1)