Public Documentation

amplitude_coefficients(λ, layers; θ=0.0, μ=1.0, sheets=nothing, method=:exp)

Complex amplitude reflection and transmission coefficients of the stack, returned as a named tuple (; rpp, rps, rsp, rss, tpp, tps, tsp, tss). Arguments match transfer, which reports the corresponding intensities.

Each coefficient is xᵢₒ = amplitude of the outgoing mode o per unit amplitude of the incident mode i, in the package's polarization basis: the first index is the input polarization and the second the output (rps is p in → s out). The amplitudes multiply the unit-norm mode vectors of ModeAmplitudes: in the incident medium the amplitude 4-vector at z = 0 is (1, 0, rpp, rps) for p incidence and (0, 1, rsp, rss) for s, and in the exit medium (tpp, tps, 0, 0) / (tsp, tss, 0, 0) at the last interface.

Sign convention

The backward p mode vector is −x̂ at normal incidence (E_modes[3,:] = (−1,0,0)), while both s mode vectors are +ŷ. At normal incidence on an isotropic stack the two polarizations are physically identical, so rpp = −rss and tpp = tss: rss is the reflection coefficient of the tangential electric field itself (it tends to −1 for a perfect conductor), and rpp carries the extra sign of the p basis vector.

Relation to intensities

Rpp = |rpp|² etc. (exact for a transparent incident medium). Transmittances are Poynting-flux ratios, T = |t|² Re(q_out)/q_in for isotropic lossless media (q = n cos θ), so Tss = |tss|² only when the exit medium equals the incident medium; see transfer.

sweep_angle(λs, θs, layers; threads=true, verbose=false, basis=:linear)

Calculate transmittance/reflectance spectra over wavelength and angle of incidence.

Returns a TransferResult with fields Tpp, Tss, Rpp, Rss, each a matrix of size (length(θs), length(λs)).

Arguments

• λs: Vector of wavelengths in μm
• θs: Vector of angles of incidence in radians
• layers: AbstractVector{<:Layer} representing the stack
• threads: Enable multithreading (default: true)
• verbose: Print thread count info (default: false)
• basis: :linear (default) → TransferResult; :circular → CircularTransferResult (see transfer)
• method: propagation backend, :exp (default) or :eig (see transfer)

Units

• Wavelengths: μm (micrometers) recommended
• Angles: radians

sweep_thickness(λs, ts, layers, t_index; θ=0.0, threads=true, verbose=false, basis=:linear)

Sweep the thickness of a specific layer and calculate transmittance/reflectance spectra.

Returns a TransferResult with fields Tpp, Tss, Rpp, Rss, each a matrix of size (length(ts), length(λs)).

Arguments

• λs: Vector of wavelengths in μm
• ts: Vector of thicknesses in μm to sweep
• layers: AbstractVector{<:Layer} representing the stack
• t_index: Index of the layer (1-based) whose thickness to vary
• θ: Angle of incidence in radians (default: 0.0, normal incidence)
• threads: Enable multithreading (default: true)
• verbose: Print thread count info (default: false)
• basis: :linear (default) → TransferResult; :circular → CircularTransferResult (see transfer)
• method: propagation backend, :exp (default) or :eig (see transfer)

Units

• Wavelengths and thicknesses: μm (micrometers) recommended
• Angle: radians

transfer(λ, layers; θ=0.0, μ=1.0, validate=false, basis=:linear, method=:exp)

Calculate the transmittance and reflectance of a layered structure.

With the default basis=:linear, returns a TransferResult with fields Tpp, Tss, Rpp, Rss (and cross-polarization terms Tps, Tsp, Rps, Rsp). With basis=:circular, returns a CircularTransferResult.

Reflectance and transmittance calculation

Reflectance uses $R = |r|^2$ from the transfer matrix coefficients (Passler & Paarmann 2017, Eq. 17). This is exact for transparent incident media. For absorbing incident media, $|r|^2$ is not a true energy ratio — Poynting vectors become non-additive due to interference cross-terms between incident and reflected waves (Ortiz & Mochan 2005, JOSA A 22, 2827). Proper treatment of that case requires the power_entering formalism (Byrnes 2016, arXiv:1603.02720), which is not yet implemented here.

Transmittance uses Poynting vectors (energy flux ratio $S_out / S_in$) 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
• θ: Angle of incidence in radians (default: 0.0, normal incidence)
• μ: Relative magnetic permeability (default: 1.0, non-magnetic)
• validate: Check energy conservation R + T ≈ 1 for non-absorbing media (default: false)
• basis: Output polarization basis — :linear (default) or :circular
• method: propagation backend — :eig (eigenmode) or :exp (matrix exponential, see layer_transfer_exp)

Numerical backend

method selects how interior layers are propagated:

• :exp (default) computes each interior layer's transfer matrix as the matrix exponential of the Berreman Δ matrix (see layer_transfer_exp). It needs no eigenmode sorting and is degeneracy-immune, so it handles near-degenerate and mixed propagating/evanescent interior layers that :eig cannot.
• :eig is the eigenmode/dynamical-matrix path, retained as a cross-check.

Both agree to ~1e-12 on all supported cases. The semi-infinite ambient and substrate use the eigenmode path in either mode, so the anisotropic-ambient (#71) and anisotropic-substrate (#107) boundary limitations are unaffected by method.

Mackay & Lakhtakia, 2020, https://doi.org/10.1007/978-3-031-02022-3 ; Higham, 2005, https://doi.org/10.1137/04061101X

Polarization basis

• basis=:linear (default) returns a TransferResult in the linear p/s basis.

• basis=:circular returns a CircularTransferResult in the right/left circular basis. R/L are fixed-lab-frame helicities under this package's exp(-iωt) convention. The Jones matrices transform as r_circ = C⁻¹ r_lin C with C = (1/√2)[1 1; i -i] (columns ordered L, R); the helicity flip on reflection is encoded automatically in the opposite signs of rpp/rss.

  Circular reflectance is |r_circ|² (a true energy ratio under the same condition as linear R: transparent incident medium, cf. issue #72). Circular 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)
• The first and last layers are treated as semi-infinite media
• θ is measured from the surface normal (z-axis)

Units

• Wavelength and thicknesses: μm by default. With using Unitful, λ may carry units — a length, or a wavenumber/frequency/photon-energy that maps to wavelength (e.g. 1550u"nm", 193u"THz", 0.8u"eV").
• Angle: radians
• Transmittance/Reflectance: dimensionless (0 to 1)

Physics Validation

When validate=true, the function checks:

1. Bounds: 0 ≤ R, T ≤ 1 (catches NaN, negative values, numerical instability)
2. Energy conservation: R + T ≈ 1 for non-absorbing media (imag(n) < 1e-10)
3. Absorption bound: R + T ≤ 1 for absorbing media

Warnings are issued for any violations.

calculate_tr(M_sys)

Read the reflection and transmission coefficients off the full-stack transfer matrix M_sys.

M_sys is the $4\times4$ transfer matrix for the whole stack, reordered into the p/s block convention (the two p modes first, then the two s modes) so the 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$ 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 generalized $4\times4$ coefficient formulas follow Passler & Paarmann (2017).

References:

• Yeh, 1979, https://doi.org/10.1364/JOSA.69.000742
• Passler & Paarmann, 2017, https://doi.org/10.1364/JOSAB.34.002128

efield(λ, layers; θ=0.0, μ=1.0, dz=0.001, sheets=nothing)

Calculate the electric field profile throughout the layered structure.

Returns an ElectricField struct containing:

• z: Position coordinates along the structure
• p: Electric field components (Ex, Ey, Ez) for p-polarized incidence
• s: Electric field components (Ex, Ey, Ez) for s-polarized incidence
• boundaries: z-positions of layer interfaces

Arguments

• λ: Wavelength in μm (must match units used for layer thicknesses)
• layers: Vector of Layer objects representing the stack
• θ: Angle of incidence in radians (default: 0.0, normal incidence)
• μ: Relative magnetic permeability (default: 1.0, non-magnetic)
• dz: Spatial step size in μm for field sampling (default: 0.001)
• sheets: Optional conductive sheets at interfaces (see Sheet and transfer); keyed by the index of the layer above each interface.

Wave Propagation Convention

• Light propagates in the +z direction (from first layer toward last layer)
• z = 0 is at the first interface (between layer 1 and layer 2)
• Negative z values are inside the first (incident) layer
• θ is measured from the surface normal (z-axis)

Units

• All lengths (λ, thickness, dz, z): μm by default; with using Unitful, λ and dz may carry units (e.g. efield(1.55u"μm", layers; dz=1u"nm")).
• Angle: radians
• Electric field: arbitrary units (normalized to incident field)

hfield(λ, layers; θ=0.0, μ=1.0, dz=0.001, sheets=nothing)

Calculate the magnetic-field profile through the structure, returning a MagneticField. Shares the sampling grid with efield (same dz and arguments), so E and H can be overlaid; the in-plane H discontinuity at a conductive sheet equals the surface current ẑ × (σ_s E∥).

Units / normalization

H is returned in impedance-normalized units H̃ = Z₀ H_SI (Z₀ = √(μ₀/ε₀)), so |E| ~ |H̃| for a plane wave. Arguments and conventions match efield.

LayerEnergies

Closed-form electric-energy integrals of the finite layers of a stack, as returned by layer_energies. Element j refers to layers[j+1] (the two semi-infinite media are excluded).

Fields

• p::Vector{Float64}: ∫ E† ε_H E dz over each finite layer for p incidence
• s::Vector{Float64}: the same for s incidence
• boundaries::Vector{Float64}: z-positions of the interfaces (as in efield)

ModeAmplitudes

Per-layer Berreman mode amplitudes of a stack for unit-amplitude p- and s-polarized incidence, as returned by mode_amplitudes.

Fields

• p::Vector{SVector{4,ComplexF64}}: amplitude 4-vector of each layer for p incidence
• s::Vector{SVector{4,ComplexF64}}: amplitude 4-vector of each layer for s incidence
• qs::Vector{SVector{4,ComplexF64}}: reduced normal wavevector q = k_z/k₀ of the four modes of each layer (k₀ = 2π/λ)
• E_modes_per_layer::Vector{SMatrix{4,3,ComplexF64}}: unit-norm electric-field polarization vector of each mode (rows) in each layer
• boundaries::Vector{Float64}: z-positions of the N-1 interfaces (boundaries[1] = 0)
• k_par::ComplexF64: reduced in-plane wavevector n₁ sin θ

Slot convention

Slots 1 and 2 are the forward (+z) modes — p-like then s-like — and slots 3 and 4 the backward (−z) modes, again p-like then s-like. For an isotropic layer at normal incidence only slots 1 and 3 are populated for p incidence (E_modes[1,:] = (1,0,0), E_modes[3,:] = (−1,0,0)) and only slots 2 and 4 for s incidence (E_modes[2,:] = E_modes[4,:] = (0,1,0)); the backward p basis vector points along −x, so a p-incidence amplitude a₃ contributes −a₃ to Ex.

Reference plane and propagation sign

The field inside layer i is

\[\mathbf{E}_i(z) = \sum_{m=1}^{4} a_m \, \mathbf{E}_m \, e^{\,i k_0 q_m (z - z_i)},\]

with the exp(−iωt) time convention (forward modes have Re q > 0, or Im q > 0 when evanescent). The reference plane z_i is the layer's entry (left) face, z_i = boundaries[i-1], for every layer i ≥ 2 — including the semi-infinite exit medium, whose amplitudes are (tpp, tps, 0, 0) for p and (tsp, tss, 0, 0) for s incidence. The semi-infinite incident medium (i = 1) has no entry face and is referenced to its exit face z₁ = boundaries[1] = 0, where its amplitudes are (1, 0, rpp, rps) for p and (0, 1, rsp, rss) for s incidence (see amplitude_coefficients).

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

Closed-form integral of the electric energy density over every finite layer of the stack, for unit-amplitude p and s incidence, returned as a LayerEnergies (vectors of length length(layers) - 2; element j is layers[j+1]).

For each finite layer with entry-face amplitudes a from mode_amplitudes the field is E(z) = Σₘ aₘ Eₘ exp(i k₀ qₘ z) on 0 ≤ z ≤ d, and

\[U = \int_0^d \mathbf{E}^\dagger\, ε_H\, \mathbf{E}\, dz = \sum_{m,m'} \bar a_m a_{m'} \left(\mathbf{E}_m^\dagger ε_H \mathbf{E}_{m'}\right) \int_0^d e^{\,i k_0 (q_{m'} - \bar q_m) z}\, dz ,\]

where the z-integral is evaluated analytically (it reduces to d when q_{m'} = conj(q_m), e.g. the self-terms of a lossless mode or the cross terms of an evanescent pair) and ε_H = (ε + ε†)/2 is the Hermitian part of the layer's dielectric tensor in the lab frame — the same tensor the transfer matrix uses (dielectric_constant(n) on each principal axis, rotated by the layer's Euler angles). For a symmetric ε this is Re ε, so for an isotropic lossless layer U = n² ∫|E|² dz and for an absorbing one U = Re(n²) ∫|E|² dz.

Units

E is normalized to unit incident amplitude (the incident mode vector has unit norm) and z is in μm, so U is in μm times the (dimensionless) relative permittivity: a vacuum layer returns ∫|E|² dz in μm. Ratios such as the electric-energy fraction of a cavity mode, U_gap / Σ U, are dimensionless.

Example

U = layer_energies(λ, layers)
f_E = U.p[gap - 1] / sum(U.p)     # gap = index of the cavity layer in `layers`

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

Per-layer Berreman mode amplitudes for unit-amplitude p and s incidence, returned as a ModeAmplitudes. These are the coefficients that efield and hfield expand into field profiles; the struct docstring states the slot ordering (forward p, forward s, backward p, backward s), the reference plane of each layer (entry face, except the incident medium which is referenced to z = 0) and the propagation sign convention (exp(+i k₀ q (z − z_ref)), exp(−iωt) time dependence).

Arguments and units match efield: λ in μm (or a Unitful quantity), θ in radians, sheets as in transfer. The eigenmode backend is used, as for efield.

Example

A = mode_amplitudes(1.0, layers)
A.p[1]           # (1, 0, rpp, rps) — incident medium at z = 0
A.p[end]         # (tpp, tps, 0, 0) — exit medium at the last interface
a = A.s[2]       # s incidence, layer 2, referenced to its entry face
# E_y(z) in layer 2:  a[2] exp(i k₀ q₂ (z - z₁)) + a[4] exp(i k₀ q₄ (z - z₁))

calculate_tr(S::Poynting)

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, 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 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), so the negative sign in Rpp = -S.refl_p[3] / S.in_p[3] yields positive reflectance.

CircularTransferResult{T}

Circular-polarization (RCP/LCP) basis transmittance/reflectance — the analogue of TransferResult. Fields are right/left co- and cross-handedness terms: Trr is T_{R←R}, Trl is T_{R←L}, etc.

As with TransferResult, T is Float64 for a single transfer(λ, layers; basis=:circular) call and Matrix{Float64} for the sweep functions.

Convention

R/L are viewer-facing helicities defined in the fixed lab frame relative to +z, under this package's exp(-iωt) time convention (basis states e_R = (x̂ - iŷ)/√2, e_L = (x̂ + iŷ)/√2). Reflection flips the helicity label, so at normal incidence on an isotropic interface the diagonal Rrr/Rll vanish and reflection appears entirely in the off-diagonal Rrl/Rlr; at oblique incidence the diagonal terms are small but nonzero. Unequal Rrl ≠ Rlr signals genuine optical activity / magneto-optic coupling, not a bug.

See transfer for energy-ratio caveats.

ElectricField

Spatial electric-field profile through a layered structure, mirrored by MagneticField.

• z: position coordinates
• p: (Ex, Ey, Ez) for p-polarized incidence
• s: (Ex, Ey, Ez) for s-polarized incidence
• boundaries: z-positions of interfaces

MagneticField

Spatial magnetic-field profile through a layered structure, mirroring ElectricField. Fields are in impedance-normalized units H̃ = Z₀ H_SI so |E| ~ |H̃| for a plane wave and E/H can be overlaid directly.

• z: position coordinates (same grid as efield)
• p: (Hx, Hy, Hz) for p-polarized incidence
• s: (Hx, Hy, Hz) for s-polarized incidence
• boundaries: z-positions of interfaces

TransferResult{T}

Container for reflectance and transmittance results from transfer matrix calculations.

Fields

• Tpp::T: p-polarized transmittance (p-in → p-out)
• Tss::T: s-polarized transmittance (s-in → s-out)
• Tps::T: cross-polarized transmittance (p-in → s-out)
• Tsp::T: cross-polarized transmittance (s-in → p-out)
• Rpp::T: p-polarized reflectance (p-in → p-out)
• Rss::T: s-polarized reflectance (s-in → s-out)
• Rps::T: cross-polarized reflectance (p-in → s-out)
• Rsp::T: cross-polarized reflectance (s-in → p-out)

For single-wavelength calculations via transfer(), T is Float64. For sweep calculations via sweep_angle() or sweep_thickness(), T is Matrix{Float64}.

Examples

# Single wavelength - returns TransferResult{Float64}
result = transfer(1.0, layers)
result.Tpp  # Float64

# Sweep - returns TransferResult{Matrix{Float64}}
result = sweep_angle(λs, θs, layers)
result.Tpp  # Matrix{Float64}

# Destructuring works
(; Tpp, Rpp) = transfer(1.0, layers)

Layer(material, thickness)
Layer(nx, ny, nz, thickness)

Construct a single layer for transfer matrix calculations.

Isotropic Layer (single refractive index)

Layer(material, thickness)

• material: a dispersion function λ -> n(λ), or a RefractiveMaterial from RefractiveIndex.jl (the material form requires using RefractiveIndex, which activates the RefractiveIndexExt extension)
• thickness: Layer thickness in the same units as wavelength (typically μm)

Anisotropic Layer (biaxial: three refractive indices)

Layer(nx, ny, nz, thickness)

• nx, ny, nz: Dispersion functions λ -> n(λ) for each principal axis
• thickness: Layer thickness in μm

For uniaxial materials, set two axes equal (e.g., nx = ny for optic axis along z).

Units Convention

Lengths default to micrometers (μm) as bare numbers (Layer(n, 0.1) = 0.1 μm). With using Unitful, thickness may carry units and is normalized to μm: Layer(n, 100u"nm"). RefractiveIndex.jl dispersion functions are parameterized in μm, so μm remains the internal unit.

Layer is parametric as Layer{F,Mf,T} where F is the dispersion function type, Mf is the permeability type (Nothing or a callable λ -> SMatrix{3,3,ComplexF64}), and T is the thickness type. For anisotropic layers, F is a Tuple of three dispersion functions.

Examples

# Isotropic: Using a RefractiveIndex.jl material (requires `using RefractiveIndex`)
n_sio2 = RefractiveMaterial("main", "SiO2", "Malitson")
layer = Layer(n_sio2, 0.1)  # 100 nm = 0.1 μm

# Isotropic: Using custom dispersion function
layer = Layer(λ -> 1.5 + 0.01im, 0.25)  # constant n = 1.5 + 0.01i

# Anisotropic: Uniaxial crystal (calcite-like, optic axis along z)
no = λ -> 1.658  # ordinary index
ne = λ -> 1.486  # extraordinary index
layer = Layer(no, no, ne, 0.5)

# Anisotropic: Biaxial crystal
layer = Layer(λ -> 1.5, λ -> 1.6, λ -> 1.7, 0.3)

Magnetic / Permeability (mu=)

All constructors accept an optional mu= keyword that sets the layer's magnetic permeability. Accepted forms: nothing (default, use the global μ= fallback), a scalar Number (isotropic μ), a constant 3×3 AbstractMatrix, or a callable λ -> 3×3 matrix for wavelength-dependent tensors.

# Gyromagnetic (Polder) ferrite film, bias ∥ z
film = Layer(λ -> 1.0, 0.4; mu = gyrotropic_tensor(2.0, 0.6))

# Uniaxial-μ slab
slab = Layer(λ -> 1.5, 0.3; mu = [2 0 0; 0 2 0; 0 0 3])

Layer(nx, ny, nz, thickness; euler=(0,0,0))

Construct an anisotropic layer with different refractive indices along each principal axis.

Arguments

• nx, ny, nz: Dispersion functions λ -> n(λ) for each principal axis
• thickness: Layer thickness in μm
• euler: Optional tuple (φ, θ, ψ) of ZYZ Euler angles in radians (default: no rotation)

The dielectric tensor is constructed as a diagonal matrix with εi = ni² in the crystal frame, then rotated to the lab frame using the Euler angles.

Euler Angle Convention (ZYZ)

• φ (phi): First rotation about z-axis
• θ (theta): Rotation about new y-axis (tilt angle from z)
• ψ (psi): Second rotation about new z-axis

Examples

# Uniaxial crystal with optic axis along z (no rotation needed)
layer = Layer(no, no, ne, 0.5)

# Same crystal with optic axis tilted 45° from z in the xz-plane
layer = Layer(no, no, ne, 0.5; euler=(0, π/4, 0))

# Optic axis in the xy-plane at 30° from x
layer = Layer(no, no, ne, 0.5; euler=(π/6, π/2, 0))

dielectric_constant(n_re::Real, n_im::Real)

Return the complex dielectric function from the real and imaginary parts of the index of refraction.

The complex index of refraction, given by

    n' = n_re + i * n_im

(in terms of nre and nim), can be used to obtain the frequency-dependent complex dielectric function

    ε_r(ω) = ε' + iε''

via the relation

    (n_re + i * n_im)^2 = ε' + iε''.

find_bounds(layers)

Return cumulative thickness positions measured from the start of layer 1, along with the total thickness.

The returned vector has one entry per layer, where entry i is the sum of thicknesses from layer 1 through layer i. These are not interface positions relative to z = 0; callers that need z-coordinates (e.g., efield) must subtract the first layer thickness.

get_euler_angles(layer::Layer)

Return the Euler angles (φ, θ, ψ) for crystal rotation.

For isotropic layers or unrotated anisotropic layers, returns (0.0, 0.0, 0.0). For rotated anisotropic layers, returns the stored angles.

get_permeability(layer::Layer, λ)

Return the 3×3 permeability tensor SMatrix{3,3,ComplexF64} for the layer at wavelength λ, or nothing if the layer is non-magnetic.

get_refractive_indices(layer::Layer, λ)

Return the refractive indices for a layer at wavelength λ.

For isotropic layers, returns (n, n, n) where n is the scalar refractive index. For anisotropic layers, returns (nx, ny, nz) for each principal axis.

isanisotropic(layer::Layer)

Return true if the layer has anisotropic optical properties (different refractive indices along principal axes), false for isotropic layers.

ismagnetic(layer::Layer)

Return true if the layer has a non-trivial magnetic permeability (mu ≠ nothing).

isrotated(layer::Layer)

Return true if the layer has non-zero Euler angles for crystal rotation.

refractive_index(λs, ns, ks)

Build a complex dispersion function λ -> n(λ) + i·k(λ) by linearly interpolating tabulated real (ns) and imaginary (ks) refractive-index data over wavelengths λs (sorted, increasing). Evaluation outside the tabulated range throws a DomainError.

A refractive_index(material::RefractiveMaterial) method that derives the dispersion function from a RefractiveIndex.jl material is provided by the RefractiveIndexExt package extension — load RefractiveIndex to enable it (along with the Layer(::RefractiveMaterial, d) and Sheet(::RefractiveMaterial, d) constructors).

dielectric_tensor(ε1, ε2, ε3)

Return the diagonal complex dielectric tensor

\[\varepsilon = \begin{pmatrix} \varepsilon_1 & 0 & 0 \\0 & \varepsilon_2 & 0 \\0 & 0 & \varepsilon_3 \end{pmatrix}\]

euler_rotation_matrix(φ, θ, ψ)

Return the 3×3 rotation matrix for ZYZ Euler angles (in radians).

This transforms vectors from the crystal frame to the lab frame: v_lab = R * v_crystal

The rotation is performed as: R = Rz(φ) * Ry(θ) * Rz(ψ)

Convention

• φ (phi): First rotation about z-axis (0 to 2π)
• θ (theta): Rotation about new y-axis (0 to π) - the tilt angle
• ψ (psi): Second rotation about new z-axis (0 to 2π)

Common cases

• Optic axis along z: (0, 0, 0) - no rotation needed
• Optic axis in xz-plane at angle θ from z: (0, θ, 0)
• Quarter-wave plate at 45°: (π/4, π/2, 0) for optic axis in xy-plane

gyrotropic_tensor(d, od; axis=:z)

Constant gyrotropic (e.g. gyromagnetic Polder) tensor with diagonal d and antisymmetric imaginary off-diagonal ±i·od, with the gyration vector along axis. For axis=:z:

[ d      i·od   0
 −i·od   d      0
  0      0      1 ]

Hermitian (hence lossless) for real d, od. The sign/handedness follows the package exp(-iωt) convention (validated against the t_ps = -t_sp non-reciprocity relation).

polder_permeability(; f0, fm, linewidth=0.0, axis=:z)

Return a function f -> μ_tensor giving the gyromagnetic (Polder) permeability tensor of a saturated ferrite at frequency f:

\[μ(f) = 1 + \frac{f_0 f_m}{f_0^2 - f^2}, \qquad κ(f) = \frac{f f_m}{f_0^2 - f^2}\]

with the gyration axis axis. f0 is the ferromagnetic-resonance frequency (= γμ₀H₀/2π), f_m = γμ₀M_s/2π, and γ = g·e/2mₑ. A nonzero linewidth (ΔH, in the same units as f0) adds loss via f0 → f0 − i·linewidth/2.

Use as a layer permeability with Layer(...; mu = λ -> polder_permeability(...)(f(λ))), mapping wavelength to frequency as appropriate for your unit system.

rotate_dielectric_tensor(ε_diag, R)

Rotate a diagonal dielectric tensor from crystal frame to lab frame.

Given a diagonal tensor ε in the crystal's principal axis frame and a rotation matrix R, returns the rotated tensor: ε_lab = R * ε * R'

Arguments

• ε_diag: Diagonal dielectric tensor in crystal frame
• R: 3×3 rotation matrix from euler_rotation_matrix

Returns

Full 3×3 SMatrix (may have off-diagonal elements after rotation)

Sheet(σ)
Sheet(; xx, yy, xy=0, yx=0)
Sheet(material, d)
Sheet(nx, ny, d)

A zero-thickness 2D conductive sheet (e.g. a TMDC monolayer) for the transfer matrix. The internal representation is a callable λ -> SMatrix{2,2,ComplexF64} returning the SI sheet conductivity tensor [σxx σxy; σyx σyy] in Siemens.

• Sheet(σ) — scalar isotropic conductivity (σ::Number or σ(λ)->Number), stored as diag(σ, σ).
• Sheet(; xx, yy, xy=0, yx=0) — anisotropic tensor; each entry a Number or λ->Number.
• Sheet(material, d) — convert a refractive index to a sheet via σ = -i ω ε₀ d (n²-1) (isotropic). material is an n(λ) function (or a RefractiveMaterial when the RefractiveIndex extension is loaded); d is the effective thickness in μm.
• Sheet(nx, ny, d) — in-plane anisotropic index conversion (σxx from nx, σyy from ny).

The factor Z₀ = √(μ₀/ε₀) is applied exactly once, later, in sheet_matrix.

drude(ω_p, γ; ε_inf=1.0)

Return a closure λ -> n for the Drude free-carrier (metal) model:

\[ε(ω) = ε_\infty - \frac{ω_p^2}{ω^2 + iγω}, \qquad n = \sqrt{ε}.\]

Arguments

• ω_p: plasma energy (eV by default; a Unitful energy/wavenumber/frequency is normalized to eV when Unitful is loaded).
• γ: damping/collision rate (eV).
• ε_inf: high-frequency dielectric constant (dimensionless, default 1.0).

The returned closure takes a vacuum wavelength λ in μm and returns the complex refractive index n (with n″ > 0 for absorption under the package's exp(-iωt) convention).

Example

# Gold-like Drude metal, 50 nm film
au = drude(9.0, 0.07)        # ω_p = 9 eV, γ = 0.07 eV
layer = Layer(au, 0.05)

See also: lorentz, drude_lorentz.

drude_lorentz(ω_p, γ_D, oscillators; ε_inf=1.0)

Return a closure λ -> n for a combined Drude + Lorentz model (the standard form for real metals such as gold):

\[ε(ω) = ε_\infty - \frac{ω_p^2}{ω^2 + iγ_D ω} + \sum_j \frac{Δε_j\, ω_{0j}^2}{ω_{0j}^2 - ω^2 - iγ_j ω}.\]

Arguments

• ω_p, γ_D: Drude plasma energy and damping (eV).
• oscillators: iterable of Lorentz (ω_0, Δε, γ) 3-tuples (may be empty for a pure Drude metal).
• ε_inf: high-frequency dielectric constant (default 1.0).

The returned closure takes a vacuum wavelength λ in μm.

See also: drude, lorentz.

lorentz(ω_0, Δε, γ; ε_inf=1.0)
lorentz(oscillators; ε_inf=1.0)

Return a closure λ -> n for the Lorentz oscillator model:

\[ε(ω) = ε_\infty + \sum_j \frac{Δε_j\, ω_{0j}^2}{ω_{0j}^2 - ω^2 - iγ_j ω}, \qquad n = \sqrt{ε}.\]

Arguments

• ω_0: resonance energy (eV by default).
• Δε: dimensionless oscillator strength (the static contribution of the mode; ε(0) = ε_inf + Σ Δε_j).
• γ: damping (eV).
• ε_inf: high-frequency dielectric constant (default 1.0).
• oscillators: an iterable of (ω_0, Δε, γ) 3-tuples for a multi-oscillator model (e.g. several phonon modes).

The returned closure takes a vacuum wavelength λ in μm.

Examples

# Single oscillator
n = lorentz(2.0, 3.0, 0.05)                       # ω_0 = 2 eV

# Two oscillators on top of ε∞ = 2.5
n = lorentz([(2.0, 1.0, 0.05), (3.5, 0.4, 0.1)]; ε_inf=2.5)

See also: drude, drude_lorentz.

airy(n0, nf, ns, d, λ; θ=0.0)

Calculate the reflectance and transmittance of a single thin film using the exact Airy formula (multiple-beam interference).

Arguments

• n0: refractive index of incident medium (can be complex)
• nf: refractive index of the film (can be complex)
• ns: refractive index of the substrate (can be complex)
• d: film thickness (same units as λ)
• λ: wavelength (same units as d)
• θ=0.0: angle of incidence in radians

Returns

(Rs, Rp, Ts, Tp) - reflectance and transmittance for s and p polarizations.

Physics

The Airy formula accounts for all multiple reflections within the film:

\[r = \frac{r_{01} + r_{12} e^{2i\delta}}{1 + r_{01} r_{12} e^{2i\delta}}\]

where δ = 2π nf d cos(θf) / λ is the phase thickness and r₀₁, r₁₂ are the Fresnel reflection coefficients at the two interfaces.

This provides an exact analytical solution for validating numerical TMM results.

Example

# Quarter-wave anti-reflection coating
n_air, n_film, n_glass = 1.0, 1.38, 1.52
λ = 0.55  # μm
d = λ / (4 * n_film)  # quarter-wave thickness
Rs, Rp, Ts, Tp = airy(n_air, n_film, n_glass, d, λ)

See also: fresnel, fresnel_coefficients

dbr_reflectivity(no, ns, n1, n2, N)

Approximate the reflectivity of a DBR structure with originating medium with refractive index no, substrate with index ns, and alternating materials with indices n1 and n2 and number of repetitions N. The repeated pair of materials are assumed to have quarter-wave thickness $nd = \lambda / 4$, where $n$ is the refractive index, $d$ is the layer thickness, and $\lambda$ is the wavelength of the light.

Distributed Bragg reflector

fresnel(θ, n1, n2)

Calculate the reflectance for s-polarized and p-polarized light given the incidence angle θ (in radians) and indices of refraction of two media n1 and n2 at a plane interface.

Returns (Rs, Rp) where Rs is s-polarized reflectance and Rp is p-polarized reflectance.

The Fresnel equations for reflectance are:

\[R_s = \left| \frac{n_1 \cos\theta_i - n_2 \cos\theta_t}{n_1 \cos\theta_i + n_2 \cos\theta_t} \right|^2\]

\[R_p = \left| \frac{n_2 \cos\theta_i - n_1 \cos\theta_t}{n_2 \cos\theta_i + n_1 \cos\theta_t} \right|^2\]

where $\theta_t$ is the transmitted angle given by Snell's law: $n_1 \sin\theta_i = n_2 \sin\theta_t$.

Special cases:

• At grazing incidence (θ → π/2), returns (1.0, 1.0)
• For total internal reflection (when n1 > n2 and θ > θ_critical), returns (1.0, 1.0)

See also: fresnel_coefficients for complex amplitude coefficients.

Fresnel equations

fresnel_coefficients(θ, n1, n2)

Calculate the Fresnel reflection and transmission amplitude coefficients for s-polarized and p-polarized light at a plane interface.

Returns (rs, rp, ts, tp) where:

• rs, rp: reflection amplitude coefficients (complex)
• ts, tp: transmission amplitude coefficients (complex)

These are the amplitude (not intensity) coefficients, so reflectance R = |r|² and transmittance requires the full expression T = (n2 cosθt)/(n1 cosθi) |t|².

For total internal reflection, the reflection coefficients have unit magnitude with a phase shift (evanescent wave), and transmission coefficients are zero.

See also: fresnel for intensity reflectances.

stopband(n1, n2)

Calculate the frequency bandwidth Δf of the photonic stopband for a distributed bragg reflector (DBR) with two alternating materials of refractive indices n1 and n2.

\[ \frac{\Delta f_0}{f_0} = \frac{4}{\pi} \arcsin \left( \frac{n_2 - n_1}{n_2 + n_1} \right)\]

Distributed Bragg reflector