04. Clear-air and ideal-model diagnostics

This notebook makes the analytic models inspectable and independently recalculates their tabulated values. It does not replace cloud radiative transfer or validate source depth from morphology. The models are useful precisely because their assumptions and failure modes can be written down.

Stable identifiers P01-P09 refer to methods/design_and_diagnostics.md. The independent code does not import the reference scripts. Poisson cut values used in the FOV/auroral conversion rows are declared inputs and must be checked in the statistics notebook.

from pathlib import Path
import sys, os, json, math, hashlib
import numpy as np
from io import BytesIO
from IPython.display import display, Image, Markdown
candidates = [Path.cwd(), *Path.cwd().parents]
ROOT = next(p for p in candidates if (p/'data/design_results.json').exists() and (p/'modules/design_review.py').exists())
if (ROOT/'runtime/deps').is_dir(): sys.path.insert(0, str(ROOT/'runtime/deps'))
sys.path.insert(0, str(ROOT))
os.environ.setdefault('MPLCONFIGDIR', str(ROOT/'checks/diagnostics_mplconfig'))
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from modules import design_review as v
from modules.common import table
plt.rcParams.update({'figure.figsize': (8.2, 4.3), 'font.size': 10, 'axes.grid': True, 'grid.alpha': .22})
def show(fig):
    b=BytesIO(); fig.savefig(b, format='png', dpi=125, bbox_inches='tight')
    display(Image(data=b.getvalue())); plt.close(fig)
design = v.load('design_results.json')
analytic = v.load('analytic_results.json')
print('Independent module:', ROOT/'modules/design_review.py')
print('SHA256:', hashlib.sha256((ROOT/'modules/design_review.py').read_bytes()).hexdigest())

Independent module: D:\Photometer_research_20261022\verification_20261022_r1\modules\design_review.py
SHA256: af4c3056a61716d48c5eb8eec3f1abb67be7dee123d78c6bbc7c6f6cc37322ca

P01-P02. A standard atmospheric column, not an invented cloud

The retained dry USSA1976 layers specify temperature and hydrostatic pressure in geopotential altitude $H=Rz/(R+z)$. Rather than repeating the original 1 m density trapezoid, use

N_{col}={R^*\over g_0 M k_B}\int {dp\over[1-H(p)/R]^2}+N_{tail},\qquad \tau_R=\sigma_RN_{col},\quad T_{direct}=e^{-\tau_R}.

The same small pressure-tail approximation above about 86 km is exposed, not replaced by a fabricated thermosphere. The original exponential CloudScat profile is a separate reference with its own sea-level number density/scale height. Neither profile is polar weather, ozone/aerosol absorption or total cloudy directional transmission. The molecular fit retains the original nominal wavelength coordinate.

atm=v.ReferenceAtmosphere();sigma,tau0=v.rayleigh_bodhaine(337.1e-9);sc,n0,Hscale=v.rayleigh_frozen(337.1e-9)
rows=[]
for z in [0.,15000.,16000.,18000.,18500.]:
 p,T=atm.state(z);col,tail=atm.pressure_column(z,96);hi,_=atm.pressure_column(z,192)
 rows.append({'z_km':z/1000,'p_Pa':p,'T_K':T,'column_m-2':col,'tail_m-2':tail,
   'quadrature96_vs192_relative':abs(hi-col)/col,'standard_tau':sigma*col,'standard_direct':math.exp(-sigma*col),
   'frozen_direct':math.exp(-sc*n0*Hscale*math.exp(-z/Hscale))})
table(rows)
fig,ax=plt.subplots();sel=rows[1:]
ax.plot([r['z_km'] for r in sel],[r['standard_direct'] for r in sel],'o-',label='standard density column')
ax.plot([r['z_km'] for r in sel],[r['frozen_direct'] for r in sel],'s-',label='frozen exponential profile')
ax.set(xlabel='Geometric source altitude (km)',ylabel='Uncollided molecular transmission',title='Reference profiles; no cloud/gas absorption or polar weather');ax.legend();show(fig)

P03. Direct dilution versus diffuse light

For an isotropic point source and small nadir aperture, $f_{direct}=A\exp[-\tau_R-\tau_{cloud}]/(4\pi r^2)$. This is the fraction of all emitted photons crossing the aperture without collision. Detector efficiency is applied separately. Scattering redirects photons; it does not automatically absorb them.

The source/satellite altitude difference sets range. A source at or above the cloud top has zero overhead cloud optical depth, but downward light can still return from the cloud below. Therefore the table is a ballistic reference, not a total cloud-transfer curve.

area=v.aperture_area(.025);rows=[]
for z in [15000.,18000.]:
 col,_=atm.pressure_column(z);tau=sigma*col;distance=550000.-z
 for tc in [0.,.1,.3,1.,3.,10.]:
  fraction=v.direct_fraction(area,distance,tau,tc)
  rows.append({'source_km':z/1000,'overlying_cloud_tau':tc,'range_km':distance/1000,
   'direct_aperture_fraction':fraction,'detected_fraction_assumed_eta0.2':.2*fraction})
table(rows,max_rows=12)

P04-P05. A cloud-below counterexample with an explicit surface model

Use an infinite plane with instantaneous Lambert reflection $\rho/\pi$, an isotropic point source height $h$ above it and an observer height $H$ above it. Integrate ring contributions with both projection factors. The code uses the independent coordinate $q=\ln(1+r/h)$.

Extra delay is $[\sqrt{h^2+r^2}+\sqrt{H^2+r^2}-(H-h)]/c$. Common direct flight time has been subtracted. At $h\to0$, directional gain is $1+2\rho$, but all-direction upward escape is only $(1+\rho)/2$. This creates no photons.

The pulse is an infinite plateau after a linear ramp. Each plotted shape is normalized to its own final plateau. There is no finite-cloud residence time, gas, physical source extent, electronics or noise. The h=0 mathematical limit is not a real embedded cloud-top boundary.

rows=[];fig,ax=plt.subplots();times=np.linspace(0,200e-6,160)
for h in [0.,100.,1000.,3000.]:
 model=v.LambertPlane(535000.,h,1.,math.pi/6)
 rise=model.rise(30e-6)*1e6
 rows.append({'height_above_plane_m':h,'nadir_vacuum_gain':1+model.B,'total_upward_fraction':1.,
              'minimum_extra_delay_us':2*h/v.C*1e6,'output10-90_us':rise})
 ax.plot(times*1e6,[model.received_ramp(t,30e-6/.8) for t in times],label=f'h={h:g}m')
table(rows);ax.set(xlabel='Retarded arrival time (µs)',ylabel='Response / own final plateau',title='Instantaneous Lambert surface, rho 1,30 µs input,60° field');ax.legend();show(fig)

P06. Diffusion: exact mathematics, limited physical scope

For a half-space with absorbing boundary, point injection at depth $L$, diffusion $D$ and uniform killing rate $\gamma$,

h(t)={L\over\sqrt{4\pi D}}t^{-3/2}\exp[-L^2/(4Dt)-\gamma t].

For $\gamma=0$, $t_{peak}=L^2/(6D)$ and the leading amplitude 10-90rise is $0.07432773536 L^2/D$. It is not onset-to-peak, a cumulative arrival quantile, or the rise after source convolution.

Exactly the same $L^2/D$ and gamma produce exactly the same entire kernel. The example family therefore demonstrates an inverse degeneracy, not measured source altitudes. Also, parabolic diffusion has nonzero response before $L/c$; it is not a causal early-time RTE benchmark.

fig,axs=plt.subplots(1,2,figsize=(11,4));tt=np.linspace(1e-6,700e-6,550)
for i,(L,D) in enumerate([(1000.,3e9),(2000.,1.2e10),(3000.,2.7e10)]):
 y=v.diffusion_kernel(tt,L,D)
 axs[0].plot(tt*1e6,y,label=f'L {L/1000:g} km, D {D:.2g}',ls=['-','--',':'][i],lw=2)
axs[0].set(xlabel='Time since ideal impulse (µs)',ylabel='First-exit density (s⁻¹)',title='Equal L²/D: identical entire kernels');axs[0].legend(fontsize=8)
D=7.9e10;coefficient=v.diffusion_rise(1.,1.)[1]-v.diffusion_rise(1.,1.)[0]
L=math.sqrt(D*5.5e-6/coefficient);t10,t90,peak=v.diffusion_rise(L,D)
tt=np.linspace(.05e-6,50e-6,600);y=v.diffusion_kernel(tt,L,D)
axs[1].plot(tt*1e6,y/y.max());axs[1].axvspan(0,L/v.C*1e6,color='tomato',alpha=.2,label='before straight light flight')
axs[1].axvline(t10*1e6,color='black',ls='--',label='diffusion 10% threshold')
axs[1].set(xlabel='Time since ideal impulse (µs)',ylabel='Density / own peak',title='High D early-time approximation can be acausal');axs[1].legend(fontsize=8);show(fig)
print('Example depth:',L,'m;10%time:',t10*1e6,'µs; straight flight:',L/v.C*1e6,'µs')

Example depth: 2417.794804432272 m;10%time: 3.1665804131838793 µs; straight flight: 8.064895363152438 µs

The all-time exit probability is $\exp[-L\sqrt{\gamma/D}]$. The finite-time closed CDF is

F(T)=\tfrac12[e^{-b}\mathrm{erfc}(a-q)+e^b\mathrm{erfc}(a+q)],\quad a=L/\sqrt{4DT},\ q=\sqrt{\gamma T},\ b=L\sqrt{\gamma/D}.

The code also integrates the kernel numerically in log-u coordinates. Gamma is a uniform model parameter, not a measured ice or ozone profile. A 3 ms integration window does not imply a 3 ms intrinsic event or molecular lifetime.

L=1000.;D=3e9;rows=[]
for gamma in [0.,.1,10.,50.]:
 numeric=v.diffusion_exit_integral(L,D,gamma,.003);closed=v.diffusion_cdf(.003,L,D,gamma)
 total=math.exp(-L*math.sqrt(gamma/D))
 rows.append({'gamma_s-1_assumed':gamma,'all_time_exit':total,'within3ms_closed':closed,
    'within3ms_independent_quadrature':numeric,'fraction_of_eventual_exit_in3ms':closed/total})
table(rows)

P07. Raw dead time can compress a rise

The nonparalyzable stationary mean is $m=r/(1+r\tau)$. For a true ramp $b+Sf$, baseline-subtracted registered fraction is $y=f(1+x)/(1+xf)$ with $x=S\tau/(1+b\tau)$. Therefore $T_{out}/T_{in}=(1+x)/(1+x+0.09x^2)$.

The independent implementation instead inverts each registered amplitude using $r=m/(1-m\tau)$. Tau=17ns is a retained study assumption, not a universal detector property. This is a quasi-static raw-mean diagnostic. It does not accuse corrected published MMIA data of the bias, simulate lost-count likelihood, or restore Poisson statistics after inversion.

signals=np.logspace(6,10,120);ratios=[v.deadtime_fractions(50.,s,17e-9)[2] for s in signals]
fig,ax=plt.subplots();ax.semilogx(signals,200*np.array(ratios));ax.axhline(200,color='black',ls='--',label='true 200 µs')
ax.set(xlabel='True peak source rate (s⁻¹)',ylabel='Uncorrected mean10-90rise (µs)',title='Assumed 17 ns nonparalyzable counter, background 50/s');ax.legend();show(fig)
table([{'true_peak_signal_s':s,'registered_peak_total_s':v.deadtime_mean(s+50,17e-9),
        'raw200us_rise_us':200*v.deadtime_fractions(50,s,17e-9)[2]} for s in [1e6,1e8,1e9,1e10]])

P08-P09. A field/population sensitivity, not a useful event rate

Use independent shell-central angle psi. The cloud cap area is $4\pi r_c^2\sin^2(\psi_{max}/2)$, and a point source needs $N_{cut}=4\pi r^2S_{required}/(\eta A\cos\theta)$.

For an assumed cumulative law $C N^{-\alpha}$, integrate $W=\int [N_{cut}/N_{pivot}]^{-\alpha}dA_c$. W is an area-weighted bright-enough population proxy. Its pivot $10^{20}$photons is arbitrary, not a source observation. Source geography, occurrence, exposure and the actual brightness law are unknown.

The supplied ideal-count $S_{required}$ is used as an input; its Poisson threshold is checked in the statistics notebook. The weight omits below-cut detections and is not an integral of the full detection probability. High counter-occupancy scenarios cannot be turned into a real17ns receiver forecast.

area=v.aperture_area(.025);rows=[]
for alpha in [.5,1.,1.5,2.,2.5,3.]:
 pair=[]
 for full in [60.,90.]:
  row=next(r for r in analytic['fov_population']['rows'] if r['assumed_received_background_W_m2_sr_nm']==3e-9 and r['diameter_mm']==25 and r['full_circular_FOV_deg']==full and r['assumed_cumulative_source_number_slope_alpha']==alpha)
  S=row['required_source_count_for_95pct_detection']
  pair.append(v.population_weight(S,area,.2,math.radians(full/2),alpha))
 rows.append({'assumed_alpha':alpha,'W60_m2':pair[0],'W90_m2':pair[1],'conditional90_to60_weight':pair[1]/pair[0]})
table(rows)
fig,ax=plt.subplots();ax.plot([r['assumed_alpha'] for r in rows],[r['conditional90_to60_weight'] for r in rows],'o-')
ax.axhline(1,color='black',ls='--');ax.set(xlabel='Assumed cumulative source-number slope alpha',ylabel='90°/60° bright-enough population weight',title='Conditional comparison; neither observed rate nor full detection integral');show(fig)
table([{'full_circular_deg':f,'footprint_km2':v.spherical_field(math.radians(f/2),6921000.,6388000.)[1]/1e6} for f in [10,30,60,90]])

Full diagnostic comparison and the retained initial tolerance failure

This cell independently checks all declared analytic rows. Original n/S count thresholds are excluded from this module's scope; their conversions are checked. The first attempt used a relative tolerance on a nearly cancelled Rayleigh transmission difference and retained four tiny failures. The final method propagates the original 1 m/2 m column refinement to that derived residual and records an independent 96/192 pressure-coordinate comparison. This is an implementation comparison, not an atmospheric confidence interval.

Separate fault controls must reject wrong aperture fractions, etendue and photocurrent; reproduce them with python checks/diagnostics_negative_controls.py.

The raw Noll spectrum and Bianchi sky cube are not present in this curated bundle, so their actual source spectral folds are explicitly not newly integrated here. No frozen calculation is silently relabeled as a new observation.

result,recomputed=v.check_diagnostics(write=False)
assert result['status']=='PASS',result['failures']
table([{'table':k,'rows_checked':n} for k,n in result['rows_checked'].items()])
print('Status:',result['status'],'scalar comparisons:',result['comparisons'])
print('Independent 10-90 coefficient:',result['diffusion_rise_coefficient'])
table(result['rayleigh_grid_tolerance_audit'])
display(Markdown('**Explicit scope exclusions:**\n\n'+'\n'.join('- '+s for s in result['scope_exclusions'])))

Status: PASS scalar comparisons: 6955
Independent 10-90 coefficient: 0.07432773536040116

<IPython.core.display.Markdown object>

Primary references and applicability

• US Standard Atmosphere1976, Part1 tables2-4, Part4 tables18-19; NASA/TM2005-213659, AppendixA equationsA1-A5, pp34-35: the actual layer constants and hydrostatic model.
• Bodhaine et al.1999, equations29-30: retained dry-air cross-section/reference-pressure fits.
• Lambert reflection model: BRDF normalization. The finite-FOV geometry and delayed convolution are explicit derivations here.
• Bai et al.2023, sections6-7: optical diffusion/source-context and rounded rise relation. The full first-passage derivation and closed CDF are given in the methods companion; no uninspected equation bitmap is claimed read.
• MMIA instrument design, section6.3: counter architecture context. Our17ns nonparalyzable toy is not a universal detector.

All row/field comparisons, input hashes, original-file:function map, tolerance audit and missing physical measurements are exposed in checks/diagnostics_* and methods/design_and_diagnostics.md. Passing these conditional tests does not validate an anvil, a timing-depth classifier, actual detector radiation/noise, source current or mission yield.