diff --git a/CLAUDE.md b/CLAUDE.md
index 1e99873..d146eaf 100644
--- a/CLAUDE.md
+++ b/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
diff --git a/docs/src/assets/examples/anisotropic_energy_budget.png b/docs/src/assets/examples/anisotropic_energy_budget.png
index be1fd9d..5d3d067 100644
Binary files a/docs/src/assets/examples/anisotropic_energy_budget.png and b/docs/src/assets/examples/anisotropic_energy_budget.png differ
diff --git a/docs/src/examples/anisotropic_energy_budget.md b/docs/src/examples/anisotropic_energy_budget.md
index 8eb42f6..732629c 100644
--- a/docs/src/examples/anisotropic_energy_budget.md
+++ b/docs/src/examples/anisotropic_energy_budget.md
@@ -15,22 +15,23 @@ 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:
@@ -38,30 +39,33 @@ 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).
diff --git a/docs/src/guide/transfer_matrix.md b/docs/src/guide/transfer_matrix.md
index ff43f72..7e1fbf8 100644
--- a/docs/src/guide/transfer_matrix.md
+++ b/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:
diff --git a/docs/src/guide/validation.md b/docs/src/guide/validation.md
index e9a44f0..fb0ddf4 100644
--- a/docs/src/guide/validation.md
+++ b/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
 
diff --git a/examples/anisotropic_energy_budget.jl b/examples/anisotropic_energy_budget.jl
index ea3a00e..642f5f9 100644
--- a/examples/anisotropic_energy_budget.jl
+++ b/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))
diff --git a/examples/anisotropic_fresnel_validation.jl b/examples/anisotropic_fresnel_validation.jl
new file mode 100644
index 0000000..31c83a5
--- /dev/null
+++ b/examples/anisotropic_fresnel_validation.jl
@@ -0,0 +1,205 @@
+# Oracle validation of the per-eigenmode cross-polarized transmittance.
+#
+# transfer() reports four transmittances per input polarization: the
+# co-polarized Tpp, Tss and the converted Tps, Tsp. Each is the Poynting flux
+# of ONE substrate eigenmode — evaluated with that mode's own wavevector —
+# divided by the incident flux. The energy-budget example
+# (anisotropic_energy_budget.jl) shows that these four channels close
+# Rpp + Rps + Tpp + Tps = 1 internally; a budget can close, however, and still
+# split the power wrongly between the channels. Here each channel is checked
+# against independent closed-form textbook physics the Berreman machinery was
+# never fitted to.
+#
+# Validation 1 — per-axis Fresnel oracle (quantitative, < 1e-10).
+#   Air over a semi-infinite uniaxial crystal with the optic axis along ŷ:
+#   Layer(n_o, n_e, n_o, ...). With the plane of incidence x–z, the two
+#   transmitted eigenmodes are PURE p and s waves: the p mode (E in x–z) sees
+#   only ε_xx = ε_zz = n_o², the s mode (E ∥ ŷ) only ε_yy = n_e². So Tpp must
+#   equal the isotropic single-interface p-polarization Fresnel transmittance
+#   with index n_o, and Tss the s-polarization one with index n_e — including
+#   the flux factor n·cosθ_t/cosθ_i, which DIFFERS between the two output
+#   modes because each refracts with its own index. Cross terms must vanish.
+#
+# Validation 2 — half-wave plate oracle (quantitative, < 1e-10).
+#   A calcite-like slab (n_o = 1.658, n_e = 1.486) with the optic axis in the
+#   layer plane at 45° to the plane of incidence (euler = (π/4, π/2, 0)) and
+#   thickness t = λ₀/(2Δn) at λ₀ = 0.59 μm: a half-wave plate. At normal
+#   incidence a p input splits equally onto the two eigenaxes; the slab
+#   retards one by π and the recombined polarization is rotated by 90°, so the
+#   transmitted power exits almost entirely s-polarized (Tps). Two formulas:
+#     (a) the classic waveplate conversion T_conv = sin²(2·45°)·sin²(πΔn·t/λ),
+#         scaled by the single-pass Fresnel surface losses (fringe-free);
+#     (b) an EXACT Jones × Airy oracle — at normal incidence each eigenmode is
+#         an independent Fabry–Pérot etalon with its own index, so
+#         Tps = |t_e − t_o|²/4 and Tpp = |t_e + t_o|²/4 in closed form.
+#   (b) must agree with transfer() to < 1e-10 across the sweep; (a) is the
+#   envelope the etalon fringes oscillate around.
+#
+# Validation 3 — the reading the fix changed (for readers migrating).
+#   In releases ≤ v3.2.1, Tpp was the TOTAL per-input transmitted flux and
+#   Tps = |t_ps|² a raw amplitude ratio: at the half-wave point
+#   Rpp + Rps + Tpp = 1 looked like a closed budget (with Tpp ≈ 0.848), while
+#   the four-channel sum double-counted the converted power to ≈ 1.85. Now
+#   Tpp is the co-polarized SHARE (≈ 0.0004 here) and the budget that closes
+#   is the four-channel one.
+#
+# Run from the examples environment. The examples project normally uses the
+# registered TransferMatrix, so point it at this checkout first:
+#   julia --project=examples -e 'using Pkg; Pkg.develop(path="."); Pkg.instantiate()'
+#   julia --project=examples examples/anisotropic_fresnel_validation.jl
+
+using TransferMatrix
+using CairoMakie
+
+function check_below(label, deviation, tol)
+    ok = deviation < tol
+    println("  ", label, " = ", deviation, ok ? "  (PASS < $tol)" : "  (FAIL, tol $tol)")
+    ok || error("validation failed: " * label)
+end
+
+no, ne = 1.658, 1.486
+λ0 = 0.59
+air = Layer(λ -> 1.0, 1.0)
+
+# --- Validation 1: air | uniaxial half-space, optic axis along ŷ ------------
+# Textbook single-interface Fresnel transmittance (incident index n₁ = 1),
+# written out for a substrate of index n:
+#   cosθ_t = √(1 − sin²θᵢ/n²)                       (Snell)
+#   t_p = 2cosθᵢ/(n·cosθᵢ + cosθ_t)                 (amplitude)
+#   t_s = 2cosθᵢ/(cosθᵢ + n·cosθ_t)
+#   T   = (n·cosθ_t/cosθᵢ)·|t|²                     (transmitted flux ratio)
+# The p mode refracts with n_o and the s mode with n_e, so each channel
+# carries its OWN n·cosθ_t flux factor.
+
+substrate_y = Layer(λ -> no, λ -> ne, λ -> no, 1.0)
+θdeg = 0.0:1.0:80.0
+θi = deg2rad.(θdeg)
+res1 = [transfer(λ0, [air, substrate_y]; θ = θ) for θ in θi]
+Tpp1 = getfield.(res1, :Tpp)
+Tss1 = getfield.(res1, :Tss)
+
+cosθp = sqrt.(1 .- (sin.(θi) ./ no) .^ 2)
+cosθs = sqrt.(1 .- (sin.(θi) ./ ne) .^ 2)
+tp = 2 .* cos.(θi) ./ (no .* cos.(θi) .+ cosθp)
+ts = 2 .* cos.(θi) ./ (cos.(θi) .+ ne .* cosθs)
+Tp_fresnel = no .* cosθp ./ cos.(θi) .* abs2.(tp)
+Ts_fresnel = ne .* cosθs ./ cos.(θi) .* abs2.(ts)
+
+cross1 = maximum(res -> max(res.Tps, res.Tsp, res.Rps, res.Rsp), res1)
+
+println("Validation 1: air | uniaxial half-space (optic axis ∥ ŷ), θ = 0–80°")
+check_below("max |Tpp − Fresnel T_p(n_o)|", maximum(abs.(Tpp1 .- Tp_fresnel)), 1e-10)
+check_below("max |Tss − Fresnel T_s(n_e)|", maximum(abs.(Tss1 .- Ts_fresnel)), 1e-10)
+check_below("max cross channel (Tps, Tsp, Rps, Rsp)", cross1, 1e-10)
+
+# --- Validation 2: calcite-like half-wave plate at 45° ----------------------
+Δn = no - ne
+t_hwp = λ0 / (2Δn)
+hwp = Layer(λ -> no, λ -> no, λ -> ne, t_hwp; euler = (π/4, π/2, 0.0))
+λs = 0.40:0.002:1.00
+res2 = [transfer(λ, [air, hwp, air]) for λ in λs]
+Tps2 = getfield.(res2, :Tps)
+Tpp2 = getfield.(res2, :Tpp)
+
+# (a) Classic waveplate conversion, scaled by surface losses.
+# a_j = t_in·t_out = [2/(1+n_j)]·[2n_j/(1+n_j)] is the two-face single-pass
+# amplitude transmission of eigenmode j; the fringe-free envelope of the
+# converted power is a_o·a_e·sin²(2β)·sin²(πΔn·t/λ) with β = 45°.
+a_o = 4no / (1 + no)^2
+a_e = 4ne / (1 + ne)^2
+Tps_classic = a_o * a_e * sin(2 * (π/4))^2 .* sin.(π * Δn * t_hwp ./ λs) .^ 2
+
+# (b) Exact Jones × Airy oracle. Each eigenmode traverses the slab as its own
+# Fabry–Pérot etalon:
+#   t_slab(n) = t_in·t_out·e^{iδ}/(1 − r²e^{2iδ}),  δ = 2πnt/λ,  r = (n−1)/(n+1)
+# With the eigenaxes at 45°, an x̂ (p) input recombines in the output air as
+# E_x = (t_e + t_o)/2 and E_y = (t_e − t_o)/2, so
+#   Tpp = |t_e + t_o|²/4,  Tps = |t_e − t_o|²/4.
+function slab_transmission(n, t, λ)
+    δ = 2π * n * t / λ
+    r = (n - 1) / (n + 1)
+    return (2 / (1 + n)) * (2n / (1 + n)) * cis(δ) / (1 - r^2 * cis(2δ))
+end
+t_o = slab_transmission.(no, t_hwp, λs)
+t_e = slab_transmission.(ne, t_hwp, λs)
+Tps_airy = abs2.(t_e .- t_o) ./ 4
+Tpp_airy = abs2.(t_e .+ t_o) ./ 4
+
+budget_p = [r.Rpp + r.Rps + r.Tpp + r.Tps for r in res2]
+budget_s = [r.Rss + r.Rsp + r.Tss + r.Tsp for r in res2]
+
+println("Validation 2: half-wave plate (t = ", round(t_hwp, digits = 5),
+        " μm), λ = 0.40–1.00 μm, normal incidence")
+check_below("max |Tps − Jones×Airy oracle|", maximum(abs.(Tps2 .- Tps_airy)), 1e-10)
+check_below("max |Tpp − Jones×Airy oracle|", maximum(abs.(Tpp2 .- Tpp_airy)), 1e-10)
+check_below("max |p-input budget − 1|", maximum(abs.(budget_p .- 1)), 1e-12)
+check_below("max |s-input budget − 1|", maximum(abs.(budget_s .- 1)), 1e-12)
+
+# The fringe-free conversion envelope peaks exactly at the design wavelength;
+# the exact curve's global maximum rides the etalon fringe nearest to it
+# (fringe spacing λ²/(2n̄t) ≈ 0.065 μm here).
+i0 = argmin(abs.(λs .- λ0))
+λpeak_env = λs[argmax(Tps_classic)]
+λpeak = λs[argmax(Tps2)]
+println("  conversion envelope peaks at λ = ", λpeak_env, " μm  (design λ₀ = ", λ0, " μm)")
+println("  exact Tps max at λ = ", λpeak, " μm (nearest etalon fringe);  Tps(λ₀) = ",
+        round(Tps2[i0], digits = 6), ",  Tpp(λ₀) = ", round(Tpp2[i0], digits = 6))
+abs(λpeak_env - λ0) < 1e-9 || error("validation failed: envelope does not peak at the design wavelength")
+abs(λpeak - λ0) < 0.05 || error("validation failed: Tps peak is not within one fringe of the design wavelength")
+Tps2[i0] > 0.8 || error("validation failed: conversion at the design wavelength should be nearly complete")
+Tpp2[i0] < 1e-3 || error("validation failed: Tpp at the half-wave point should be ≈ 0")
+
+# --- Validation 3: the reading the per-mode fix changed ---------------------
+# Migration note (releases ≤ v3.2.1 → this version): Tpp used to be the TOTAL
+# transmitted Poynting flux for a p input (both output polarizations, one
+# shared wavevector) and Tps = |t_ps|² a raw amplitude ratio. Rpp + Rps + Tpp
+# was therefore the sum that read 1.0, and adding Tps double-counted the
+# converted power (≈ 1.85 for this stack). With the per-eigenmode
+# decomposition, Tpp at the half-wave point drops from ≈ 0.848 (the total) to
+# ≈ 0.0004 (the co-polarized share), Rpp + Rps + Tpp is correctly ≪ 1, and
+# the budget that closes is the four-channel one.
+res_hw = transfer(λ0, [air, hwp, air])
+old_reading = res_hw.Rpp + res_hw.Rps + res_hw.Tpp
+println("Validation 3: half-wave point, per-mode decomposition vs the pre-fix reading")
+println("  Tpp = ", round(res_hw.Tpp, digits = 6),
+        "   (co-polarized share; ≤ v3.2.1 put the TOTAL ≈ 0.848 here)")
+println("  Tps = ", round(res_hw.Tps, digits = 6),
+        "   (converted share; ≤ v3.2.1 reported |t_ps|² instead)")
+println("  Rpp + Rps + Tpp       = ", round(old_reading, digits = 6),
+        "   (the old \"closed\" budget — correctly ≠ 1 now)")
+println("  Rpp + Rps + Tpp + Tps = ", res_hw.Rpp + res_hw.Rps + res_hw.Tpp + res_hw.Tps,
+        "   (the budget that closes)")
+old_reading < 0.9 || error("validation failed: Tpp looks like a total flux again (pre-fix semantics?)")
+
+# --- Figure ------------------------------------------------------------------
+fig = Figure(size = (980, 430))
+
+axa = Axis(fig[1, 1];
+    xlabel = "Incidence angle (°)", ylabel = "Transmittance",
+    title = "Air | uniaxial half-space: per-axis Fresnel oracle")
+lines!(axa, θdeg, Tp_fresnel; color = :black, linewidth = 2.5,
+    label = "Fresnel T_p(n_o) analytic")
+lines!(axa, θdeg, Ts_fresnel; color = :gray55, linewidth = 2.5,
+    label = "Fresnel T_s(n_e) analytic")
+scatter!(axa, θdeg[1:4:end], Tpp1[1:4:end]; color = :orangered3, markersize = 8,
+    label = "transfer() Tpp")
+scatter!(axa, θdeg[1:4:end], Tss1[1:4:end]; color = :dodgerblue3, marker = :utriangle,
+    markersize = 8, label = "transfer() Tss")
+axislegend(axa; position = :lb, framevisible = true)
+
+axb = Axis(fig[1, 2];
+    xlabel = "Wavelength (μm)", ylabel = "Transmittance (p input)",
+    title = "Half-wave plate at 45°: waveplate oracle")
+limits!(axb, 0.38, 1.02, -0.04, 1.38)
+lines!(axb, λs, Tps_classic; color = :gray55, linewidth = 2, linestyle = :dash,
+    label = "classic sin²(πΔn·t/λ) × surface loss")
+lines!(axb, λs, Tps2; color = :firebrick3, linewidth = 2.5, label = "transfer() Tps")
+lines!(axb, λs, Tpp2; color = :dodgerblue3, linewidth = 1.8, label = "transfer() Tpp")
+scatter!(axb, λs[1:15:end], Tps_airy[1:15:end]; color = :black, markersize = 6,
+    label = "Jones × Airy exact")
+vlines!(axb, λ0; color = :black, linestyle = :dot)
+text!(axb, λ0 + 0.015, 0.06; text = "design λ₀", align = (:left, :center), fontsize = 11)
+axislegend(axb; position = :rt, framevisible = true)
+
+save(joinpath(@__DIR__, "anisotropic_fresnel_validation.png"), fig)
+println("saved ", joinpath(@__DIR__, "anisotropic_fresnel_validation.png"))
diff --git a/src/coefficients.jl b/src/coefficients.jl
index 0f35ddd..23d4d6f 100644
--- a/src/coefficients.jl
+++ b/src/coefficients.jl
@@ -5,8 +5,9 @@
 const _C_CIRC = SMatrix{2,2,ComplexF64}(1, im, 1, -im) / sqrt(2)      # [1 1; im -im]
 const _C_CIRC_INV = SMatrix{2,2,ComplexF64}(1, 1, -im, im) / sqrt(2)  # [1 -im; 1 im]
 
-# Transmission-null tolerance: below this the transmitted Poynting flux is taken
-# as zero, so circular transmittance is defined as 0 (avoids 0/0 at a perfect mirror).
+# Transmission-null tolerance: below this the transmitted amplitude weight
+# |t_c[L,j]|² + |t_c[R,j]|² is taken as zero, so circular transmittance is
+# defined as 0 (avoids 0/0 at a perfect mirror).
 const _CIRC_T_TOL = 1e-12
 
 # Transform a 2×2 linear Jones matrix ([out,in] ordered p,s) into the circular
@@ -14,19 +15,49 @@ const _CIRC_T_TOL = 1e-12
 _jones_to_circular(J_lin) = _C_CIRC_INV * J_lin * _C_CIRC
 
 # Assemble circular-basis R/T from the linear complex Jones coefficients
-# r = (rpp,rps,rss,rsp), t = (tpp,tps,tsp,tss) and the Poynting diagonal
-# transmittances Tpp, Tss. Reflectance is |r_circ|^2; transmittance is N·|t_circ|^2
-# with the single polarization-independent Poynting scalar N (exact for an isotropic
-# substrate; see `transfer` docstring for caveats).
-function _circular_result(r, t, Tpp, Tss)
-    r_lin = SMatrix{2,2,ComplexF64}(r[1], r[4], r[2], r[3])  # [rpp rps; rsp rss]
-    t_lin = SMatrix{2,2,ComplexF64}(t[1], t[3], t[2], t[4])  # [tpp tps; tsp tss]
+# r = (rpp,rps,rss,rsp), t = (tpp,tps,tsp,tss) and the Poynting flux factors:
+# f_out1/f_out2 are the unit-amplitude z-fluxes of the two substrate eigenmodes
+# (p-like / s-like), f_in_p/f_in_s those of the incident modes.
+#
+# Reflectance is |r_circ|^2. Transmittance for circular input j is computed by
+# mapping the input Jones vector through t_lin to the substrate (p,s) amplitude
+# pair `a`, whose per-mode flux `f_out1|a₁|² + f_out2|a₂|²` is the EXACT total
+# transmitted power (distinct lossless eigenmodes carry no cross z-flux). That
+# total is then split between the L/R outputs in proportion to |t_circ|². For an
+# isotropic substrate (f_out1 = f_out2) the split itself is exact, because the
+# circular unit vectors are orthonormal combinations of degenerate modes; for an
+# anisotropic substrate the two eigenmodes are not circular, so only the per-input
+# TOTAL is exact and the L/R split is the natural |t_circ|²-weighted estimate
+# (see `transfer` docstring).
+function _circular_result(r, t, f_out1, f_out2, f_in_p, f_in_s)
+    # Package coefficients are indexed in-first (rps = p input → s output), so
+    # the [out,in] Jones matrices place the p→s conversion in row 2, column 1.
+    # (Using the transposed orientation is invisible for isotropic stacks —
+    # their Jones matrices are diagonal — but breaks flux unitarity, and hence
+    # the circular energy budget, for polarization-converting stacks.)
+    r_lin = SMatrix{2,2,ComplexF64}(r[1], r[2], r[4], r[3])  # [rpp rsp; rps rss]
+    t_lin = SMatrix{2,2,ComplexF64}(t[1], t[2], t[3], t[4])  # [tpp tsp; tps tss]
     r_c = _jones_to_circular(r_lin)
     t_c = _jones_to_circular(t_lin)
     Rc = abs2.(r_c)
-    denom = abs2(t[1]) + abs2(t[4])  # |tpp|^2 + |tss|^2
-    N = denom > _CIRC_T_TOL ? (Tpp + Tss) / denom : zero(denom)
-    Tc = N .* abs2.(t_c)
+
+    # Unit-amplitude circular incident flux: |χ_p|² f_in_p + |χ_s|² f_in_s with
+    # χ = (1, ±i)/√2 (incident p/s modes carry no cross z-flux).
+    f_in = (f_in_p + f_in_s) / 2
+
+    Tc = MMatrix{2,2,Float64}(undef)
+    for j in 1:2
+        a = t_lin * _C_CIRC[:, j]                 # substrate (p,s) amplitudes for input j (1=L, 2=R)
+        w = abs2(t_c[1, j]) + abs2(t_c[2, j])     # = |a₁|² + |a₂|² (C is unitary)
+        if w > _CIRC_T_TOL
+            T_tot = (f_out1 * abs2(a[1]) + f_out2 * abs2(a[2])) / f_in
+            Tc[1, j] = T_tot * abs2(t_c[1, j]) / w
+            Tc[2, j] = T_tot * abs2(t_c[2, j]) / w
+        else
+            Tc[1, j] = 0.0
+            Tc[2, j] = 0.0
+        end
+    end
     # circular matrices are [out,in] ordered (L,R): [1,1]=LL [1,2]=LR [2,1]=RL [2,2]=RR
     return CircularTransferResult(Tc[2, 2], Tc[1, 1], Tc[2, 1], Tc[1, 2],
                                   Rc[2, 2], Rc[1, 1], Rc[2, 1], Rc[1, 2])
@@ -45,9 +76,10 @@ reflection and transmission Jones blocks can be read straight off its elements.
 The function returns the co- and cross-polarized amplitude coefficients
 (`rpp, rss, rps, rsp`, `tpp, tss, tps, tsp`) and their squared magnitudes.
 
-Reflectance is ``R = |r|^2`` directly. The returned ``|t|^2`` is *not* the
-physical transmittance — the transmitted wave is in a different medium — so the
-diagonal `Tpp`/`Tss` reported by `transfer` come from the Poynting-vector
+Reflectance is ``R = |r|^2`` directly. The returned ``|t|^2`` values (all four,
+co- and cross-polarized) are *not* physical transmittances — the transmitted
+wave is in a different medium — so every transmittance reported by `transfer`
+(`Tpp`, `Tss`, `Tps`, `Tsp`) comes from the per-mode Poynting-vector
 calculation (`calculate_tr(S::Poynting)`) instead.
 
 The ``2\\times2`` transfer-matrix relations originate with Yeh (1979); the
diff --git a/src/matrix_constructors.jl b/src/matrix_constructors.jl
index 1f88f1b..2ff700e 100644
--- a/src/matrix_constructors.jl
+++ b/src/matrix_constructors.jl
@@ -531,7 +531,10 @@ E-field eigenvectors for a layer with a full 3×3 permeability tensor `μ`, as t
 null space of `W(q) = K μ⁻¹ K + ε` for each mode `q`. Degenerate pairs (qᵢ ≈ qⱼ,
 e.g. isotropic media) share a 2-D null space; both orthogonal null vectors are
 assigned, avoiding the singular dynamical matrix a naive per-mode null space
-produces. Rows are unit-Hermitian-normalized, matching [`calculate_E_modes`].
+produces. Degenerate pairs are oriented to the fixed p/s polarization references
+of the scalar path (see [`_orient_degenerate_pair`](@ref)) so the r/t coefficient
+labels and per-mode transmittance channels keep their p/s meaning. Rows are
+unit-Hermitian-normalized, matching [`calculate_E_modes`].
 """
 function calculate_E_modes_tensor(k_par, q, ε, μ; rtol=1e-7)
     μinv = inv(SMatrix{3,3,ComplexF64}(μ)); εm = SMatrix{3,3,ComplexF64}(ε)
@@ -545,15 +548,68 @@ function calculate_E_modes_tensor(k_par, q, ε, μ; rtol=1e-7)
                 partner = n; break
             end
         end
-        v1 = SVector{3,ComplexF64}(F.V[:,3]...); E_modes[m,:] = v1 ./ norm(v1); assigned[m] = true
-        if partner != 0
-            v2 = SVector{3,ComplexF64}(F.V[:,2]...); E_modes[partner,:] = v2 ./ norm(v2); assigned[partner] = true
+        if partner == 0
+            v1 = SVector{3,ComplexF64}(F.V[:,3]...)
+            E_modes[m,:] = v1 ./ norm(v1)
+            assigned[m] = true
+        else
+            v1 = SVector{3,ComplexF64}(F.V[:,3]...)
+            v2 = SVector{3,ComplexF64}(F.V[:,2]...)
+            w_p, w_s = _orient_degenerate_pair(v1, v2, m)
+            E_modes[m,:] = w_p
+            E_modes[partner,:] = w_s
+            assigned[m] = true
+            assigned[partner] = true
         end
     end
     return SMatrix(E_modes)
 end
 
 
+"""
+    _orient_degenerate_pair(v1, v2, slot)
+
+Rotate an orthonormal basis `(v1, v2)` of a degenerate 2-D eigenmode space so
+the pair matches the fixed p/s polarization references of the scalar-μ path
+(`calculate_E_modes`): the first returned vector is p-like (no Ey component;
+Ex real-positive for the forward pair `slot ≤ 2`, real-negative for the
+backward pair) and the second is s-like (Ey real-positive).
+
+Any orthonormal basis of the degenerate space is a valid eigenvector pair, but
+the transfer-matrix coefficient labels (`tpp`, `tps`, ...) and the per-mode
+transmittance channels assume this orientation; an arbitrary (e.g. raw SVD)
+basis silently relabels or mixes the p/s channels of a magnetic ambient or
+substrate. If the space has no ŷ component the input basis is returned as-is.
+"""
+function _orient_degenerate_pair(v1, v2, slot)
+    # Project ŷ onto span{v1, v2} to extract the s-like direction.
+    w_s = v1 * conj(v1[2]) + v2 * conj(v2[2])
+    n_s = norm(w_s)
+    if n_s > 1e-9
+        w_s = w_s / n_s
+        # p-like: the orthogonal complement of w_s within the span.
+        w_p = v1 - (w_s' * v1) * w_s
+        if norm(w_p) < 1e-9
+            w_p = v2 - (w_s' * v2) * w_s
+        end
+        w_p = w_p / norm(w_p)
+    else
+        w_p, w_s = v1, v2
+    end
+    # Fix the arbitrary complex phases to the scalar-path sign conventions:
+    # s mode Ey real-positive; p mode Ex real-positive (forward) or
+    # real-negative (backward), matching E_modes[3,1] = -1 in calculate_E_modes.
+    if abs(w_s[2]) > 1e-12
+        w_s = w_s * (conj(w_s[2]) / abs(w_s[2]))
+    end
+    if abs(w_p[1]) > 1e-12
+        phase = conj(w_p[1]) / abs(w_p[1])
+        w_p = slot <= 2 ? w_p * phase : w_p * (-phase)
+    end
+    return w_p, w_s
+end
+
+
 """
     dynamical_matrix(k_par, q, E_modes, μ::AbstractMatrix)
 
diff --git a/src/poynting.jl b/src/poynting.jl
index ffd3ea8..1b3d309 100644
--- a/src/poynting.jl
+++ b/src/poynting.jl
@@ -115,17 +115,27 @@ function poynting(k_par, q_in, q_out, E_modes_in, E_modes_out, t_coefs, r_coefs,
     k_out ./= c_0
     k_out = SMatrix(k_out)
 
-    # Transmitted power is the sum of the two substrate eigenmodes' Poynting
-    # vectors, EACH evaluated with its OWN wavevector (k_out modes 1, 2). The
-    # two modes carry different wavevectors when the substrate is anisotropic,
-    # and the cross-interference between distinct propagating modes transports
-    # no net z-flux, so the per-mode sum is the physical transmitted flux.
-    # (For an isotropic substrate k_out[1]=k_out[2] and this reduces to the
-    # single-wavevector form, since the p/s cross term carries no z-power.)
-    S_out_p = real(0.5 * E_out_p1 × conj(μout_inv * (k_out[1, :] × E_out_p1))) +
-              real(0.5 * E_out_p2 × conj(μout_inv * (k_out[2, :] × E_out_p2)))
-    S_out_s = real(0.5 * E_out_s1 × conj(μout_inv * (k_out[1, :] × E_out_s1))) +
-              real(0.5 * E_out_s2 × conj(μout_inv * (k_out[2, :] × E_out_s2)))
+    # Transmitted power is resolved PER substrate eigenmode, EACH evaluated
+    # with its OWN wavevector (k_out modes 1, 2). The two modes carry different
+    # wavevectors when the substrate is anisotropic, and the cross-interference
+    # between distinct eigenmodes of a lossless substrate transports no net
+    # z-flux, so each per-mode flux is the physical power in that output
+    # polarization channel and their sum is the total transmitted flux.
+    # (For an isotropic substrate k_out[1]=k_out[2] and the cross term vanishes
+    # by p/s mode orthogonality.)
+    S_out_pp = real(0.5 * E_out_p1 × conj(μout_inv * (k_out[1, :] × E_out_p1)))
+    S_out_ps = real(0.5 * E_out_p2 × conj(μout_inv * (k_out[2, :] × E_out_p2)))
+    S_out_sp = real(0.5 * E_out_s1 × conj(μout_inv * (k_out[1, :] × E_out_s1)))
+    S_out_ss = real(0.5 * E_out_s2 × conj(μout_inv * (k_out[2, :] × E_out_s2)))
+    S_out_p = S_out_pp + S_out_ps
+    S_out_s = S_out_sp + S_out_ss
+
+    # Unit-amplitude z-flux of each substrate eigenmode (flux-per-|amplitude|²
+    # conversion factors, needed by the circular-basis transmittance).
+    E_out_1 = E_modes_out[1, :]
+    E_out_2 = E_modes_out[2, :]
+    f_out1 = real(0.5 * E_out_1 × conj(μout_inv * (k_out[1, :] × E_out_1)))[3]
+    f_out2 = real(0.5 * E_out_2 × conj(μout_inv * (k_out[2, :] × E_out_2)))[3]
 
     # Reflected waves: use incident-medium wavevectors (k_in modes 3,4),
     # not substrate wavevectors, because reflected light propagates backward
@@ -133,7 +143,8 @@ function poynting(k_par, q_in, q_out, E_modes_in, E_modes_out, t_coefs, r_coefs,
     S_refl_p = real(0.5 * E_ref_p × conj(μin_inv * (k_in[3, :] × E_ref_p)))
     S_refl_s = real(0.5 * E_ref_s × conj(μin_inv * (k_in[4, :] × E_ref_s)))
 
-    return Poynting(S_out_p, S_in_p, S_out_s, S_in_s, S_refl_p, S_refl_s)
+    return Poynting(S_out_p, S_in_p, S_out_s, S_in_s, S_refl_p, S_refl_s,
+                    S_out_pp, S_out_ps, S_out_sp, S_out_ss, f_out1, f_out2)
 end
 
 
@@ -225,7 +236,19 @@ Calculate transmittance and reflectance from the Poynting vector struct,
 which contains incident, transmitted, and reflected energy flux for both
 p-polarized and s-polarized waves.
 
-Returns `(Tpp, Tss, Rpp, Rss)`.
+Returns `(Tpp, Tss, Rpp, Rss, Tps, Tsp)`.
+
+The transmittances are per-output-mode flux ratios: `Tpp` is the power carried
+by the p-like substrate eigenmode alone (co-polarized channel) and `Tps` the
+power carried by the s-like eigenmode (cross-polarized channel), each divided
+by the incident p flux; likewise `Tss`/`Tsp` for s input. For a lossless stack
+the per-input budgets close exactly: `Rpp + Rps + Tpp + Tps = 1`.
+
+`Rpp`/`Rss` here are the TOTAL reflected flux per input polarization (the
+reflected field is a coherent sum of both backward modes, so for converting
+stacks this equals `|rpp|² + |rps|²`, not the co-polarized piece alone).
+[`transfer`](@ref) takes its per-channel reflectances from ``|r|^2`` instead —
+see the `calculate_tr(M_sys)` method.
 
 # Sign Convention
 The reflected Poynting vector z-component is negative (pointing in -z direction),
@@ -233,12 +256,14 @@ so the negative sign in `Rpp = -S.refl_p[3] / S.in_p[3]` yields positive reflect
 """
 function calculate_tr(S::Poynting)
 
-    Tpp = S.out_p[3] / S.in_p[3]
-    Tss = S.out_s[3] / S.in_s[3]
+    Tpp = S.out_pp[3] / S.in_p[3]
+    Tss = S.out_ss[3] / S.in_s[3]
+    Tps = S.out_ps[3] / S.in_p[3]
+    Tsp = S.out_sp[3] / S.in_s[3]
 
     # Reflected Poynting vector points in -z, so negate to get positive R
     Rpp = -S.refl_p[3] / S.in_p[3]
     Rss = -S.refl_s[3] / S.in_s[3]
 
-    return Tpp, Tss, Rpp, Rss
+    return Tpp, Tss, Rpp, Rss, Tps, Tsp
 end
diff --git a/src/results.jl b/src/results.jl
index 767e4c7..94783ba 100644
--- a/src/results.jl
+++ b/src/results.jl
@@ -1,3 +1,27 @@
+"""
+    Poynting
+
+Time-averaged Poynting vectors for the incident, transmitted, and reflected
+waves of a stack, with the transmitted flux resolved per substrate eigenmode.
+
+The transmitted field for each input polarization is a superposition of the two
+forward substrate eigenmodes (mode 1 is p-like, mode 2 is s-like). Each
+per-channel vector is that mode's own Poynting vector, evaluated with its own
+wavevector:
+
+- `out_pp` / `out_ps`: p input → substrate mode 1 / mode 2
+- `out_sp` / `out_ss`: s input → substrate mode 1 / mode 2
+
+`out_p = out_pp + out_ps` and `out_s = out_sp + out_ss` are the transmitted
+totals per input polarization. Distinct eigenmodes of a lossless substrate
+carry no cross z-flux (their interference term either vanishes by mode
+orthogonality or oscillates in z and must vanish because the total flux is
+z-independent), so the per-mode split is exact wherever a transmittance into
+the substrate is well defined.
+
+`f_out1` / `f_out2` are the z-fluxes of the two substrate eigenmodes at unit
+amplitude, used to convert circular-basis ``|t|^2`` into flux ratios.
+"""
 struct Poynting
     out_p::SVector{3, Float64}
     in_p::SVector{3, Float64}
@@ -5,10 +29,12 @@ struct Poynting
     in_s::SVector{3, Float64}
     refl_p::SVector{3, Float64}
     refl_s::SVector{3, Float64}
-
-    function Poynting(out_p::T, in_p::T, out_s::T, in_s::T, refl_p::T, refl_s::T) where {T<:SVector{3, Float64}}
-        new(out_p, in_p, out_s, in_s, refl_p, refl_s)
-    end
+    out_pp::SVector{3, Float64}
+    out_ps::SVector{3, Float64}
+    out_sp::SVector{3, Float64}
+    out_ss::SVector{3, Float64}
+    f_out1::Float64
+    f_out2::Float64
 end
 
 """
@@ -82,6 +108,11 @@ For sweep calculations via `sweep_angle()` or `sweep_thickness()`, `T` is `Matri
     The cross-polarization terms (`Tps`, `Tsp`, `Rps`, `Rsp`) are zero for isotropic media
     and become non-zero for anisotropic (birefringent) materials.
 
+    Every entry is a true energy flux ratio: transmittances are per-output-mode
+    Poynting fluxes (see [`transfer`](@ref)) and reflectances are ``|r|^2``.
+    For a lossless stack the per-input-polarization budgets close exactly:
+    `Rpp + Rps + Tpp + Tps = 1` and `Rss + Rsp + Tss + Tsp = 1`.
+
 # Examples
 ```julia
 # Single wavelength - returns TransferResult{Float64}
diff --git a/src/transfer.jl b/src/transfer.jl
index 1fac701..d6bafca 100644
--- a/src/transfer.jl
+++ b/src/transfer.jl
@@ -23,6 +23,17 @@ rather than ``|t|^2``, because the transmitted wave propagates in a different
 medium than the incident wave. As noted in the 2019 erratum (JOSAB 36, 3246):
 ``T ≠ |t|^2`` in general; only when the substrate is vacuum does ``T = |t|^2``.
 
+All four linear transmittances are per-output-mode flux ratios: the transmitted
+field is decomposed into the two forward substrate eigenmodes (p-like and
+s-like), and each channel is that mode's own Poynting z-flux — with its own
+wavevector — divided by the incident flux. Distinct eigenmodes of a lossless
+substrate carry no cross z-flux, so this decomposition is exact for any
+transparent exit medium, isotropic or birefringent, at any incidence angle.
+The per-input-polarization energy budgets therefore close exactly for a
+lossless stack:
+
+``R_{pp} + R_{ps} + T_{pp} + T_{ps} = 1``, and likewise for s input.
+
 # Arguments
 - `λ`: Wavelength in μm (must match units used for layer thicknesses)
 - `layers`: Vector of `Layer` objects representing the stack
@@ -59,17 +70,19 @@ Higham, 2005, https://doi.org/10.1137/04061101X
 
   Circular **reflectance** is `|r_circ|²` (a true energy ratio under the same
   condition as linear `R`: transparent incident medium, cf. issue #72). Circular
-  **transmittance** is `N·|t_circ|²` with a single Poynting normalization scalar
-  `N = (Tpp+Tss)/(|tpp|²+|tss|²)`; this is exact and energy-conserving for an
-  isotropic substrate (and reduces to `|t_circ|²` for a vacuum/index-matched
-  substrate). For an **anisotropic substrate** the single scalar `N` is only
-  approximate: the two transmitted eigenmodes carry different wavevectors, so
-  their Poynting-to-`|t|²` conversion factors differ, and one polarization-
-  independent scalar cannot capture both. (The linear basis avoids this by
-  summing each transmitted eigenmode with its own wavevector in `poynting()`,
-  so linear `T` — including the cross-pol `Tps/Tsp` — is per-mode exact.)
-  Circular `T` for an anisotropic substrate is therefore not guaranteed to be
-  energy-conserving to machine precision. `validate` applies to the linear basis only.
+  **transmittance** maps each circular input through the linear Jones matrix to
+  the substrate (p,s) amplitude pair, whose per-eigenmode Poynting flux gives
+  the exact total transmitted power for that input; the total is then split
+  between the L/R outputs in proportion to `|t_circ|²`. For an **isotropic
+  substrate** (degenerate eigenmodes) the split itself is exact, so circular
+  `T` is energy-conserving to machine precision (and reduces to `|t_circ|²`
+  for a vacuum/index-matched substrate). For an **anisotropic substrate** the
+  two transmitted eigenmodes carry different wavevectors and are not circular,
+  so the L/R *split* is only the natural `|t_circ|²`-weighted estimate — but
+  the per-input total `Trr+Tlr` (resp. `Tll+Trl`) remains exact, so the energy
+  budget still closes. (The linear basis is per-mode exact in all channels,
+  because the substrate eigenmodes ARE the linear output channels.)
+  `validate` applies to the linear basis only.
 
 # Wave Propagation Convention
 - Light propagates in the **+z direction** (from first layer toward last layer)
@@ -103,19 +116,20 @@ end
 
 function _assemble(::Val{:linear}, M_sys, S, λ, layers, sd, validate)
     r, R, t, T = calculate_tr(M_sys)
-    Tpp, Tss, _, _ = calculate_tr(S)
 
-    # Diagonal terms from matrix calculation; cross terms remapped from the
-    # SVector packing (R = (Rpp,Rss,Rsp,Rps), T = (Tpp,Tss,Tsp,Tps)).
+    # ALL four transmittances are per-output-mode Poynting flux ratios (the
+    # |t|² values from calculate_tr(M_sys) are amplitude ratios, not fluxes).
+    # Reflectances are |r|² per channel, exact for a transparent isotropic
+    # incident medium. R SVector packing: R = (Rpp,Rss,Rsp,Rps).
+    Tpp, Tss, _, _, Tps, Tsp = calculate_tr(S)
     Rpp = R[1]
     Rss = R[2]
     Rps = R[4]
     Rsp = R[3]
-    Tps = T[4]
-    Tsp = T[3]
 
     if validate
-        _validate_physics(λ, layers, Tpp, Tss, Rpp, Rss; sheets=sd)
+        _validate_physics(λ, layers, Tpp, Tss, Rpp, Rss;
+                          sheets=sd, Tps=Tps, Tsp=Tsp, Rps=Rps, Rsp=Rsp)
     end
 
     return TransferResult(Tpp, Tss, Tps, Tsp, Rpp, Rss, Rps, Rsp)
@@ -123,8 +137,7 @@ end
 
 function _assemble(::Val{:circular}, M_sys, S, λ, layers, sd, validate)
     r, _, t, _ = calculate_tr(M_sys)
-    Tpp, Tss, _, _ = calculate_tr(S)
-    return _circular_result(r, t, Tpp, Tss)
+    return _circular_result(r, t, S.f_out1, S.f_out2, S.in_p[3], S.in_s[3])
 end
 
 _assemble(::Val{B}, M_sys, S, λ, layers, sd, validate) where {B} =
@@ -132,17 +145,23 @@ _assemble(::Val{B}, M_sys, S, λ, layers, sd, validate) where {B} =
 
 
 """
-    _validate_physics(λ, layers, Tpp, Tss, Rpp, Rss; atol=1e-6, k_threshold=1e-10)
+    _validate_physics(λ, layers, Tpp, Tss, Rpp, Rss; Tps=0.0, Tsp=0.0, Rps=0.0, Rsp=0.0, atol=1e-6, k_threshold=1e-10)
 
 Validate physical constraints on R and T values:
 1. Bounds check: 0 ≤ R, T ≤ 1 (always)
 2. Energy conservation: R + T ≈ 1 (for non-absorbing media only)
 
+The energy budget per input polarization includes the cross-polarized channels
+(zero unless the caller passes them): `Rpp + Rps + Tpp + Tps ≈ 1` for p input
+and `Rss + Rsp + Tss + Tsp ≈ 1` for s input.
+
 Issues warnings if constraints are violated.
 
 Internal function called by `transfer` when `validate=true`.
 """
-function _validate_physics(λ, layers, Tpp, Tss, Rpp, Rss; sheets=nothing, atol=1e-6, k_threshold=1e-10)
+function _validate_physics(λ, layers, Tpp, Tss, Rpp, Rss;
+                           sheets=nothing, Tps=0.0, Tsp=0.0, Rps=0.0, Rsp=0.0,
+                           atol=1e-6, k_threshold=1e-10)
     # Check for NaN values (indicates numerical failure)
     if any(isnan, (Tpp, Tss, Rpp, Rss))
         @warn "NaN detected in R/T values" Tpp Tss Rpp Rss
@@ -177,28 +196,27 @@ function _validate_physics(λ, layers, Tpp, Tss, Rpp, Rss; sheets=nothing, atol=
 
     is_lossless = layers_lossless && sheets_lossless
 
+    # Per-input-polarization budgets include the cross-polarized channels
+    # (nonzero only for converting anisotropic/gyrotropic stacks).
+    sum_p = Tpp + Tps + Rpp + Rps
+    sum_s = Tss + Tsp + Rss + Rsp
+
     if is_lossless
         # Energy conservation: R + T = 1 for lossless media
-        sum_p = Tpp + Rpp
-        sum_s = Tss + Rss
-
         if !isapprox(sum_p, 1.0; atol=atol)
-            @warn "Energy conservation violated for p-polarization" Tpp Rpp sum=sum_p expected=1.0 deviation=abs(sum_p - 1.0)
+            @warn "Energy conservation violated for p-polarization" Tpp Tps Rpp Rps sum=sum_p expected=1.0 deviation=abs(sum_p - 1.0)
         end
 
         if !isapprox(sum_s, 1.0; atol=atol)
-            @warn "Energy conservation violated for s-polarization" Tss Rss sum=sum_s expected=1.0 deviation=abs(sum_s - 1.0)
+            @warn "Energy conservation violated for s-polarization" Tss Tsp Rss Rsp sum=sum_s expected=1.0 deviation=abs(sum_s - 1.0)
         end
     else
         # For lossy media: R + T ≤ 1
-        sum_p = Tpp + Rpp
-        sum_s = Tss + Rss
-
         if sum_p > 1.0 + atol
-            @warn "Absorption violation: R + T > 1 for p-polarization" Tpp Rpp sum=sum_p
+            @warn "Absorption violation: R + T > 1 for p-polarization" Tpp Tps Rpp Rps sum=sum_p
         end
         if sum_s > 1.0 + atol
-            @warn "Absorption violation: R + T > 1 for s-polarization" Tss Rss sum=sum_s
+            @warn "Absorption violation: R + T > 1 for s-polarization" Tss Tsp Rss Rsp sum=sum_s
         end
     end
 
diff --git a/test/circular.jl b/test/circular.jl
index 90c4e48..8d7903b 100644
--- a/test/circular.jl
+++ b/test/circular.jl
@@ -49,10 +49,12 @@ end
     # r = (rpp, rps, rss, rsp), t = (tpp, tps, tsp, tss).
     r = SVector{4,ComplexF64}(0.2, 0.0, -0.2, 0.0)
     t = SVector{4,ComplexF64}(0.8, 0.0, 0.0, 0.8)
-    Tpp = 0.96  # Poynting transmittance (N = (0.96+0.96)/(0.64+0.64) = 1.5)
-    Tss = 0.96
+    # Poynting flux factors: unit-amplitude substrate-mode flux 1.5, unit
+    # incident flux 1.0, so each |t|² converts with factor 1.5 (0.64 -> 0.96).
+    f_out1 = f_out2 = 1.5
+    f_in_p = f_in_s = 1.0
 
-    res = TransferMatrix._circular_result(r, t, Tpp, Tss)
+    res = TransferMatrix._circular_result(r, t, f_out1, f_out2, f_in_p, f_in_s)
     @test res isa TransferMatrix.CircularTransferResult{Float64}
 
     # Reflection flips: diagonal ~0, cross = 0.04.
@@ -67,8 +69,9 @@ end
     @test isapprox(res.Trl, 0.0; atol=1e-12)
     @test isapprox(res.Tlr, 0.0; atol=1e-12)
 
-    # N guard: zero transmission amplitude -> all circular T are 0 (not NaN).
-    res0 = TransferMatrix._circular_result(r, SVector{4,ComplexF64}(0, 0, 0, 0), 0.0, 0.0)
+    # Null guard: zero transmission amplitude -> all circular T are 0 (not NaN).
+    res0 = TransferMatrix._circular_result(r, SVector{4,ComplexF64}(0, 0, 0, 0),
+                                           f_out1, f_out2, f_in_p, f_in_s)
     @test res0.Trr == 0.0
     @test res0.Tll == 0.0
     @test !isnan(res0.Trr)
@@ -138,6 +141,38 @@ end
     @test isapprox(res.Rll + res.Rrl + res.Tll + res.Trl, 1.0; atol=1e-6)
 end
 
+@testset "circular budget closes for a polarization-converting stack" begin
+    # Half-wave plate at 45°: strong p<->s conversion in transmission. The old
+    # scalar normalization N = (Tpp+Tss)/(|tpp|²+|tss|²) double-counted the
+    # converted flux (Tpp/Tss were totals) and blew up when the co-polarized
+    # amplitudes vanished at the half-wave point. The per-mode formulation
+    # keeps every per-input circular budget exactly 1 for a lossless stack.
+    no = λ -> 1.658 + 0.0im
+    ne = λ -> 1.486 + 0.0im
+    t_hw = 0.59 / (2 * (1.658 - 1.486))
+    amb = Layer(λ -> 1.0 + 0.0im, 1.0)
+    glass = Layer(λ -> 1.5 + 0.0im, 1.0)
+    wp45 = Layer(no, no, ne, t_hw; euler=(π/4, π/2, 0.0))
+
+    for (sub, θ) in ((amb, 0.0), (amb, 0.3), (glass, 0.3))
+        res = transfer(0.59, [amb, wp45, sub]; θ=θ, basis=:circular)
+        @test isapprox(res.Rrr + res.Rlr + res.Trr + res.Tlr, 1.0; atol=1e-10)
+        @test isapprox(res.Rll + res.Rrl + res.Tll + res.Trl, 1.0; atol=1e-10)
+        for v in (res.Trr, res.Tll, res.Trl, res.Tlr, res.Rrr, res.Rll, res.Rrl, res.Rlr)
+            @test 0.0 <= v <= 1.0
+        end
+    end
+
+    # A half-wave plate flips helicity in transmission: at the half-wave point
+    # the co-handed Trr/Tll are small and the cross-handed Trl/Tlr carry the
+    # transmitted power (linear-basis analogue: Tps ≈ 0.85).
+    res = transfer(0.59, [amb, wp45, amb]; basis=:circular)
+    @test res.Tlr > 0.8
+    @test res.Trl > 0.8
+    @test res.Trr < 0.05
+    @test res.Tll < 0.05
+end
+
 @testset "circular TIR guard" begin
     # ξ = n_inc sinθ > n_substrate -> total internal reflection, no real flux.
     inc = Layer(λ -> 1.8 + 0im, 0.0)
diff --git a/test/exp_propagator.jl b/test/exp_propagator.jl
index 3502374..4558ed6 100644
--- a/test/exp_propagator.jl
+++ b/test/exp_propagator.jl
@@ -104,13 +104,15 @@ using .MuReference
     end
 
     @testset "energy conservation under :exp (lossless)" begin
+        # Both films convert p↔s, so the budgets include the cross-polarized
+        # transmitted channels (per-mode Poynting fluxes).
         amb = Layer(λ -> 1.0, 1.0); sub = Layer(λ -> 1.0, 1.0)
         films = (Layer(λ -> 1.6, λ -> 1.6, λ -> 1.9, 0.4; euler=(π/6, π/4, 0)),  # out-of-plane tilt (#70)
                  Layer(λ -> 1.5, 0.3; mu=gyrotropic_tensor(2.0, 0.6)))           # gyromagnetic μ
         for film in films, θ in (0.0, 0.3, 0.6, 0.9)
             r = transfer(1.0, [amb, film, sub]; θ=θ, method=:exp)
-            @test isapprox(r.Rpp + r.Rps + r.Tpp, 1.0; atol=1e-9)
-            @test isapprox(r.Rss + r.Rsp + r.Tss, 1.0; atol=1e-9)
+            @test isapprox(r.Rpp + r.Rps + r.Tpp + r.Tps, 1.0; atol=1e-9)
+            @test isapprox(r.Rss + r.Rsp + r.Tss + r.Tsp, 1.0; atol=1e-9)
         end
     end
 
diff --git a/test/functions.jl b/test/functions.jl
index ad62cba..7b8ff58 100644
--- a/test/functions.jl
+++ b/test/functions.jl
@@ -464,7 +464,15 @@ end
     sub = Layer(λ->sqrt(2.0), 1.0; mu = 3.0)
     for θ in (0.0, 0.4)
         r = transfer(1.0, [air, film, sub]; θ=θ)
-        @test isapprox(r.Rpp + r.Rps + r.Tpp, 1.0; atol=1e-7)
-        @test isapprox(r.Rss + r.Rsp + r.Tss, 1.0; atol=1e-7)
+        @test isapprox(r.Rpp + r.Rps + r.Tpp + r.Tps, 1.0; atol=1e-7)
+        @test isapprox(r.Rss + r.Rsp + r.Tss + r.Tsp, 1.0; atol=1e-7)
+        # Everything is isotropic (scalar μ), so nothing converts p↔s: the
+        # cross channels must vanish. Guards the degenerate-pair orientation
+        # in calculate_E_modes_tensor — an arbitrarily rotated substrate mode
+        # basis leaks the entire transmitted power into the "cross" label.
+        @test isapprox(r.Tps, 0.0; atol=1e-10)
+        @test isapprox(r.Tsp, 0.0; atol=1e-10)
+        @test isapprox(r.Rps, 0.0; atol=1e-10)
+        @test isapprox(r.Rsp, 0.0; atol=1e-10)
     end
 end
diff --git a/test/integration.jl b/test/integration.jl
index 1c91b0a..27d1241 100644
--- a/test/integration.jl
+++ b/test/integration.jl
@@ -1057,20 +1057,20 @@ end
 
     # Energy conservation for the rotated lossless crystal (issue #70).
     # The (π/6, π/4, 0) rotation tilts the optic axis OUT of the plane of
-    # incidence, coupling p↔s. The correct per-input-polarization budget uses
-    # the reflected cross-pol term and the Poynting transmittance (which
-    # already carries the cross-transmitted power — do NOT also add Tps/Tsp):
-    #   p-input: Rpp + Rps + Tpp = 1     s-input: Rss + Rsp + Tss = 1
+    # incidence, coupling p↔s. The complete per-input-polarization budget
+    # includes BOTH cross-pol channels (each transmittance is the per-mode
+    # Poynting flux of its own output channel):
+    #   p-input: Rpp + Rps + Tpp + Tps = 1   s-input: Rss + Rsp + Tss + Tsp = 1
     # This holds to ~machine precision. The "R+T ≈ 0.998" reported as #70 was
-    # Rpp+Tpp, which simply OMITS the reflected cross-pol Rps (≈0.0022 here) —
-    # a budget-accounting artifact, not a physics error. (Transmission into an
-    # anisotropic *substrate* was a genuine Poynting error, fixed separately in
-    # poynting(); see the "anisotropic substrate energy conservation" testset.)
+    # Rpp+Tpp, which simply OMITS the converted channels — a budget-accounting
+    # artifact, not a physics error. (Transmission into an anisotropic
+    # *substrate* was a genuine Poynting error, fixed separately in poynting();
+    # see the "anisotropic substrate energy conservation" testset.)
     result = transfer(λ, layers)
-    @test isapprox(result.Rpp + result.Rps + result.Tpp, 1.0; atol=1e-9)
-    @test isapprox(result.Rss + result.Rsp + result.Tss, 1.0; atol=1e-9)
-    # The naive co-pol-only budget falls short by exactly the omitted Rps.
-    @test isapprox(result.Rpp + result.Tpp, 1.0 - result.Rps; atol=1e-9)
+    @test isapprox(result.Rpp + result.Rps + result.Tpp + result.Tps, 1.0; atol=1e-9)
+    @test isapprox(result.Rss + result.Rsp + result.Tss + result.Tsp, 1.0; atol=1e-9)
+    # The naive co-pol-only budget falls short by exactly the omitted channels.
+    @test isapprox(result.Rpp + result.Tpp, 1.0 - result.Rps - result.Tps; atol=1e-9)
     @test result.Rps > 1e-4   # cross-pol genuinely nonzero for this out-of-plane tilt
 end
 
@@ -1083,8 +1083,10 @@ end
     # over-counts the flux and breaks energy conservation (≈1.7% at θ=0).
     #
     # Correct lossless budget per input polarization (reflection into the
-    # isotropic ambient is exact |r|²; transmission via Poynting):
-    #   p-input: Rpp + Rps + Tpp = 1     s-input: Rss + Rsp + Tss = 1
+    # isotropic ambient is exact |r|²; each transmittance is its own substrate
+    # eigenmode's Poynting flux, so the two transmitted channels have DIFFERENT
+    # flux-per-|amplitude|² factors here — the case a shared rescale gets wrong):
+    #   p-input: Rpp + Rps + Tpp + Tps = 1   s-input: Rss + Rsp + Tss + Tsp = 1
     λ = 1.0
     air = Layer(λ -> 1.0, 1.0)
     sub = Layer(λ -> 1.5, λ -> 1.7, λ -> 2.0, 1.0; euler=(π/5, π/3, π/7))
@@ -1092,12 +1094,69 @@ end
     for layers in ([air, sub], [air, Layer(λ -> 1.4, 0.4), sub])
         for θ in (0.0, 0.3, 0.6, 0.9, 1.2)
             r = transfer(λ, layers; θ=θ)
-            @test isapprox(r.Rpp + r.Rps + r.Tpp, 1.0; atol=1e-9)
-            @test isapprox(r.Rss + r.Rsp + r.Tss, 1.0; atol=1e-9)
+            @test isapprox(r.Rpp + r.Rps + r.Tpp + r.Tps, 1.0; atol=1e-9)
+            @test isapprox(r.Rss + r.Rsp + r.Tss + r.Tsp, 1.0; atol=1e-9)
         end
     end
 end
 
+@testset "cross-polarized transmittance is a per-mode flux (energy budget)" begin
+    # A half-wave plate with its optic axis in the layer plane at 45° to the
+    # plane of incidence converts nearly all transmitted p light to s. Tps must
+    # be the s-like substrate eigenmode's own Poynting flux ratio — NOT |tps|²,
+    # which is an amplitude ratio. Before the fix, Tpp was the TOTAL transmitted
+    # flux and Tps=|tps|², so the p-input "budget" summed to ~1.85 here.
+    no = λ -> 1.658 + 0.0im
+    ne = λ -> 1.486 + 0.0im
+    λ_hw = 0.59
+    t_hw = λ_hw / (2 * (1.658 - 1.486))     # half-wave thickness
+    amb = Layer(λ -> 1.0 + 0.0im, 1.0)
+    vac_sub = Layer(λ -> 1.0 + 0.0im, 1.0)
+    glass = Layer(λ -> 1.5 + 0.0im, 1.0)
+    wp45 = Layer(no, no, ne, t_hw; euler=(π/4, π/2, 0.0))
+
+    # (1) Energy closure, air/air, normal AND oblique incidence, both methods.
+    # (2) Energy closure with a DIFFERENT exit medium (air -> glass): both
+    #     transmitted channels pick up the substrate flux factor, which a raw
+    #     |t|² misses entirely.
+    for sub in (vac_sub, glass), θ in (0.0, 0.3), method in (:exp, :eig)
+        r = transfer(λ_hw, [amb, wp45, sub]; θ=θ, method=method)
+        @test isapprox(r.Rpp + r.Rps + r.Tpp + r.Tps, 1.0; atol=1e-10)
+        @test isapprox(r.Rss + r.Rsp + r.Tss + r.Tsp, 1.0; atol=1e-10)
+        for v in (r.Rpp, r.Rss, r.Rps, r.Rsp, r.Tpp, r.Tss, r.Tps, r.Tsp)
+            @test 0.0 <= v <= 1.0
+        end
+    end
+
+    # (3) Conversion physics preserved: at the half-wave point nearly all
+    # transmitted light converts; the co-polarized channel is tiny.
+    r = transfer(λ_hw, [amb, wp45, vac_sub])
+    @test isapprox(r.Tps, 0.847; atol=0.01)
+    @test r.Tpp < 1e-3
+    @test isapprox(r.Tsp, r.Tps; atol=1e-10)   # normal incidence symmetry
+
+    # c-cut (optic axis ∥ z): no conversion, cross terms vanish.
+    ccut = Layer(no, no, ne, t_hw; euler=(0.0, 0.0, 0.0))
+    rc = transfer(λ_hw, [amb, ccut, vac_sub]; θ=0.3)
+    @test isapprox(rc.Tps, 0.0; atol=1e-12)
+    @test isapprox(rc.Tsp, 0.0; atol=1e-12)
+    @test isapprox(rc.Rps, 0.0; atol=1e-12)
+    @test isapprox(rc.Rsp, 0.0; atol=1e-12)
+
+    # (4) Lossy converting stack: every channel in [0, 1] and each per-input
+    # budget strictly below 1 (absorption), never above.
+    no_l = λ -> 1.658 + 0.02im
+    ne_l = λ -> 1.486 + 0.02im
+    wp_lossy = Layer(no_l, no_l, ne_l, t_hw; euler=(π/4, π/2, 0.0))
+    rl = transfer(λ_hw, [amb, wp_lossy, glass]; θ=0.3)
+    for v in (rl.Rpp, rl.Rss, rl.Rps, rl.Rsp, rl.Tpp, rl.Tss, rl.Tps, rl.Tsp)
+        @test 0.0 <= v <= 1.0
+    end
+    @test rl.Rpp + rl.Rps + rl.Tpp + rl.Tps < 1.0 - 1e-3
+    @test rl.Rss + rl.Rsp + rl.Tss + rl.Tsp < 1.0 - 1e-3
+    @test rl.Tps > 0.01    # still genuinely converting
+end
+
 @testset "TIR into anisotropic substrate throws (mode-sorting limitation)" begin
     # Pins a KNOWN, intentionally-unsupported edge case (adjacent to documented
     # limitations #71 and #72): total internal reflection from a higher-index
diff --git a/test/mu_reference.jl b/test/mu_reference.jl
index 831b3aa..6a2920f 100644
--- a/test/mu_reference.jl
+++ b/test/mu_reference.jl
@@ -117,10 +117,10 @@ function transfer_exp(stack::Vector{GLayer}, λ; θ=0.0)
     Γs = SMatrix{4,4,ComplexF64}(Γ)
     r, R, t, T = TM.calculate_tr(Γs)
     S = TM.poynting(ξ, q1, qN, γ1, γN, t, r)
-    Tpp, Tss, _, _ = TM.calculate_tr(S)
+    Tpp, Tss, _, _, Tps, Tsp = TM.calculate_tr(S)
     return (; Γ=Γs, r=r, t=t,
             Rpp=R[1], Rss=R[2], Rsp=R[3], Rps=R[4],
-            Tpp=Tpp, Tss=Tss, Tps=T[4], Tsp=T[3])
+            Tpp=Tpp, Tss=Tss, Tps=Tps, Tsp=Tsp)
 end
 
 # Full stack solve using a candidate γ for interior layers (boundaries: package).
diff --git a/test/mu_tensor.jl b/test/mu_tensor.jl
index 7882fb6..1634df8 100644
--- a/test/mu_tensor.jl
+++ b/test/mu_tensor.jl
@@ -122,16 +122,17 @@ const E_TOL  = 1e-7      # energy / duality tolerance
 
     # ---------------------------------------------------------------------
     # (6) Energy conservation for lossless μ-tensor stacks (vacuum/vacuum):
-    # Rpp+Rps+Tpp = 1 and Rss+Rsp+Tss = 1 (Tpp/Tss are TOTAL per-input-pol
-    # Poynting transmittances). Includes gyromagnetic (Hermitian Polder) μ.
+    # Rpp+Rps+Tpp+Tps = 1 and Rss+Rsp+Tss+Tsp = 1 (each transmittance is the
+    # per-output-mode Poynting flux of its own channel). Includes gyromagnetic
+    # (Hermitian Polder) μ, which converts p↔s in transmission.
     # ---------------------------------------------------------------------
     @testset "energy conservation, lossless μ-tensor" begin
         for θ in (0.0, 0.2, 0.5, 0.9),
             μt in (μ_diag(1.5,2.5,3.5), μ_rot(1.5,2.5,3.5,0.4,0.6,0.3),
                    μ_polder(2.0,0.6), μ_polder(2.5,1.2))
             g = transfer_exp([vac(), glay(μ_iso(2.25), μt, 0.4), vac()], 1.0; θ=θ)
-            @test isapprox(g.Rpp + g.Rps + g.Tpp, 1.0; atol=E_TOL)
-            @test isapprox(g.Rss + g.Rsp + g.Tss, 1.0; atol=E_TOL)
+            @test isapprox(g.Rpp + g.Rps + g.Tpp + g.Tps, 1.0; atol=E_TOL)
+            @test isapprox(g.Rss + g.Rsp + g.Tss + g.Tsp, 1.0; atol=E_TOL)
         end
     end
 
@@ -192,3 +193,18 @@ end
     H = hfield(1.0, layers; θ=0.3, dz=0.01)
     @test all(isfinite, H.p) && all(isfinite, H.s)
 end
+
+@testset "_orient_degenerate_pair: span without ŷ passes through" begin
+    # A degenerate 2-D eigenmode space with no Ey content has no s-like
+    # direction to orient to, so the input basis must come back unchanged —
+    # up to the p-mode sign convention (Ex real-negative for the backward
+    # pair, matching E_modes[3,1] = -1 in calculate_E_modes).
+    x̂ = SVector{3,ComplexF64}(1, 0, 0)
+    ẑ = SVector{3,ComplexF64}(0, 0, 1)
+    w_p, w_s = TransferMatrix._orient_degenerate_pair(x̂, ẑ, 1)
+    @test w_p == x̂
+    @test w_s == ẑ
+    w_p, w_s = TransferMatrix._orient_degenerate_pair(x̂, ẑ, 3)
+    @test w_p == -x̂
+    @test w_s == ẑ
+end