02 Transport: equations, native kernels and a worked photon

Goal

Verify the conditional scalar calculation, rather than treating an existing plot as evidence. This light notebook independently recomputes all 56 direct rays and native batch fluences/errors, checks HG and a first-collision integral, and reconstructs a pulse from response mass. The 24 nonzero-GHM cases retain an expected warning: point-nadir rare-path convergence is incomplete. No full RT run, raw ice LUT or new physical cloud validation is performed. See the complete SI derivation.

Setup and assumptions

The model is an isotropic point source inside/at/above a homogeneous finite cylinder: base 14 km,top 17 km,radius 30 km; sources15/16/17/18km; observer 550 km. SSA is assumed 1. Exterior is vacuum; gases,ice absorption,polarization and real cloud structure are absent. The receiver is a small flat aperture with 60 degree full circular FOV. For photon radiance $L$,the scalar RTE is

(c^{-1}\partial_t+\hat\Omega\cdot\nabla)L=-\beta_eL+\beta_s\int p(\hat\Omega'\cdot\hat\Omega)L(\hat\Omega')\,d\Omega'
+\Phi_s(t)\delta^3(\mathbf r-\mathbf r_s)/(4\pi).

$\beta_e,\beta_s$:m⁻¹; $L$:photons m⁻²s⁻¹sr⁻¹; $p$:sr⁻¹; $\Phi_s$:photons s⁻¹. This notebook's times are excess flight delays unless labeled absolute.

from pathlib import Path
import sys,os,json,math
ROOT=next(p for p in [Path.cwd().resolve(),*Path.cwd().resolve().parents] if (p/'modules').is_dir() and (p/'data').is_dir())
for p in [ROOT/'runtime/deps',ROOT/'modules']:
    if p.is_dir():sys.path.insert(0,str(p))
os.environ.setdefault('MPLCONFIGDIR',str(ROOT/'checks/transport_mpl_cache'))
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display,Markdown
from common import table
import transport_review as tr
plt.rcParams.update({'font.size':10,'figure.figsize':(8,4),'axes.grid':True,'grid.alpha':.25})
print('Package:',ROOT,'| NumPy:',np.__version__)
root,data,inputs,input_paths=tr.load_inputs(ROOT)
results=inputs['results'];geometry=results['model']['geometry']
table([{'input':k,'sha256':tr.sha(p)} for k,p in input_paths.items()],max_rows=3)

Package: D:\Photometer_research_20261022\verification_20261022_r1 | NumPy: 2.3.5

1 Direct light: compute from geometry

\beta_e=\tau_{\rm total}/(z_t-z_b),\quad\ell_{\rm above}=\max[0,\min(z_t-z_s,z_t-z_b)],
\quad K_0=\frac{\exp(-\beta_e\ell_{\rm above})}{4\pi R_0^2},\quad R_0=z_o-z_s .

K0 is m⁻² per emitted photon. Multiply by aperture area once. Total directional gain can exceed1; it is not an escape probability. Total tau changes the whole cloud,including below an internal source.

direct_rows=[]
for case in results['cases']:
    z=case['source_altitude_km']*1000
    H=geometry['top_m']-geometry['base_m']
    beta=case['total_vertical_cloud_tau']/H
    R0=geometry['observer_z_m']-z
    tau_above=beta*max(0,min(geometry['top_m']-z,H))
    Kvac=1/(4*np.pi*R0**2)
    gain0=np.exp(-tau_above)
    direct_rows.append(dict(case_id=case['case_id'],tau_above=tau_above,
        direct_gain=float(gain0),direct_m_inverse2=float(Kvac*gain0),
        error=float(gain0-case['direct_relative_to_vacuum'])))
assert max(abs(x['error'])for x in direct_rows)<1e-13
table([x for x in direct_rows if x['case_id'] in ['hg_z15_tau10p0','hg_z16_tau10p0','hg_z17_tau10p0','hg_z18_tau10p0']],max_rows=4)
print('All direct rays checked:',len(direct_rows))

All direct rays checked: 56

2 Native response mass and its uncertainty

For collision j,the next-event summand per unit aperture area is

w_j\omega\,p(\hat\Omega_{\rm in}\cdot\hat u_o)\,
c_A r_j^{-2}\exp(-\beta_e\ell_{j\to o})\,\chi_{\rm FOV}.

The launch weight is physical isotropic density divided by proposal density. The area cosine appears once. Main-grid omega1 needs no additional absorption weight. For excess delay $\Delta_j=(\ell_j+r_j-R_0)/c$,the direct light is an atom at0. NPZ values are scattered bin mass,not density. The last element is overflow,not an ordinary time bin.

cid='hg_z15_tau10p0'
case=next(c for c in results['cases']if c['case_id']==cid)
kernel_path,info_path=tr.native_paths(data,cid)
with np.load(kernel_path,allow_pickle=False)as z:
    masses=z['batch_mass_relative'];native_dt=float(z['dt_s'])
direct0=next(x['direct_gain']for x in direct_rows if x['case_id']==cid)
batch_gain=direct0+masses.sum(axis=1)
mean_gain=float(batch_gain.mean())
SEM=float(batch_gain.std(ddof=1)/np.sqrt(len(batch_gain)))
table([{'batch':i,'native_scattered_mass':float(masses[i].sum()),'total_gain':float(v)}for i,v in enumerate(batch_gain)],max_rows=6)
print('Recomputed gain ± batch SEM:',mean_gain,SEM,'; stored:',case['total_relative_to_vacuum'],case['total_mc_standard_error'])
# Covariance of bin means is needed; never assume native bins are independent.
dev=masses-masses.mean(axis=0);B=len(masses)
variance_sum=float(np.sum(dev.sum(axis=1)**2)/(B*(B-1)))
incorrect_diagonal_variance=float(np.sum(dev**2)/(B*(B-1)))
print('Correct variance of summed mass:',variance_sum)
print('Incorrect sum of bin variances:',incorrect_diagonal_variance)
assert abs(variance_sum-SEM**2)<1e-12

Recomputed gain ± batch SEM: 1.9837274106214 0.008724990325018198 ; stored: 1.9837274106214002 0.008724990325018238
Correct variance of summed mass: 7.612545617166103e-05
Incorrect sum of bin variances: 9.431541538139586e-05

centers_ms=(np.arange(masses.shape[1]-1)+.5)*native_dt*1e3
density=masses[:,:-1].mean(axis=0)/native_dt
density_SE=masses[:,:-1].std(axis=0,ddof=1)/np.sqrt(B)/native_dt
fig,ax=plt.subplots(layout='constrained')
ax.plot(centers_ms,density,label='HG mean scattered response')
ax.fill_between(centers_ms,np.maximum(0,density-density_SE),density+density_SE,alpha=.25,label='±1 batch SEM; bins correlated')
ax.set(xlim=(0,.4),xlabel='Excess arrival delay (ms)',ylabel='Vacuum-normalized response density (s⁻¹)',
       title=cid+' — scattered component only')
ax.text(.98,.95,f'Direct atom at0: mass={direct0:.6g}',ha='right',va='top',transform=ax.transAxes)
ax.legend();plt.show()

3 Reconcile every native case

For independent equal-size batches,

SE(\bar G)^2=\frac{\sum_b(G_b-\bar G)^2}{B(B-1)}
=\sum_{ij}\mathrm{Cov}(\bar m_i,\bar m_j).

The following module reads every native NPZ and batch JSON. Energy/causality metadata are inspectable,but a histogram cannot replay missing individual 3D histories. An analytic vacuum row is not an MCsample.

report,all_rows,inputs=tr.run_review(ROOT,write=True)
table([{'check':k,'value':v}for k,v in report['summary'].items()],max_rows=20)
table([r for r in all_rows if r['case_id']in ['ghm_z15_tau10p0','hg_z15_tau10p0']],
      ['case_id','gain_recomputed','SEM_recomputed','naive_independent_bin_SEM','quality'],max_rows=2)

4 Phase probability: an independent HG check

f(\mu)=2\pi p(\mu),\quad\int_{-1}^{1}f\,d\mu=1,\quad g=\int_{-1}^{1}\mu f\,d\mu,
\qquad f_{\rm HG}=\frac{1-g^2}{2(1+g^2-2g\mu)^{3/2}} .

Equal g does not determine the full phase function. The executable check below integrates HG,checks its CDF and seeded moments. It does not reconstruct GHM without its absent raw LUT.

hg=tr.hg_check(results['model']['phase_input']['derived_g'])
table([{'quantity':k,'value':v}for k,v in hg.items()if k!='scope'],max_rows=12)
phase_receipt=json.loads((ROOT/'data/cloud_phase_QA.json').read_text())
table([{'phase':k,'retained_norm':v['normalization'],'retained_g':v['g_integral'],
        'retained_P2':v['second_Legendre_moment'],'status':'retained QA,not recomputed raw GHM'}
       for k,v in phase_receipt['phases'].items()],max_rows=2)

5 A deterministic first-collision integral

For isotropic scattering and an internal on-axis source,azimuth integrates out:

G^{(1)}=\frac{R_0^2}{2}\int_{-1}^{1}d\mu\int_0^{s_e(\mu)}
 \beta_e e^{-\beta_e[s+\ell_{\rm out}(s,\mu)]}\frac{c_A}{r^2}\,ds .

This independent quadrature splits the side/top/base boundary transitions. It checks the first-collision equation,not a full independent RT solver or the GHM phase.

single=[{'Gauss_order':n,'relative_first_scatter_gain':tr.isotropic_first_scatter_gain(n)}for n in [24,48,96]]
table(single,max_rows=3)
mc=inputs['validation']['single_scatter_reference']['MC']
print('Retained first-collision MC:',mc['mean_first_scatter_gain'],'±',mc['first_scatter_batch_standard_error'])
print('Difference / retained MC SEM:',(mc['mean_first_scatter_gain']-single[-1]['relative_first_scatter_gain'])/mc['first_scatter_batch_standard_error'])
assert tr.isotropic_first_scatter_gain(tau=0)==0

Retained first-collision MC: 0.2981522566870474 ± 0.0009498907194338359
Difference / retained MC SEM: 1.7725941290478089

6 A worked packet you can inspect

This is a new fixed-seed mechanical example,not an original production history. The true source is isotropic; its 50:50 isotropic/HG launch proposal has an exact likelihood weight. The ledger lists free optical depth,path,position,phase probability,observer survival and excess delay. Each NEE value is a raw summand before division by a history count. These summands illustrate the algorithm; they are not used as a reported event-fluence estimate.

ledger=tr.worked_hg_packet(hg['g'])
print('Seed:',ledger['seed'],'| launch direction:',ledger['initial_direction'],
      '| source weight:',ledger['initial_source_weight'],'| end:',ledger['end_status'])
table(ledger['collisions'],['collision','free_optical_depth','free_path_m','z_m','cos_to_observer',
      'p_per_sr','observer_tau','delay_us','NEE_summand_relative_vacuum'],max_rows=10)
display(Markdown('Full vectors and all SI fields: [transport_packet_ledger.json](../checks/transport_packet_ledger.json).'))

Seed: 2026102202 | launch direction: [-0.3433989582523385, 0.760648603053623, 0.5509000436955697] | source weight: 1.226302883552434 | end: escaped_before_next_collision

<IPython.core.display.Markdown object>

7 Reconstruct a pulse from mass,with its own peak

Use retarded time $u=t_{\rm absolute}-t_{\rm onset}-R_0/c$.

s(u)\propto e^{-u/3ms}(1-e^{-u/\tau_r})H(u),\qquad y=s*h.

The continuous-source calibration below is independent of the reference sampled calibration. It places the direct atom at0 and scattered masses at native bin centers. These assumed intrinsic10–90%rises are not fitted heritage observations. Onset-to-peak and10–90%rise remain separate. The 3 ms is a decay time,not a measured duration.

t,y,ts,source,pulse=tr.independent_pulse(masses.mean(axis=0),direct0,80.,native_dt)
stored=next(x for x in case['timing']if x['intrinsic_rise_10_90_us']==80 and x['analog_one_pole_fc_hz']==0)
pulse_check={**pulse,'reference_10_90_us':stored['rise_10_90_us'],
             'difference_10_90_us':pulse['rise_10_90_us']-stored['rise_10_90_us'],
             'reference_onset_to_peak_us':stored['onset_to_peak_us'],
             'scope':'Independent continuous calibration; ideal electronics; no noise/classifier validation.'}
table([pulse_check],max_rows=1)
(ROOT/'checks/transport_convolution_check.json').write_text(json.dumps(pulse_check,indent=2))
assert abs(pulse_check['difference_10_90_us'])<.2
fig,ax=plt.subplots(layout='constrained')
ax.plot(ts*1e3,source,label='Assumed source; unit own peak')
ax.plot(t*1e3,y/y.max(),label='Transported HG pulse; unit own peak')
ax.set(xlim=(0,.8),xlabel='Retarded time after known model onset (ms)',ylabel='Intensity / each curve’s own peak',
       title='Independent mass convolution — shape normalization hides photometric gain')
ax.legend();plt.show()
print('Separate observer peak gain relative to vacuum unit-peak source:',pulse['peak_gain'])

Separate observer peak gain relative to vacuum unit-peak source: 1.9809545738273258

8 Convergence is a result,not a checkbox

The plotted uncertainty is±1finite-batch SEM,not a guaranteed confidence interval. A rare forward path absent from a batch cannot be represented by that batch’s variance. GHM remains unresolved; differences below must not be promoted to a physical GHM/HG rise-time bias.

comparison=tr.convergence_comparison(inputs['validation'])
table(comparison,['phase','primary_packets','primary_gain','primary_SE','larger_packets','larger_gain','larger_SE','interpretation'],max_rows=2)
fig,axes=plt.subplots(1,2,figsize=(9,3.5),sharey=True,layout='constrained')
for ax,row in zip(axes,comparison):
    ax.errorbar([row['primary_packets']/1e6,row['larger_packets']/1e6],
                [row['primary_gain'],row['larger_gain']],
                yerr=[row['primary_SE'],row['larger_SE']],fmt='o',capsize=4)
    ax.set(title=row['phase'].upper()+(' — convergence incomplete'if row['phase']=='ghm'else' — conditional stable control'),
           xlabel='Transport histories (million)')
axes[0].set_ylabel('Point-nadir relative fluence; ±1 batch SEM')
plt.show()
table(inputs['validation']['matched_CloudScat_HG']['cases'],
      ['reference_case','reference_gain_mean','new_gain_mean','difference_in_combined_SE','qualification'],max_rows=2)

Code and limits

Function locations and hashes below map the equations to inspectable source. Original engine/runner snapshots are listed separately from later diagnostic versions. Primary CloudScat method: https://gmd.copernicus.org/articles/13/5549/2020/ . A full RT rerun is optional and requires lawful raw ice input plus the separate reproduction launcher. It is deliberately not triggered here. Physical microphysics/structure,atmospheric paths,source spectrum/waveform and noisy labelled-event inference remain open.

table(report['code_map'],['file','function','line','end_line','sha256'],max_rows=30)
print(report['status'])
print('Not recomputed:',report['not_recomputed'])
display(Markdown('Derivation: [transport.md](../methods/transport.md). Native reconciliation: [CSV](../checks/transport_native_reconciliation.csv). Machine receipt: [JSON](../checks/transport_review.json).'))

LIGHT_RECOMPUTATION_PASS_WITH_EXPECTED_GHM_CONVERGENCE_WARNING
Not recomputed: ['GHM phase normalization without raw LUT', 'Physical cloud/atmosphere/absorption validity', 'Noisy onset or population classifier', 'Per-packet geometry,not retained in kernels']

<IPython.core.display.Markdown object>