Validation of Double Sphere Formulas
Scott Prahl
April 2026
This notebook replaces the older exploratory draft with a focused validation workflow for the two-sphere model in iadpython.
It checks three things:
representative two-sphere detected powers agree with Monte Carlo,
the direct two-sphere
M_RandM_Tformulas match the C implementation inlibiad,Python and C inversions agree in measurement space for representative boundary conditions.
The parity sections expect a repo-local shared library built with make libiad.dylib in iad/src.
[1]:
%config InlineBackend.figure_format='retina'
import pathlib
import sys
from pprint import pprint
import numpy as np
def find_repo_root(start=None):
if start is None:
start = pathlib.Path.cwd().resolve()
for candidate in (start, *start.parents):
if (candidate / "iadpython").is_dir() and (candidate / "docs").is_dir():
return candidate
raise RuntimeError("Run this notebook from within the iadpython repository.")
REPO_ROOT = find_repo_root()
if str(REPO_ROOT) not in sys.path:
sys.path.insert(0, str(REPO_ROOT))
import iadpython
import iadpython.double as iad_double
import iadpython.iadc as iadc
print(f"repo root: {REPO_ROOT}")
print(f"libiad: {iadc.libiad_path}")
repo root: /Users/prahl/Documents/Code/git/iadpython
libiad: /Users/prahl/Documents/Code/git/iadpython/iad/src/libiad.dylib
[2]:
SPHERE_CONFIG = {
"d_sphere": 150.0,
"d_sample": 10.0,
"d_third": 10.0,
"d_detector": 5.0,
"r_detector": 0.04,
"r_wall": 0.98,
"r_std": 0.99,
}
MC_CASES = (
{
"name": "unbaffled moderate coupling",
"setup": {
"ur1": 0.50,
"ut1": 0.40,
"uru": 0.30,
"utu": 0.50,
"r_baffle": False,
"t_baffle": False,
"f_r": 0.0,
},
"n_photons": 10_000,
"seed": 20260411,
},
{
"name": "baffled mixed coupling",
"setup": {
"ur1": 0.55,
"ut1": 0.30,
"uru": 0.40,
"utu": 0.35,
"r_baffle": True,
"t_baffle": True,
"f_r": 0.0,
},
"n_photons": 10_000,
"seed": 20260411,
},
)
FORWARD_CASES = (
{
"name": "matched unbaffled",
"optics": {"a": 0.68, "b": 1.70, "g": 0.25},
"sample": {"n": 1.0, "n_above": 1.0, "n_below": 1.0, "quad_pts": 16},
"config": {
"r_baffle": False,
"t_baffle": False,
"fraction_of_rc_in_mr": 1.0,
"fraction_of_tc_in_mt": 1.0,
"f_r": 0.0,
},
},
{
"name": "slab in air mixed collection",
"optics": {"a": 0.82, "b": 4.10, "g": 0.72},
"sample": {"n": 1.4, "n_above": 1.0, "n_below": 1.0, "quad_pts": 16},
"config": {
"r_baffle": True,
"t_baffle": False,
"fraction_of_rc_in_mr": 0.78,
"fraction_of_tc_in_mt": 0.83,
"f_r": 0.15,
},
},
{
"name": "glass slides baffled",
"optics": {"a": 0.66, "b": 2.30, "g": 0.38},
"sample": {"n": 1.4, "n_above": 1.5, "n_below": 1.5, "quad_pts": 16},
"config": {
"r_baffle": True,
"t_baffle": True,
"fraction_of_rc_in_mr": 0.64,
"fraction_of_tc_in_mt": 0.71,
"f_r": 0.22,
},
},
{
"name": "matched baffled partial collection",
"optics": {"a": 0.74, "b": 2.60, "g": 0.55},
"sample": {"n": 1.0, "n_above": 1.0, "n_below": 1.0, "quad_pts": 16},
"config": {
"r_baffle": True,
"t_baffle": True,
"fraction_of_rc_in_mr": 0.62,
"fraction_of_tc_in_mt": 0.47,
"f_r": 0.10,
},
},
)
INVERSE_CASES = (
{
"name": "matched unbaffled",
"optics": {"a": 0.68, "b": 1.70, "g": 0.25},
"sample": {"n": 1.0, "n_above": 1.0, "n_below": 1.0, "quad_pts": 16},
"config": {
"r_baffle": False,
"t_baffle": False,
"fraction_of_rc_in_mr": 1.0,
"fraction_of_tc_in_mt": 1.0,
"f_r": 0.0,
},
"default_g": 0.25,
},
{
"name": "slab in air mixed collection",
"optics": {"a": 0.82, "b": 4.10, "g": 0.72},
"sample": {"n": 1.4, "n_above": 1.0, "n_below": 1.0, "quad_pts": 16},
"config": {
"r_baffle": True,
"t_baffle": False,
"fraction_of_rc_in_mr": 0.78,
"fraction_of_tc_in_mt": 0.83,
"f_r": 0.15,
},
"default_g": 0.72,
},
{
"name": "glass slides baffled",
"optics": {"a": 0.66, "b": 2.30, "g": 0.38},
"sample": {"n": 1.4, "n_above": 1.5, "n_below": 1.5, "quad_pts": 16},
"config": {
"r_baffle": True,
"t_baffle": True,
"fraction_of_rc_in_mr": 0.64,
"fraction_of_tc_in_mt": 0.71,
"f_r": 0.22,
},
"default_g": 0.38,
},
)
def make_spheres(config):
r_sphere = iadpython.Sphere(refl=True, **SPHERE_CONFIG)
t_sphere = iadpython.Sphere(refl=False, **SPHERE_CONFIG)
r_sphere.baffle = config.get("r_baffle", False)
t_sphere.baffle = config.get("t_baffle", False)
return r_sphere, t_sphere
def default_double(**setup):
r_sphere, t_sphere = make_spheres(setup)
double = iadpython.DoubleSphere(r_sphere, t_sphere)
double.ur1 = setup["ur1"]
double.ut1 = setup["ut1"]
double.uru = setup["uru"]
double.utu = setup["utu"]
double.f_r = setup.get("f_r", 0.0)
return double
def make_sample(optics, sample_kwargs):
params = dict(sample_kwargs)
params.update(optics)
return iadpython.Sample(**params)
def make_experiment(sample, config, measured=None, default_g=None):
r_sphere, t_sphere = make_spheres(config)
exp = iadpython.Experiment(
sample=sample,
num_spheres=2,
r_sphere=r_sphere,
t_sphere=t_sphere,
)
exp.fraction_of_rc_in_mr = config["fraction_of_rc_in_mr"]
exp.fraction_of_tc_in_mt = config["fraction_of_tc_in_mt"]
exp.f_r = config["f_r"]
if measured is not None:
exp.m_r, exp.m_t = measured
exp.default_g = default_g
return exp
def measurement_inputs(sample):
ur1, ut1, uru, utu = sample.rt()
nu_inside = iadpython.cos_snell(1.0, sample.nu_0, sample.n)
ru, tu = iadpython.specular_rt(sample.n_above, sample.n, sample.n_below, sample.b, nu_inside)
return ur1, ut1, uru, utu, ru, tu
def assert_repo_libiad():
lib_path = pathlib.Path(str(iadc.libiad_path)).resolve()
expected_dir = REPO_ROOT / "iad" / "src"
assert (
lib_path.parent == expected_dir.resolve()
), "Build the repo-local library with `make libiad.dylib` in iad/src before running parity checks."
return lib_path
LIBIAD_PATH = assert_repo_libiad()
def run_monte_carlo_case(case, z_value=5.0):
iad_double.time.time = lambda: case.get("seed", 20260411)
dd = default_double(**case["setup"])
analytic_r = float(dd.two_sphere_r(dd.ur1, dd.uru, dd.ut1, dd.utu))
analytic_t = float(dd.two_sphere_t(dd.ur1, dd.uru, dd.ut1, dd.utu))
gain_11, gain_22 = dd.gain()
ave_r, stderr_r, ave_t, stderr_t = dd.do_N_photons(case["n_photons"])
result = {
"name": case["name"],
"gain_11": float(gain_11),
"gain_22": float(gain_22),
"analytic_r": analytic_r,
"analytic_t": analytic_t,
"mc_r": float(ave_r),
"mc_r_stderr": float(stderr_r),
"mc_t": float(ave_t),
"mc_t_stderr": float(stderr_t),
}
assert abs(result["analytic_r"] - result["mc_r"]) <= z_value * result["mc_r_stderr"], result
assert abs(result["analytic_t"] - result["mc_t"]) <= z_value * result["mc_t_stderr"], result
return result
def run_forward_parity(case):
sample = make_sample(case["optics"], case["sample"])
exp = make_experiment(sample, case["config"])
ur1, ut1, uru, utu, ru, tu = measurement_inputs(sample)
python_mr, python_mt = exp.measured_rt()
c_mr, c_mt = iadc._c_calculate_measured_rt_from_rt(exp, ur1, ut1, uru, utu, ru, tu)
result = {
"name": case["name"],
"python_mr": float(python_mr),
"c_mr": float(c_mr),
"python_mt": float(python_mt),
"c_mt": float(c_mt),
}
assert abs(result["python_mr"] - result["c_mr"]) <= 1e-10, result
assert abs(result["python_mt"] - result["c_mt"]) <= 1e-10, result
return result
def run_inverse_parity(case):
forward_sample = make_sample(case["optics"], case["sample"])
forward_exp = make_experiment(forward_sample, case["config"])
measured_r, measured_t = forward_exp.measured_rt()
inverse_sample = iadpython.Sample(**case["sample"])
inverse_exp = make_experiment(
inverse_sample,
case["config"],
measured=(float(measured_r), float(measured_t)),
default_g=case["default_g"],
)
py_a, py_b, py_g = inverse_exp.invert_rt()
c_a, c_b, c_g, c_error = iadc._c_invert_experiment(inverse_exp)
assert c_error == 0, {"name": case["name"], "c_error": c_error}
py_sample = make_sample({"a": py_a, "b": py_b, "g": py_g}, case["sample"])
c_sample = make_sample({"a": c_a, "b": c_b, "g": c_g}, case["sample"])
py_forward = make_experiment(py_sample, case["config"]).measured_rt()
c_forward = make_experiment(c_sample, case["config"]).measured_rt()
result = {
"name": case["name"],
"target_mr": float(measured_r),
"py_mr": float(py_forward[0]),
"c_mr": float(c_forward[0]),
"target_mt": float(measured_t),
"py_mt": float(py_forward[1]),
"c_mt": float(c_forward[1]),
}
assert abs(result["py_mr"] - result["target_mr"]) <= 6e-2, result
assert abs(result["py_mt"] - result["target_mt"]) <= 6e-2, result
assert abs(result["c_mr"] - result["target_mr"]) <= 6e-2, result
assert abs(result["c_mt"] - result["target_mt"]) <= 6e-2, result
return result
Monte Carlo Validation
DoubleSphere.do_N_photons() is still a noisy stochastic check, so this section keeps the Monte Carlo validation focused on two representative cases that finish in a reasonable amount of time. Each case uses a fixed random seed and a loose 5-sigma gate so the notebook stays reproducible while still catching large regressions.
The exact analytic quantities being checked are the detected powers from two_sphere_r() and two_sphere_t().
[3]:
mc_results = [run_monte_carlo_case(case) for case in MC_CASES]
for result in mc_results:
pprint(result)
print()
r_average detected = 0.006 ± 0.000
average gain = 21.367 ± 1.406
calculated gain = 45.266
t_average detected = 0.004 ± 0.000
average gain = 14.420 ± 1.064
calculated gain = 45.266
r_average detected = 0.006 ± 0.001
average gain = 24.726 ± 2.499
calculated gain = 45.460
t_average detected = 0.003 ± 0.000
average gain = 11.441 ± 1.063
calculated gain = 45.460
{'analytic_r': 0.006275470683819339,
'analytic_t': 0.00507488370815709,
'gain_11': 45.265881491096685,
'gain_22': 45.265881491096685,
'mc_r': 0.005697754577160003,
'mc_r_stderr': 0.00037490699176383294,
'mc_t': 0.0038452814553210915,
'mc_t_stderr': 0.000283619067553569,
'name': 'unbaffled moderate coupling'}
{'analytic_r': 0.006863009841612849,
'analytic_t': 0.00382608942454521,
'gain_11': 45.46027327734377,
'gain_22': 45.46027327734377,
'mc_r': 0.006454462277271494,
'mc_r_stderr': 0.0006523053794636887,
'mc_t': 0.002986694621985705,
'mc_t_stderr': 0.0002774186764735002,
'name': 'baffled mixed coupling'}
Direct Parity With libiad
The next check uses the same representative cases as the local regression tests. Here the Python and C implementations should agree to machine precision because both are evaluating the same two-sphere formulas from the same UR1, UT1, URU, UTU, R_u, and T_u inputs.
[4]:
print(f"Using repo-local libiad: {LIBIAD_PATH}")
forward_results = [run_forward_parity(case) for case in FORWARD_CASES]
for result in forward_results:
pprint(result)
print()
Using repo-local libiad: /Users/prahl/Documents/Code/git/iadpython/iad/src/libiad.dylib
{'c_mr': 0.13310815289954323,
'c_mt': 0.3442006283596898,
'name': 'matched unbaffled',
'python_mr': 0.13310815289954317,
'python_mt': 0.3442006283596898}
{'c_mr': 0.0703981788519456,
'c_mt': 0.17532911241059623,
'name': 'slab in air mixed collection',
'python_mr': 0.07039817885194567,
'python_mt': 0.1753291124105963}
{'c_mr': 0.07527311700128111,
'c_mt': 0.1539783139434328,
'name': 'glass slides baffled',
'python_mr': 0.07527311700128111,
'python_mt': 0.15397831394343278}
{'c_mr': 0.10203337889756095,
'c_mt': 0.24848590199102127,
'name': 'matched baffled partial collection',
'python_mr': 0.10203337889756095,
'python_mt': 0.24848590199102125}
Inverse Parity In Measurement Space
These cases verify the practical inversion gate: both implementations should recover optical properties that re-forward to the same measured two-sphere M_R and M_T values.
[5]:
inverse_results = [run_inverse_parity(case) for case in INVERSE_CASES]
for result in inverse_results:
pprint(result)
print()
{'c_mr': 0.13310526801503655,
'c_mt': 0.3441927771276344,
'name': 'matched unbaffled',
'py_mr': 0.13310526801503655,
'py_mt': 0.3441927771276344,
'target_mr': 0.13310815289954317,
'target_mt': 0.3442006283596898}
{'c_mr': 0.0703967517210653,
'c_mt': 0.17533091363743505,
'name': 'slab in air mixed collection',
'py_mr': 0.0703967517210653,
'py_mt': 0.17533091363743505,
'target_mr': 0.07039817885194567,
'target_mt': 0.1753291124105963}
{'c_mr': 0.07526848931776685,
'c_mt': 0.1539773851130476,
'name': 'glass slides baffled',
'py_mr': 0.07526848931776678,
'py_mt': 0.1539773851130476,
'target_mr': 0.07527311700128111,
'target_mt': 0.15397831394343278}
Result
If all cells above run without an assertion failure, then the two-sphere notebook is exercising the current working model rather than the older exploratory draft.