Add oracle-validation example for cross-polarized transmittance

CLAUDE.md

@@ -96,7 +96,9 @@ Pass via the `sheets=` kwarg (Dict or iterable of `i => sheet`) on
 - **Time convention**: `exp(-iωt)` (Berreman / Passler & Paarmann), opposite of
   Yeh. Flips the sign in the propagation matrix — matters when comparing formulas.
 - **Transmittance ≠ |t|²** in general: the transmitted wave is in a different
-  medium, so T needs a Poynting-vector ratio (`T = S_out,z / S_inc,z`).
+  medium, so T needs a Poynting-vector ratio (`T = S_out,z / S_inc,z`). This holds
+  for **all four** channels: `Tpp/Tss/Tps/Tsp` are each the Poynting flux of one
+  substrate eigenmode (evaluated with its own wavevector) over the incident flux.
   Reflectance uses `R = |r|²` directly.
 - **γᵢ₃₃ sign correction**: Eq. 13 of Xu et al. (2000) / Eq. 20 of Passler (2017)
   has a sign error in γᵢ₃₃ as printed. All γⱼ₃ must satisfy the z-constraint
@@ -111,15 +113,19 @@ conservation): `.claude/rules/berreman-4x4-equations.md`.
 
 - **#70** (resolved): rotated anisotropic crystals conserve energy. The quoted
   `Tpp+Rpp ≈ 0.998` (uniaxial `euler=(π/6,π/4,0)`) was a **budget artifact** — it
-  omits the reflected cross-pol `Rps` that an out-of-plane tilt produces. The correct
-  per-input-polarization budget is `Rpp+Rps+Tpp = 1` and `Rss+Rsp+Tss = 1` (holds to
-  ~1e-14; the Poynting `Tpp` already includes the cross-*transmitted* power, so do
-  **not** also add `Tps`/`Tsp`). Convention: `r_{in,out}` — p-input cross-reflection
-  is `Rps`, not `Rsp` (they coincide only at normal incidence). Separately, a genuine
-  bug — transmission into an **anisotropic substrate** (`Rpp+Rps+Tpp ≈ 1.017`) — was
+  omits the converted channels that an out-of-plane tilt produces. The correct
+  per-input-polarization budget is `Rpp+Rps+Tpp+Tps = 1` and `Rss+Rsp+Tss+Tsp = 1`
+  (holds to ~1e-14). Every transmittance is a **per-mode Poynting flux**: `Tpp` is
+  the p-like substrate eigenmode's power only, `Tps` the s-like eigenmode's, each
+  evaluated with its own wavevector. (Historical bug: `Tpp` used to be the TOTAL
+  per-input flux while `Tps = |tps|²` was an amplitude ratio — converting stacks
+  then double-counted, e.g. a lossless half-wave plate at 45° gave a p-input sum
+  of 1.85.) Convention: `r_{in,out}` — p-input cross-reflection is `Rps`, not
+  `Rsp` (they coincide only at normal incidence). Separately, a genuine bug —
+  transmission into an **anisotropic substrate** (`Rpp+Rps+Tpp ≈ 1.017`) — was
   fixed in `poynting()`: the two transmitted substrate eigenmodes carry different
-  wavevectors, so their Poynting vectors are now summed per-mode instead of using one
-  wavevector for the combined field.
+  wavevectors, so their Poynting vectors are evaluated per-mode instead of using
+  one wavevector for the combined field.
 - **#71**: anisotropic ambient at oblique incidence → NaN (uses nx for k_par).
 - **#72**: absorbing incident medium — |r|² is not a true energy reflectance. The
   Poynting vector is non-additive for absorbing incident media (interference

docs/src/examples/anisotropic_energy_budget.md

@@ -15,53 +15,57 @@ air  = Layer(λ -> 1.0, 1.0)
 slab = Layer(λ -> 1.658, λ -> 1.658, λ -> 1.486, 0.5; euler = (π/6, π/4, 0.0))
 
 r = transfer(1.0, [air, slab, air])          # normal incidence
-r.Rpp + r.Tpp        # ≈ 0.9978  — looks like 0.2% has vanished
+r.Rpp + r.Tpp        # ≈ 0.9873  — looks like 1.3% has vanished
 ```
 
-The slab is transparent, yet `Rpp + Tpp ≈ 0.998`. No energy is actually lost — the
+The slab is transparent, yet `Rpp + Tpp ≈ 0.987`. No energy is actually lost — the
 budget is simply incomplete.
 
-## The missing piece is polarization conversion
+## The missing pieces are the converted channels
 
 A rotated anisotropic crystal does something an isotropic film never does: it
 **converts polarization**. Send in a purely p‑polarized wave and a small
-s‑polarized component comes back, the cross‑polarized reflectance `Rps`. That
-reflected light is real, carries energy, and is exactly what the naive sum omits.
-Adding it closes the budget to machine precision:
+s‑polarized component appears in *both* directions — the cross‑polarized
+reflectance `Rps` and the cross‑polarized transmittance `Tps`. That converted
+light is real, carries energy, and is exactly what the naive sum omits. Adding
+both channels closes the budget to machine precision:
 
 ```julia
-r.Rpp + r.Rps + r.Tpp        # = 1.0   (to ~1e-15)
+r.Rpp + r.Rps + r.Tpp + r.Tps        # = 1.0   (to ~1e-15)
 ```
 
 A few conventions make this add up cleanly:
 
   - The coefficient is indexed `r_{in,out}`, so the **p‑input** cross‑reflection is
     `Rps` (not `Rsp`, which is the *s*-input term — the two coincide only at normal
     incidence).
-  - `Tpp` here is the Poynting‑flux transmittance, which already includes the
-    cross‑transmitted (p→s) light, so you do **not** also add `Tps`.
+  - Every transmittance is a **per‑mode Poynting flux**: `Tpp` is the power carried
+    by the p‑like substrate eigenmode alone, `Tps` the power carried by the s‑like
+    eigenmode, each evaluated with its own wavevector. The four channels therefore
+    partition the transmitted energy exactly — nothing is counted twice.
 
 The per‑input‑polarization statements of energy conservation are therefore
 
 ```math
-R_{pp} + R_{ps} + T_{pp} = 1, \qquad R_{ss} + R_{sp} + T_{ss} = 1 .
+R_{pp} + R_{ps} + T_{pp} + T_{ps} = 1, \qquad R_{ss} + R_{sp} + T_{ss} + T_{sp} = 1 .
 ```
 
 ![Energy budget of a rotated anisotropic slab](../assets/examples/anisotropic_energy_budget.png)
 
 **Panel (a)** sweeps the incidence angle with the optic axis held out of plane. The
 naive `Rpp + Tpp` (blue) dips below one at every angle; the complete
-`Rpp + Rps + Tpp` (black) is flat at one. The shaded gap between them *is* `Rps`.
-(It pinches shut near 64° where this particular geometry happens to stop converting,
-then reopens — `Rps` is itself angle‑dependent.)
+`Rpp + Rps + Tpp + Tps` (black) is flat at one. The shaded gap between them *is*
+the converted power `Rps + Tps`.
 
 **Panel (b)** explains *when* the conversion happens. Holding the incidence angle
-fixed and rotating the optic axis in azimuth, `Rps` is exactly zero when the axis
-lies in the plane of incidence (`α = 0`): there the mirror symmetry of the problem
-keeps p and s decoupled, just like an isotropic medium. Tilt the axis out of that
-plane and the conversion — and the apparent energy "leak" — switches on.
+fixed and rotating the optic axis in azimuth, the converted power is exactly zero
+when the axis lies in the plane of incidence (`α = 0`): there the mirror symmetry
+of the problem keeps p and s decoupled, just like an isotropic medium. Tilt the
+axis out of that plane and the conversion — and the apparent energy "leak" —
+switches on, mostly in transmission (`Tps`) with a smaller reflected share (`Rps`).
 
 The lesson generalizes: for any anisotropic stack, always account for the
-cross‑polarized channels before concluding that energy is or isn't conserved.
+cross‑polarized channels — both reflected *and* transmitted — before concluding
+that energy is or isn't conserved.
 
 The full runnable script is [`examples/anisotropic_energy_budget.jl`](https://github.com/garrekstemo/TransferMatrix.jl/blob/main/examples/anisotropic_energy_budget.jl).

docs/src/guide/transfer_matrix.md

@@ -68,9 +68,11 @@ r_{sp} &= \frac{M_{11}M_{23} - M_{21}M_{13}}{d}.
 \end{aligned}
 ```
 
-Reflectance is then ``R = |r|^2`` directly, while transmittance ``T`` is a
-Poynting-vector flux ratio rather than ``|t|^2`` (the transmitted wave lives in a
-different medium) — see [Physics Validation](validation.md). These formulas are exactly
+Reflectance is then ``R = |r|^2`` directly, while every transmittance ``T`` —
+co- and cross-polarized alike — is a Poynting-vector flux ratio rather than
+``|t|^2`` (the transmitted wave lives in a different medium): each channel is the
+flux of one substrate eigenmode, evaluated with its own wavevector — see
+[Physics Validation](validation.md). The amplitude formulas above are exactly
 what [`calculate_tr`](@ref) evaluates.
 
 The full ``16``-element map:

docs/src/guide/validation.md

@@ -70,13 +70,23 @@ Dividing by ``S_{\text{incident}}``:
 T = \frac{n_2 \cos\theta_2}{n_1 \cos\theta_1} |t|^2
 ```
 
-This is why TransferMatrix.jl uses Poynting vector calculations rather than simply squaring the transmission coefficient.
+This is why TransferMatrix.jl uses Poynting vector calculations rather than simply squaring the transmission coefficient — for **every** channel, including the cross-polarized ones (each transmittance is the Poynting flux of its own substrate eigenmode).
+
+**Anisotropic (polarization-converting) stacks:** a rotated birefringent or gyrotropic layer converts p ↔ s, so part of the incident power leaves through the cross-polarized channels. The budget checked per input polarization is therefore
+
+```math
+R_{pp} + R_{ps} + T_{pp} + T_{ps} = 1, \qquad R_{ss} + R_{sp} + T_{ss} + T_{sp} = 1 ,
+```
+
+which reduces to ``R + T = 1`` for isotropic stacks, where the cross terms vanish.
 
 **Example warning:**
 ```
 ┌ Warning: Energy conservation violated for p-polarization
 │   Tpp = 0.72
+│   Tps = 0.0
 │   Rpp = 0.25
+│   Rps = 0.0
 │   sum = 0.97
 │   expected = 1.0
 │   deviation = 0.03
@@ -142,7 +152,7 @@ Validation does **not** catch all possible errors:
 - **Field continuity:** Tangential E and H fields should be continuous across interfaces. This is tested in the test suite but not validated at runtime.
 - **Incorrect but consistent results:** If there's a systematic error that affects both R and T equally, energy conservation might still hold.
 - **Material data errors:** Validation can't detect if your refractive index data is wrong—only that the calculation is self-consistent.
-- **Anisotropic media issues:** The validation assumes isotropic materials. For birefringent media, additional checks may be needed.
+- **Per-channel errors that cancel in the budget:** The energy check sums all four channels per input polarization, so a systematic redistribution between co- and cross-polarized channels that preserves the total would not be flagged.
 
 ## Tolerance Parameters
 

examples/anisotropic_energy_budget.jl

@@ -3,22 +3,25 @@
 # A lossless dielectric slab must obey R + T = 1. But take a uniaxial crystal and
 # rotate its optic axis OUT of the plane of incidence, and the naive budget
 #
-#     R + T  =  Rpp + Tpp  ≈  0.998
+#     R + T  =  Rpp + Tpp  ≈  0.987
 #
 # falls short of 1 — even though nothing absorbs. The slab is not leaking energy:
 # a rotated anisotropic crystal CONVERTS polarization. A purely p-polarized input
-# comes back with a small s-polarized component (Rps), and that reflected
-# cross-polarized light is exactly the piece the naive budget forgets. The
+# comes back with a small s-polarized component in BOTH directions — the
+# cross-polarized reflectance Rps and the cross-polarized transmittance Tps.
+# Those converted channels are exactly what the naive budget forgets. The
 # complete, polarization-resolved budget
 #
-#     Rpp + Rps + Tpp  =  1
+#     Rpp + Rps + Tpp + Tps  =  1
 #
-# closes to machine precision. (Tpp from the Poynting flux already includes the
-# cross-transmitted light, so it is not added twice.)
+# closes to machine precision. (Each transmittance is the Poynting flux of its
+# own substrate eigenmode, so the four channels partition the energy exactly.)
 #
-# Panel (a): the naive vs. complete budget vs. incidence angle — the shaded gap IS Rps.
-# Panel (b): Rps vs. the optic-axis azimuth — conversion vanishes when the axis lies
-#            in the plane of incidence and grows as it tilts out.
+# Panel (a): the naive vs. complete budget vs. incidence angle — the shaded gap
+#            IS the converted power Rps + Tps.
+# Panel (b): the converted channels vs. the optic-axis azimuth — conversion
+#            vanishes when the axis lies in the plane of incidence and grows as
+#            it tilts out.
 #
 # Run from the examples environment:
 #   julia --project=examples examples/anisotropic_energy_budget.jl
@@ -41,46 +44,53 @@ res_θ = [transfer(λ, [air, slab(α0, β0), air]; θ = deg2rad(θ)) for θ in 
 Rpp = getfield.(res_θ, :Rpp)
 Tpp = getfield.(res_θ, :Tpp)
 Rps = getfield.(res_θ, :Rps)
-naive    = Rpp .+ Tpp            # forgets the p→s reflection
-complete = Rpp .+ Rps .+ Tpp     # = 1
+Tps = getfield.(res_θ, :Tps)
+naive    = Rpp .+ Tpp                    # forgets both converted channels
+complete = Rpp .+ Rps .+ Tpp .+ Tps      # = 1
 
 println("Panel (a): rotated slab euler=(π/6, π/4, 0)")
 println("  naive  Rpp+Tpp:      min=", round(minimum(naive), digits=5),
         "  max=", round(maximum(naive), digits=5), "  (should dip below 1)")
-println("  complete Rpp+Rps+Tpp: max |budget-1| = ", maximum(abs.(complete .- 1)))
+println("  complete Rpp+Rps+Tpp+Tps: max |budget-1| = ", maximum(abs.(complete .- 1)))
 
-# Panel (b): sweep the optic-axis azimuth at fixed incidence; watch Rps switch on.
+# Panel (b): sweep the optic-axis azimuth at fixed incidence; watch conversion switch on.
 αdeg = 0.0:1.0:90.0
 θ_fixed = deg2rad(30)
-Rps_α = [transfer(λ, [air, slab(deg2rad(α), β0), air]; θ = θ_fixed).Rps for α in αdeg]
+res_α = [transfer(λ, [air, slab(deg2rad(α), β0), air]; θ = θ_fixed) for α in αdeg]
+Rps_α = getfield.(res_α, :Rps)
+Tps_α = getfield.(res_α, :Tps)
+conv_α = Rps_α .+ Tps_α
 println("Panel (b): azimuth sweep at θ=30°")
-println("  Rps(α=0°, optic axis in plane) = ", round(Rps_α[1], digits=8), "  (≈ 0: no conversion)")
-println("  Rps peaks at α≈", round(αdeg[argmax(Rps_α)]), "°  =", round(maximum(Rps_α), digits=5))
+println("  Rps+Tps(α=0°, optic axis in plane) = ", round(conv_α[1], digits=8), "  (≈ 0: no conversion)")
+println("  Rps+Tps peaks at α≈", round(αdeg[argmax(conv_α)]), "°  =", round(maximum(conv_α), digits=5))
 
 # --- Figure ---------------------------------------------------------------
 fig = Figure(size = (980, 430))
 
 axa = Axis(fig[1, 1];
     xlabel = "Incidence angle (°)", ylabel = "energy budget  R + T",
-    title = "A lossless slab that seems to lose ~0.2%")
-limits!(axa, 0, 85, 0.99, 1.002)
+    title = "A lossless slab that seems to lose ~1%")
+limits!(axa, 0, 85, 0.96, 1.005)
 band!(axa, θdeg, naive, complete; color = (:firebrick, 0.18))
 lines!(axa, θdeg, complete; color = :black, linewidth = 2.5,
-    label = "Rpp + Rps + Tpp  (complete) = 1")
+    label = "Rpp + Rps + Tpp + Tps  (complete) = 1")
 lines!(axa, θdeg, naive; color = :dodgerblue3, linewidth = 2.5,
     label = "Rpp + Tpp  (naive)")
-text!(axa, 42, 0.9988; text = "shaded gap = Rps\n(p→s reflection)",
+text!(axa, 42, 0.978; text = "shaded gap = Rps + Tps\n(converted p→s power)",
     align = (:center, :center), fontsize = 11, color = :firebrick4)
 axislegend(axa; position = :rb, framevisible = true)
 
 axb = Axis(fig[1, 2];
-    xlabel = "Optic-axis azimuth α (°)", ylabel = "cross-pol reflectance  Rps",
+    xlabel = "Optic-axis azimuth α (°)", ylabel = "converted power  (p input)",
     title = "Conversion needs an out-of-plane optic axis  (θ = 30°)")
-lines!(axb, αdeg, Rps_α; color = :firebrick3, linewidth = 2.5)
-scatter!(axb, [0], [Rps_α[1]]; color = :black, markersize = 9)
-text!(axb, 3, Rps_α[1]; text = "axis in plane of incidence → Rps = 0",
+lines!(axb, αdeg, conv_α; color = :firebrick3, linewidth = 2.5, label = "Rps + Tps")
+lines!(axb, αdeg, Tps_α; color = :darkorange2, linewidth = 1.8, linestyle = :dash, label = "Tps")
+lines!(axb, αdeg, Rps_α; color = :purple3, linewidth = 1.8, linestyle = :dot, label = "Rps")
+scatter!(axb, [0], [conv_α[1]]; color = :black, markersize = 9)
+text!(axb, 3, conv_α[1]; text = "axis in plane of incidence → no conversion",
     align = (:left, :bottom), fontsize = 11, color = :black)
-limits!(axb, -2, 90, -maximum(Rps_α) * 0.06, maximum(Rps_α) * 1.18)
+limits!(axb, -2, 90, -maximum(conv_α) * 0.06, maximum(conv_α) * 1.18)
+axislegend(axb; position = :rt, framevisible = true)
 
 outpath = joinpath(@__DIR__, "..", "docs", "src", "assets", "examples", "anisotropic_energy_budget.png")
 mkpath(dirname(outpath))