Internals

_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.

_propagate_core_exp(λ, layers; θ=0.0, μ=1.0, sheets=nothing)

Matrix-exponential propagation core. Interior layers propagate via layer_transfer_exp (no eigenmode sorting); the semi-infinite ambient and substrate keep the eigenmode treatment in layer_matrices, which is needed for the r/t coefficients and the Poynting transmittance. Conductive sheets are injected at their respective interfaces, mirroring _propagate_core. Returns (M_sys, S) like _propagate_core.

_propagate_full(λ, layers; θ=0.0, μ=1.0, sheets=nothing)

Internal full transfer-matrix pass. Returns (M_sys, S, Ds, Ps, E_modes_per_layer, qs) — like propagate but also returns the per-layer eigenvalue vectors qs, needed for magnetic-field reconstruction. Supports conductive sheets (Task: sheets).

propagate(λ, layers; θ=0.0, μ=1.0, sheets=nothing)

Calculate the transfer matrix and Poynting vector for the structure, plus the per-layer D, P, and E_modes matrices used for field reconstruction. Returns the 5-tuple (M_sys, S, Ds, Ps, E_modes_per_layer). See transfer for the public R/T API.

Ratio of the absolution squares of two components used to evaluate if a material is birefringent.

poynting(k_par, q_in, q_out, E_modes_in, E_modes_out, t_coefs, r_coefs[, μ_in, μ_out])

Calculate the Poynting vector from wavevectors $q$, components of the electric field E_modes, and transmission and reflection coefficients.

Transmitted Poynting vectors use substrate wavevectors (q_out), while reflected Poynting vectors use incident-medium wavevectors (k_in[3,:], k_in[4,:]), since reflected waves propagate in the incident medium.

The optional μ_in and μ_out arguments are the 3×3 permeability tensors of the ambient (incident) and substrate media respectively (default: identity, i.e. non-magnetic). They are used to form $H = μ^{-1}(k̄ \times E)$ when computing the Poynting flux, so that energy conservation holds for a magnetic substrate. A magnetic ambient is only correctly handled at normal incidence; at oblique incidence the conserved in-plane wavevector k_par is still computed from the ambient permittivity alone (see issue #71), so results may be inaccurate.

poynting(Ψ, a)

Calculates the Poynting vector for the structure from Ψ and matrix $a$.

From Berreman, 1972, Ψ is the column matrix:

\[\Psi = \begin{pmatrix} Ex \\Hy \\Ey \\-Hx \end{pmatrix}\]

for a right-handed Cartesian coordinate system with the z-axis along the normal to the multilayer structure.

Berreman, 1972, https://doi.org/10.1364/JOSA.62.000502

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.

Base.broadcastable(layer::Layer)

Treat a single Layer as a scalar under broadcasting so helpers can be broadcast over one layer across many wavelengths without an explicit Ref, e.g. get_refractive_indices.(layer, λs). Without this, broadcasting a bare Layer hits Julia's collect-based fallback and throws a confusing MethodError: no method matching length(::Layer).

A Vector{Layer} stack still broadcasts element-wise, which is the desired behavior for per-layer operations. For threaded wavelength/angle/thickness sweeps, prefer the dedicated sweep_angle / sweep_thickness API.

_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.

calculate_E_modes(k_par, q, ε, μ)

The 4 x 3 matrix E_modes contains vector components that belong to the electric field calculated such that singularities are identified and removed.

q[1] and q[2] are forward-traveling modes and q[3] and q[4] are backward-traveling modes.

Each row v of E_modes is normalized to unit length using the Hermitian norm $‖v‖ = √(Σ|vᵢ|²)$ as prescribed by the 2019 erratum (Passler & Paarmann, JOSAB 36, 3246). The Hermitian norm is the correct choice because each row represents an electric field polarization direction: its physical "length" is the field amplitude, which is $√(E⋅E*) = √(Σ|Eᵢ|²)$. The bilinear form $√(Σ vᵢ²)$ has no physical meaning for complex fields and can vanish for nonzero vectors (e.g. circular polarization), making it unsuitable as a normalization factor.

This is based on the work in Xu, et al., 2000, https://doi.org/10.1103/PhysRevB.61.1740

calculate_E_modes_tensor(k_par, q, ε, μ)

E-field eigenvectors for a layer with a full 3×3 permeability tensor μ, as the 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. Degenerate pairs are oriented to the fixed p/s polarization references of the scalar path (see _orient_degenerate_pair) 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].

calculate_q(Δ, a)

The four eigenvalues of $q$ may be obtained from the 4x4 matrix Δ and then eigenvectors may be found for each eigenvalue. Here the eigenvalues must be sorted appropriately to avoid potentially discontinuous solutions. This extends from the work in

Li et al, 1988, https://doi.org/10.1364/AO.27.001334

construct_a(k_par, C)

Construct the elements of the intermediate 6×6 matrix $a$ in terms of the elements of matrix $C$ (the 6×6 constitutive matrix holding the material dielectric and permeability tensors) and propagation vector k_par. This is implemented as described in

Berreman, 1972, https://doi.org/10.1364/JOSA.62.000502

construct_constitutive(ε, μ, ρ1, ρ2)

Construct the 6×6 constitutive matrix C = [ε ρ1; ρ2 μ] from the dielectric and permeability tensors.

construct_Δ(k_par, C, a)

Construct the reordered matrix Δ in terms of the elements of the two matrices, C and a, and the in-plane reduced wavevector kpar = ``kx / k_0``. The matrix Δ is involved in the relation

\[ \frac{\delta}{\delta z}\Psi = \frac{i \omega}{c}\Delta \Psi\]

and Δ is the reordered S matrix in Berreman's formulation.

Berreman, 1972, https://doi.org/10.1364/JOSA.62.000502

dynamical_matrix(k_par, q, E_modes, μ)

The dynamical matrix relating two layers at the interface where matrix $A_i$ for layer $i$ relates the field $E_i$ to the field in the previous layer $i - 1$ via

\[ A_{i-1}E_{i-1} = A_{i}E_{i}\]

Xu, et al., 2000, https://doi.org/10.1103/PhysRevB.61.1740

dynamical_matrix(k_par, q, E_modes, μ::AbstractMatrix)

Dynamical matrix for a tensor μ. Rows are (Eₓ, Eᵧ, Hᵧ, −Hₓ) per mode with H = μ⁻¹ (k̄ × E). Reduces to the scalar-μ method when μ = μ·I.

layer_matrices(ω, k_par, layer, μ)

Calculate all parameters for a single layer, particularly the propagation matrix and dynamical matrix so that the overall transfer matrix can be calculated.

Supports both isotropic layers (single refractive index) and anisotropic layers (different refractive indices along principal axes), with optional crystal rotation via Euler angles.

layer_transfer_exp(layer, λ, k_par, ω, μ_i)

Interior-layer 4×4 transfer matrix in the dynamical-matrix field basis (Eₓ, Eᵧ, Hᵧ, −Hₓ), computed as the matrix exponential of the Berreman Δ matrix rather than by eigenmode decomposition:

\[T = swap23\, \exp\!\left(-i\frac{ω}{c}\,Δ\,d\right) swap23.\]

Δ = construct_Δ(k_par, construct_constitutive(ε, μ), a) is built from the full 3×3 ε (with any Euler rotation) and μ (the scalar fallback μ_i·I, or the layer's tensor mu). Because Δ's eigenvalues are the mode wavevectors q, exp(-i(ω/c)Δd) equals the eigenmode propagator D·P(d)·D⁻¹ without diagonalizing, so this path needs no eigenvalue sorting and is degeneracy-immune. swap23 reorders between Berreman's field vector (Eₓ, Hᵧ, Eᵧ, −Hₓ) and the dynamical-matrix basis; the exp(-iωt) sign matches propagation_matrix.

The matrix exponential is evaluated by scaling-and-squaring with degree-13 Padé approximants (StaticArrays' exp).

Berreman, 1972, https://doi.org/10.1364/JOSA.62.000502 Mackay & Lakhtakia, 2020, https://doi.org/10.1007/978-3-031-02022-3 Higham, 2005, https://doi.org/10.1137/04061101X

permeability_tensor(μ1, μ2, μ3)

This produces the diagonal permeability tensor, which is identical to the way we build the dielectric_tensor, and we include this function simply for completeness.

propagation_matrix(ω, q)

Returns a callable object that propagates the electromagnetic field a distance z through a material for a frequency ω and wavevector q.

The propagation matrix is diagonal with elements:

P(z) = diag(exp(-iωq₁z/c), exp(-iωq₂z/c), exp(-iωq₃z/c), exp(-iωq₄z/c))

This uses the exp(-iωt) time convention, consistent with Berreman (1972) and Passler & Paarmann (2017). Note that Yeh uses exp(+iωt), which would flip the sign in the exponent.

sheet_matrix(sheet::Sheet, λ)

Return the 4×4 interface matrix G for a conductive sheet at wavelength λ, in the dynamical-matrix field basis (Eₓ, Eᵧ, H_y, −Hₓ). This is the inverse of the physical field-jump matrix, G = G_phys⁻¹, so its off-diagonal (H-row) entries are +σ̃ with the dimensionless conductivity σ̃ = Z₀ σ (Z₀ = √(μ₀/ε₀)). It is injected at an interface as Dᵢ⁻¹ G Dᵢ₊₁. Tangential E (rows 1,2) is continuous; tangential H (rows 3,4) jumps.

Z₀ is applied here and nowhere else (see Sheet).

_to_eV(x)

Convert a dispersion-model frequency/energy parameter (ω_p, ω_0, γ) to electron-volts (the internal unit for the built-in dispersion models).

No-op for plain Real inputs (already assumed to be eV). The UnitfulExt extension adds methods for unit-bearing energy (eV, J), spectroscopic wavenumber (cm^-1), and ordinary frequency (Hz), each mapped linearly to energy so resonance positions and damping widths convert identically.

_to_radians(x)

Convert an angle to radians (the package's internal angle unit).

No-op for plain Real inputs (already assumed to be radians). The UnitfulExt extension adds a method for unit-bearing angles (e.g. 45u"°") that strips to radians. Used to normalize the incidence angle θ at every public entry point.

_to_um(x)

Convert a length quantity to micrometers (the package's internal length unit).

No-op for plain Real inputs (already assumed to be μm). The UnitfulExt extension adds a method for Unitful.Length that strips to μm. Used to normalize layer thickness and the field-sampling step dz.

_to_wavelength_um(x)

Convert a spectral input to wavelength in micrometers.

No-op for plain Real inputs (already assumed to be μm). The UnitfulExt extension adds a method that accepts a Unitful length, spectroscopic wavenumber (ν̃ = 1/λ), ordinary frequency (λ = c/f), or photon energy (λ = hc/E). Used to normalize the wavelength λ at every public entry point.

HC_EV_UM

Product hc/e in eV·μm (CODATA 2018). Converts a vacuum wavelength in μm to photon energy in eV via E = HC_EV_UM / λ.