fix: repair OLED9 TMM plane-wave sweep and add tests

.github/workflows/oled-tmm-tests.yml

@@ -0,0 +1,37 @@
+name: OLED TMM numerical checks
+
+on:
+  pull_request:
+    paths:
+      - 'tmm_OLED9.py'
+      - 'tests/test_tmm_oled9.py'
+      - 'requirements-tmm-sweep.txt'
+      - '.github/workflows/oled-tmm-tests.yml'
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+jobs:
+  numerical-tests:
+    runs-on: ubuntu-latest
+    timeout-minutes: 15
+    strategy:
+      fail-fast: false
+      matrix:
+        python-version: ['3.11', '3.12']
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+          cache: pip
+          cache-dependency-path: requirements-tmm-sweep.txt
+      - name: Install dependencies
+        run: python -m pip install -r requirements-tmm-sweep.txt
+      - name: Compile Python sources
+        run: python -m py_compile tmm_OLED9.py tests/test_tmm_oled9.py
+      - name: Run numerical smoke and regression tests
+        run: python -m pytest -q tests/test_tmm_oled9.py
+      - name: Run headless plotting smoke test
+        run: python tmm_OLED9.py --wavelength-count 3 --angle-count 3 --output /tmp/oled-transmission-smoke.png

requirements-tmm-sweep.txt

@@ -0,0 +1,4 @@
+numpy>=1.24,<3
+matplotlib>=3.7,<4
+tmm>=0.2.0,<1
+pytest>=7,<10

tests/test_tmm_oled9.py

@@ -0,0 +1,38 @@
+import importlib.util
+from pathlib import Path
+import numpy as np
+import pytest
+from tmm import coh_tmm
+
+spec = importlib.util.spec_from_file_location("oled9", Path(__file__).resolve().parents[1] / "tmm_OLED9.py")
+module = importlib.util.module_from_spec(spec)
+spec.loader.exec_module(module)
+
+
+def test_stack_has_matching_lengths():
+    n, d = module.stack_at(550.0, 75.0)
+    assert len(n) == len(d) == 4
+    assert np.isinf(d[0]) and np.isinf(d[-1])
+
+
+def test_transmission_matches_reference_tmm():
+    n, d = module.stack_at(550.0, 75.0)
+    expected = coh_tmm("s", n, d, np.deg2rad(30.0), 550.0)["T"]
+    assert module.transmission(550.0, 30.0, 75.0, "s") == pytest.approx(expected)
+
+
+def test_sweep_shape_and_finite_values():
+    result = module.sweep([450.0, 550.0], [0.0, 30.0], 75.0)
+    assert result.shape == (2, 2)
+    assert np.isfinite(result).all()
+    assert np.all(result >= 0)
+
+
+@pytest.mark.parametrize("wavelength,angle,thickness", [(0, 0, 50), (550, 90, 50), (550, 0, -1)])
+def test_invalid_inputs(wavelength, angle, thickness):
+    with pytest.raises(ValueError):
+        module.transmission(wavelength, angle, thickness, "s")
+
+
+def test_no_plot_on_import():
+    assert not hasattr(module, "plt")

tmm_OLED9.py

@@ -1,80 +1,77 @@
-import numpy as np
-import matplotlib
-matplotlib.use('Agg')  # Use Agg backend for headless environments
-import matplotlib.pyplot as plt
-
-# Dispersion relationships for materials
-def n_alq3(wavelength):
-    A = 1.65
-    B = 0.5  # Adjust as necessary
-    return A + B / (wavelength ** 2)
-
-def n_silver(wavelength):
-    n_0 = 0.16  # Real part
-    kappa = 3.16  # Imaginary part
-    C = 1.5  # Coefficient for the dispersion
-    return n_0 - (C / (wavelength ** 2)) + 1j * kappa
-
-# Wavelength and angle ranges
-wavelengths = np.linspace(400, 700, 1000)  # Wavelength range (in nm)
-theta_list_deg = np.linspace(0, 90, 180)  # Angles from 0 to 90 degrees
-theta_list_rad = np.radians(theta_list_deg)  # Convert to radians
-
-# Set thickness range for the organic layer (e.g., Alq3) in nanometers
-thicknesses = np.linspace(25, 150, 10)  # Example thickness values, from 50nm to 200nm
-
-# Loop over different thickness values of the organic layer
-for thickness in thicknesses:
-    # Initialize the outcoupling efficiency array for this thickness
-    outcoupling_efficiency_2D = np.zeros((len(wavelengths), len(theta_list_deg)))
-
-    # Create the n_list for the current wavelength range and thickness
-    n_air = 1.0  # Refractive index of air
-    n_list = [n_air]  # Start with air
-    for wavelength in wavelengths:
-        n_list.append(n_alq3(wavelength))  # Alq3 layer refractive index
-    n_list.append(n_silver(wavelengths[-1]))  # Silver refractive index at last wavelength
-
-    # Set the layer thickness values
-    d_alq3 = thickness * 1e-9  # Convert nm to meters
-    d_silver = 50e-9  # 50 nm for Silver
-    d_list = [np.inf, d_alq3, d_silver, np.inf]  # Air, Alq3, Silver, and back air
-
-    # Loop over each wavelength and angle to compute the outcoupling efficiency
-    for i, wavelength in enumerate(wavelengths):
-        R_TE, T_TE, R_TM, T_TM = [], [], [], []
-
-        for theta in theta_list_rad:
-            # Ensure angle is valid for first layer (Alq3)
-            if theta < np.arcsin(n_air / np.abs(n_list[1])):
-                # TE polarization (s-polarized light)
-                coh_TE = coh_tmm('s', n_list, d_list, wavelength * 1e-9, theta)
-                R_TE.append(coh_TE['R'])
-                T_TE.append(coh_TE['T'])
-
-                # TM polarization (p-polarized light)
-                coh_TM = coh_tmm('p', n_list, d_list, wavelength * 1e-9, theta)
-                R_TM.append(coh_TM['R'])
-                T_TM.append(coh_TM['T'])
-            else:
-                R_TE.append(1.0)  # Total reflection
-                T_TE.append(0.0)  # No transmission
-                R_TM.append(1.0)  # Total reflection
-                T_TM.append(0.0)  # No transmission
-
-        # Combine TE and TM transmission efficiencies to compute the outcoupling efficiency
-        outcoupling_efficiency = np.array(T_TE) + np.array(T_TM)
-        outcoupling_efficiency_2D[i, :] = outcoupling_efficiency
-
-    # Plotting the outcoupling efficiency for this thickness as a heatmap
-    plt.figure(figsize=(10, 6))
-    plt.contourf(theta_list_deg, wavelengths, outcoupling_efficiency_2D, levels=50, cmap='viridis')
-    plt.colorbar(label='Outcoupling Efficiency')
-    plt.xlabel('Angle of Incidence (degrees)')
-    plt.ylabel('Wavelength (nm)')
-    plt.title(f'Outcoupling Efficiency (Alq3 Thickness = {thickness:.1f} nm)')
-
-    # Save the plot as an image file
-    plt.savefig(f'outcoupling_efficiency_thickness_{thickness:.1f}nm.png')
-    plt.close()  # Close the figure to free up memory
-
+"""Plane-wave multilayer transmission sweep (not OLED dipole outcoupling)."""
+from pathlib import Path
+import argparse
+import numpy as np
+from tmm import coh_tmm
+
+
+def n_alq3(wavelength_nm):
+    """Illustrative dispersion only; not calibrated material data."""
+    return 1.65 + 0.5 / wavelength_nm**2
+
+
+def n_silver(wavelength_nm):
+    """Illustrative complex refractive index; replace with measured data."""
+    return 0.16 - 1.5 / wavelength_nm**2 + 3.16j
+
+
+def stack_at(wavelength_nm, organic_nm, silver_nm=50.0):
+    if wavelength_nm <= 0 or organic_nm < 0 or silver_nm < 0:
+        raise ValueError("Wavelength must be positive and thicknesses nonnegative")
+    return ([1.0, n_alq3(wavelength_nm), n_silver(wavelength_nm), 1.0],
+            [np.inf, organic_nm, silver_nm, np.inf])
+
+
+def transmission(wavelength_nm, angle_deg, organic_nm, polarization):
+    if not 0 <= angle_deg < 90:
+        raise ValueError("Angle must be in [0, 90) degrees")
+    if polarization not in ("s", "p"):
+        raise ValueError("Polarization must be s or p")
+    indices, thicknesses = stack_at(wavelength_nm, organic_nm)
+    result = coh_tmm(polarization, indices, thicknesses,
+                     np.deg2rad(angle_deg), wavelength_nm)
+    return float(result["T"])
+
+
+def sweep(wavelengths_nm, angles_deg, organic_nm):
+    wavelengths_nm = np.asarray(wavelengths_nm, dtype=float)
+    angles_deg = np.asarray(angles_deg, dtype=float)
+    if wavelengths_nm.ndim != 1 or angles_deg.ndim != 1:
+        raise ValueError("Sweep coordinates must be 1-D")
+    if not np.all(np.isfinite(wavelengths_nm)) or not np.all(np.isfinite(angles_deg)):
+        raise ValueError("Sweep coordinates must be finite")
+    values = np.empty((len(wavelengths_nm), len(angles_deg)))
+    for i, wavelength in enumerate(wavelengths_nm):
+        for j, angle in enumerate(angles_deg):
+            values[i, j] = 0.5 * (
+                transmission(wavelength, angle, organic_nm, "s")
+                + transmission(wavelength, angle, organic_nm, "p"))
+    return values
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--wavelength-count", type=int, default=31)
+    parser.add_argument("--angle-count", type=int, default=19)
+    parser.add_argument("--thickness-nm", type=float, default=75.0)
+    parser.add_argument("--output", type=Path, default=Path("transmission.png"))
+    args = parser.parse_args()
+    if args.wavelength_count < 2 or args.angle_count < 2:
+        parser.error("sample counts must be >= 2")
+    import matplotlib
+    matplotlib.use("Agg")
+    import matplotlib.pyplot as plt
+    wavelengths = np.linspace(400.0, 700.0, args.wavelength_count)
+    angles = np.linspace(0.0, 89.0, args.angle_count)
+    values = sweep(wavelengths, angles, args.thickness_nm)
+    fig, ax = plt.subplots()
+    plot = ax.contourf(angles, wavelengths, values, levels=30)
+    fig.colorbar(plot, ax=ax, label="Unpolarized plane-wave transmittance")
+    ax.set(xlabel="Incidence angle (degrees)", ylabel="Vacuum wavelength (nm)",
+           title=f"Illustrative stack: organic {args.thickness_nm:g} nm")
+    fig.savefig(args.output, dpi=150)
+    plt.close(fig)
+
+
+if __name__ == "__main__":
+    main()