Double Integrating Sphere Derivation and Parity Checks

This notebook documents the exact two-sphere forward model implemented in iadpython.double.DoubleSphere and the CWEB/C formulas in iad/src/iad_calc.c.

The goals are:

  • derive the exact single-sphere and two-sphere gains used by the code

  • show how unscattered collection and lost-light corrections are applied before two-sphere normalization

  • verify the notebook formulas against iadpython and the repo-local libiad

Assumptions and scope

  • The command-line interfaces treat sphere measurements as substitution mode by default, and they reject comparison mode with two spheres.

  • The implementation supports baffled and unbaffled reflection and transmission spheres. Baffled geometry is still the recommended experimental case.

  • For two spheres, M_R and M_T are computed from corrected UR1, UT1, URU, and UTU. Any unscattered collection adjustment and lost-light subtraction happens before the two-sphere normalization.

  • This notebook validates exact algebra against iadpython.double.DoubleSphere and against the repo-local libiad via iadpython.iadc.

[1]:

%config InlineBackend.figure_format = "retina"

import pathlib
import numpy as np

import iadpython
import iadpython.iadc as iadc

repo = pathlib.Path.cwd().resolve()
if not (repo / "iad").exists():
    repo = repo.parent
libiad_path = pathlib.Path(str(iadc.libiad_path))
expected_lib_dir = repo / "iad" / "src"

assert expected_lib_dir.exists(), f"repo root not found from {pathlib.Path.cwd()}"
assert libiad_path.exists(), f"libiad not found: {libiad_path}"
assert (
    libiad_path.resolve().parent == expected_lib_dir.resolve()
), f"expected repo-local libiad in {expected_lib_dir}, got {libiad_path.resolve()}"

print(f"Using libiad: {libiad_path}")

Using libiad: /Users/prahl/Documents/Code/git/iadpython/iad/src/libiad.dylib

Exact gain formulas

For one sphere, the code uses the exact Gain() algebra from CWEB.

For an unbaffled sphere,

\[G = \frac{1}{1 - a_w r_w - a_d r_d - a_s r_s - a_t r_t}.\]

For a baffled sphere,

\[G = \frac{1}{1 - \left(r_w + \frac{a_t}{a_w} r_t\right) \left[a_w + (1-a_t)(a_d r_d + a_s r_s)\right]}.\]

For two spheres, the exact coupling term is

\[X = a_s^r a_s^t a_w^r a_w^t (1-a_t^r)(1-a_t^t) G_r G_t \, UTU^2.\]

The implemented gains are then

\[G_{11} = \frac{G_r}{1 - X}, \qquad G_{22} = \frac{G_t}{1 - X}.\]

This exact coupling is stronger than the old shorthand a_s^2 G_r G_t UTU^2; the shorthand omits wall-area and third-port factors and is therefore not the implemented model.

[2]:

def gain_single_formula(sphere, uru_sample, uru_third):
    as_x = sphere.sample.a
    ad_x = sphere.detector.a
    at_x = sphere.third.a
    aw_x = sphere.a_wall
    rd_x = sphere.detector.uru
    rw_x = sphere.r_wall

    if sphere.baffle:
        inv_gain = rw_x + (at_x / aw_x) * uru_third
        inv_gain *= aw_x + (1 - at_x) * (ad_x * rd_x + as_x * uru_sample)
        inv_gain = 1.0 - inv_gain
    else:
        inv_gain = 1.0 - aw_x * rw_x - ad_x * rd_x - as_x * uru_sample - at_x * uru_third

    return 1.0 / inv_gain


def gain_11_formula(double_sphere, uru, utu):
    g_r = gain_single_formula(double_sphere.r_sphere, uru, 0.0)
    g_t = gain_single_formula(double_sphere.t_sphere, uru, 0.0)
    coupling = (
        double_sphere.r_sphere.sample.a
        * double_sphere.t_sphere.sample.a
        * double_sphere.r_sphere.a_wall
        * double_sphere.t_sphere.a_wall
        * (1 - double_sphere.r_sphere.third.a)
        * (1 - double_sphere.t_sphere.third.a)
        * g_r
        * g_t
        * utu
        * utu
    )
    return g_r / (1.0 - coupling)


def gain_22_formula(double_sphere, uru, utu):
    g_r = gain_single_formula(double_sphere.r_sphere, uru, 0.0)
    g_t = gain_single_formula(double_sphere.t_sphere, uru, 0.0)
    coupling = (
        double_sphere.r_sphere.sample.a
        * double_sphere.t_sphere.sample.a
        * double_sphere.r_sphere.a_wall
        * double_sphere.t_sphere.a_wall
        * (1 - double_sphere.r_sphere.third.a)
        * (1 - double_sphere.t_sphere.third.a)
        * g_r
        * g_t
        * utu
        * utu
    )
    return g_t / (1.0 - coupling)


def two_sphere_r_formula(double_sphere, ur1, uru, ut1, utu):
    g_t = gain_single_formula(double_sphere.t_sphere, uru, 0.0)
    power = double_sphere.r_sphere.detector.a
    power *= (1 - double_sphere.r_sphere.third.a) * double_sphere.r_sphere.r_wall
    power *= gain_11_formula(double_sphere, uru, utu)
    power *= (
        (1 - double_sphere.f_r) * ur1
        + double_sphere.r_sphere.r_wall * double_sphere.f_r
        + (1 - double_sphere.f_r)
        * double_sphere.t_sphere.sample.a
        * (1 - double_sphere.t_sphere.third.a)
        * double_sphere.t_sphere.r_wall
        * ut1
        * utu
        * g_t
    )
    return power


def two_sphere_t_formula(double_sphere, ur1, uru, ut1, utu):
    g_r = gain_single_formula(double_sphere.r_sphere, uru, 0.0)
    power = double_sphere.t_sphere.detector.a
    power *= (1 - double_sphere.t_sphere.third.a) * double_sphere.t_sphere.r_wall
    power *= gain_22_formula(double_sphere, uru, utu)
    power *= (1 - double_sphere.f_r) * ut1 + (
        1 - double_sphere.r_sphere.third.a
    ) * double_sphere.r_sphere.r_wall * double_sphere.r_sphere.sample.a * utu * (
        double_sphere.f_r * double_sphere.r_sphere.r_wall + (1 - double_sphere.f_r) * ur1
    ) * g_r
    return power


def measured_rt_formula(double_sphere, ur1, uru, ut1, utu):
    r_0 = two_sphere_r_formula(double_sphere, 0.0, 0.0, 0.0, 0.0)
    t_0 = two_sphere_t_formula(double_sphere, 0.0, 0.0, 0.0, 0.0)

    m_r = double_sphere.r_sphere.r_std
    m_r *= two_sphere_r_formula(double_sphere, ur1, uru, ut1, utu) - r_0
    m_r /= (
        two_sphere_r_formula(
            double_sphere,
            double_sphere.r_sphere.r_std,
            double_sphere.r_sphere.r_std,
            0.0,
            0.0,
        )
        - r_0
    )

    m_t = two_sphere_t_formula(double_sphere, ur1, uru, ut1, utu) - t_0
    m_t /= two_sphere_t_formula(double_sphere, 0.0, 0.0, 1.0, 1.0) - t_0
    return m_r, m_t


def raw_rt_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 corrected_two_sphere_inputs(experiment):
    ur1, ut1, uru, utu, ru, tu = raw_rt_inputs(experiment.sample)
    uru_calc = max(uru - experiment.uru_lost, 0.0)
    utu_calc = max(utu - experiment.utu_lost, 0.0)
    ur1_calc = max(ur1 - (1 - experiment.fraction_of_rc_in_mr) * ru - experiment.ur1_lost, 0.0)
    ut1_calc = max(ut1 - (1 - experiment.fraction_of_tc_in_mt) * tu - experiment.ut1_lost, 0.0)
    return (ur1_calc, uru_calc, ut1_calc, utu_calc), (ur1, ut1, uru, utu, ru, tu)


def make_spheres(sphere_config, config):
    r_sphere = iadpython.Sphere(refl=True, **sphere_config)
    t_sphere = iadpython.Sphere(refl=False, **sphere_config)
    r_sphere.baffle = config["r_baffle"]
    t_sphere.baffle = config["t_baffle"]
    return r_sphere, t_sphere


def make_experiment(sphere_config, optics, sample_kwargs, config, lost=None):
    sample = iadpython.Sample(**sample_kwargs, **optics)
    r_sphere, t_sphere = make_spheres(sphere_config, config)
    exp = iadpython.Experiment(
        sample=sample,
        num_spheres=2,
        r_sphere=r_sphere,
        t_sphere=t_sphere,
    )
    exp.method = "substitution"
    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 lost is not None:
        exp.ur1_lost = lost["ur1_lost"]
        exp.ut1_lost = lost["ut1_lost"]
        exp.uru_lost = lost["uru_lost"]
        exp.utu_lost = lost["utu_lost"]
    return exp

The exact coupling matters

The old shorthand coupling term,

a_s^2 * G_r * G_t * UTU^2,

is not the implemented model. The exact code uses

a_s^r * a_s^t * a_w^r * a_w^t * (1-a_t^r) * (1-a_t^t) * G_r * G_t * UTU^2.

For realistic sphere geometries these are close, but not identical.

[3]:

demo_sphere_config = {
    "d_sphere": 100.0,
    "d_sample": 40.0,
    "d_third": 10.0,
    "d_detector": 5.0,
    "r_detector": 0.04,
    "r_wall": 0.98,
    "r_std": 0.99,
}
demo_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,
}

r_sphere, t_sphere = make_spheres(demo_sphere_config, demo_config)
demo_double = iadpython.DoubleSphere(r_sphere, t_sphere)
uru_demo = 0.5
utu_demo = 0.2

simple_alpha = r_sphere.sample.a**2 * r_sphere.gain(uru_demo) * t_sphere.gain(uru_demo) * utu_demo**2
exact_alpha = (
    demo_double.r_sphere.sample.a
    * demo_double.t_sphere.sample.a
    * demo_double.r_sphere.a_wall
    * demo_double.t_sphere.a_wall
    * (1 - demo_double.r_sphere.third.a)
    * (1 - demo_double.t_sphere.third.a)
    * gain_single_formula(demo_double.r_sphere, uru_demo, 0.0)
    * gain_single_formula(demo_double.t_sphere, uru_demo, 0.0)
    * utu_demo**2
)

print(f"simplified alpha = {simple_alpha:.12f}")
print(f"exact alpha      = {exact_alpha:.12f}")
print(f"difference       = {simple_alpha - exact_alpha:.12f}")

simplified alpha = 0.035874310639
exact alpha      = 0.032682840103
difference       = 0.003191470536

Pre-normalization correction path

For Experiment.measured_rt() with two spheres, the notebook and the implementation use the same order:

  1. Compute raw UR1, UT1, URU, and UTU from the slab.

  2. Compute R_u and T_u from the slab boundaries.

  3. Subtract diffuse lost light from URU and UTU.

  4. Subtract the uncollected direct terms and direct lost light from UR1 and UT1.

  5. Clip the corrected values at zero.

  6. Feed the corrected values into the two-sphere normalization.

That matches the Python forward path in Experiment.measured_rt() and the CWEB path in Calculate_Distance_With_Corrections().

[4]:

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,
}

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,
        },
        "lost": {"ur1_lost": 0.0, "ut1_lost": 0.0, "uru_lost": 0.0, "utu_lost": 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,
        },
        "lost": {"ur1_lost": 0.012, "ut1_lost": 0.017, "uru_lost": 0.006, "utu_lost": 0.010},
    },
    {
        "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,
        },
        "lost": {"ur1_lost": 0.018, "ut1_lost": 0.013, "uru_lost": 0.009, "utu_lost": 0.011},
    },
)

for case in CASES:
    exp = make_experiment(SPHERE_CONFIG, case["optics"], case["sample"], case["config"], case["lost"])
    corrected, raw = corrected_two_sphere_inputs(exp)
    ur1_calc, uru_calc, ut1_calc, utu_calc = corrected
    ur1, ut1, uru, utu, ru, tu = raw

    ds = iadpython.DoubleSphere(exp.r_sphere, exp.t_sphere)
    ds.f_r = exp.f_r

    g11_formula = gain_11_formula(ds, uru_calc, utu_calc)
    g22_formula = gain_22_formula(ds, uru_calc, utu_calc)
    g11_python, g22_python = ds.gain(uru_calc, utu_calc)

    assert abs(g11_formula - g11_python) < 1e-12
    assert abs(g22_formula - g22_python) < 1e-12

    c_gain_r = iadc._c_gain(exp.r_sphere, uru_calc, 0.0)
    c_gain_t = iadc._c_gain(exp.t_sphere, uru_calc, 0.0)
    assert abs(gain_single_formula(exp.r_sphere, uru_calc, 0.0) - c_gain_r) < 1e-12
    assert abs(gain_single_formula(exp.t_sphere, uru_calc, 0.0) - c_gain_t) < 1e-12

    mr_formula, mt_formula = measured_rt_formula(ds, ur1_calc, uru_calc, ut1_calc, utu_calc)
    mr_python, mt_python = exp.measured_rt()
    mr_c, mt_c = iadc._c_calculate_measured_rt_from_rt(exp, ur1, ut1, uru, utu, ru, tu)

    assert abs(float(mr_formula) - float(mr_python)) < 1e-10
    assert abs(float(mt_formula) - float(mt_python)) < 1e-10
    assert abs(float(mr_formula) - float(mr_c)) < 1e-10
    assert abs(float(mt_formula) - float(mt_c)) < 1e-10

    print(case["name"])
    print(f"  corrected inputs : UR1={ur1_calc:.6f} URU={uru_calc:.6f} UT1={ut1_calc:.6f} UTU={utu_calc:.6f}")
    print(f"  gains            : G11={g11_formula:.12f} G22={g22_formula:.12f}")
    print(f"  measured RT      : MR={float(mr_formula):.12f} MT={float(mt_formula):.12f}")
    print()

print("All parity checks passed.")

matched_unbaffled
  corrected inputs : UR1=0.134858 URU=0.191912 UT1=0.340226 UTU=0.226070
  gains            : G11=44.998759905498 G22=44.998759905498
  measured RT      : MR=0.133108152900 MT=0.344200628360

slab_in_air_mixed_collection
  corrected inputs : UR1=0.059863 URU=0.124042 UT1=0.157493 UTU=0.123036
  gains            : G11=44.836266714048 G22=44.842581457067
  measured RT      : MR=0.058752142282 MT=0.158162418362

glass_slides_baffled
  corrected inputs : UR1=0.058594 URU=0.136279 UT1=0.139992 UTU=0.123752
  gains            : G11=44.863059323607 G22=44.863059323607
  measured RT      : MR=0.057948603099 MT=0.140730648364

All parity checks passed.