Fields API

fields

Numerical transmission, layer-resolved absorption, and field profiles.

Where :mod:hyperbolic_optics.structure stops at reflection coefficients, this module reconstructs the field through the stack and derives power quantities from it directly (no closed-form transmission formulas). The tangential field F = [Ex, Ey, Hx, Hy] is continuous across interfaces and propagated by the stored layer matrices (Gᵢ = Mᵢ · G_{i+1}); the normal energy flux is S_z = ½ Re(Ex·Hy* − Ey·Hx*). For unit incident power this gives reflectance R = 1 − S_z(G₁)/S_z^inc, transmittance T = S_z(G_exit)/S_z^inc, and layer absorption Aᵢ = [S_z(Gᵢ) − S_z(G_{i+1})]/S_z^inc, which conserve energy exactly (R + T + ΣAᵢ = 1). For a single semi-infinite anisotropic exit there are no interior layers and :meth:FieldProfile.field_profile shows T absorbed with depth (decay length set by Im(kz)).

References:

• Passler, Jeannin & Paarmann, JOSA B 37, 1060 (2020); arXiv:2002.03832
• pyGTM: System.calculate_Efield / calculate_Poynting_Absorption_vs_z
• Passler & Paarmann, JOSA B 34, 2128-2139 (2017)

Classes:

• FieldProfile –

  Numerical transmission, layer-resolved absorption, and field profiles.

FieldProfile

FieldProfile(structure)

Numerical transmission, layer-resolved absorption, and field profiles.

Consumes an executed :class:~hyperbolic_optics.structure.Structure (mirrors how :class:~hyperbolic_optics.mueller.Mueller consumes one) and exposes the power quantities and field reconstruction described in the module docstring.

The engine is fully batched over the canonical [A, B, F] axes, so :meth:transmittance, :meth:reflectance and :meth:layer_absorption return arrays in the same presentation layout as structure.r_pp. :meth:field_profile adds a depth axis; it broadcasts too, but a full angle/frequency sweep × n_points × 6 field components is multi-gigabyte, so it is intended for the Simple scenario (or a single (angle, frequency) slice).

Examples:

>>> structure = Structure()
>>> structure.execute(payload)
>>> fp = FieldProfile(structure)
>>> fp.transmittance("p")          # power transmittance T
>>> fp.layer_absorption("p")       # per-interior-layer absorption
>>> fp.check_conservation("p")     # max |R + T + ΣA − 1|
>>> prof = fp.field_profile("p")   # dict: z, Ex..Hz, Sz, absorption_*

Cache the executed structure's transfer matrices and propagation constants.

Parameters:

• structure (Structure) –

  A :class:Structure on which execute has already run.

Raises:

• ValueError –

  If the structure has not been executed (no transfer matrix).

Methods:

• reflectance –

  Total power reflectance R for the given incident polarization.

• transmittance –

  Power transmittance T -- the flux crossing into the exit medium.

• transmission_coefficients –

  Amplitude transmission coefficients t_pp, t_ss, t_ps, t_sp.

• layer_absorption –

  Layer-resolved absorptance for each finite interior layer.

• summary –

  One-call R, T, per-layer absorption, total, and conservation residual.

• check_conservation –

  Return max|R + T + ΣAᵢ − 1| over the batch (≈ 0 when correct).

• polarization_resolved –

  Experimental. Split R and T into co- and cross-polarized channels.

• field_profile –

  Reconstruct E(z), H(z) and S_z(z) through the whole stack.

• stokes_from_field_profile –

  Polarization state of the transverse field vs depth (Stokes + ellipse).

Source code in hyperbolic_optics/fields.py

def __init__(self, structure: Structure) -> None:
    """Cache the executed structure's transfer matrices and propagation constants.

    Args:
        structure: A :class:`Structure` on which ``execute`` has already run.

    Raises:
        ValueError: If the structure has not been executed (no transfer matrix).
    """
    if structure.transfer_matrix is None:
        raise ValueError("Structure has not been executed; call structure.execute(payload).")
    self.structure = structure
    self.layers = structure.layers
    self.n_layers = len(self.layers)
    self.gamma = structure.transfer_matrix  # canonical [A, B, F, 4, 4]
    self.k_x = np.asarray(structure.k_x, dtype=np.complex128)  # [A, 1, 1]
    self.k_0 = np.asarray(structure.k_0, dtype=np.complex128)  # [1, 1, F]
    # Incident dynamical matrix A_inc = inv(prism transfer matrix); maps mode
    # amplitudes -> tangential field in the prism. [A, 1, 1, 4, 4].
    self.a_inc = np.linalg.inv(self.layers[0].matrix)

reflectance

reflectance(polarization='p')

Total power reflectance R for the given incident polarization.

R = 1 − S_z(G₁)/S_z^inc. For pure p/s incidence this equals |r_pp|²+|r_sp|² / |r_ss|²+|r_ps|² (a useful normalization check). Returned in the scenario's presentation shape (scalar for Simple).

Source code in hyperbolic_optics/fields.py

def reflectance(self, polarization: str | tuple[complex, complex] = "p") -> np.ndarray:
    """Total power reflectance ``R`` for the given incident polarization.

    ``R = 1 − S_z(G₁)/S_z^inc``. For pure p/s incidence this equals
    ``|r_pp|²+|r_sp|²`` / ``|r_ss|²+|r_ps|²`` (a useful normalization check).
    Returned in the scenario's presentation shape (scalar for ``Simple``).
    """
    _, _, _, s_inc, fields = self._solve(polarization)
    return present(1.0 - _poynting_z(fields[1]) / s_inc)

transmittance

transmittance(polarization='p')

Power transmittance T -- the flux crossing into the exit medium.

T = S_z(G_exit)/S_z^inc. For a transparent substrate this is the power that propagates away; for a lossy semi-infinite exit it is the power delivered into (and ultimately absorbed by) the bulk -- see :meth:field_profile for the depth distribution. Presentation shape.

Source code in hyperbolic_optics/fields.py

def transmittance(self, polarization: str | tuple[complex, complex] = "p") -> np.ndarray:
    """Power transmittance ``T`` -- the flux crossing into the exit medium.

    ``T = S_z(G_exit)/S_z^inc``. For a transparent substrate this is the power
    that propagates away; for a lossy semi-infinite exit it is the power
    delivered into (and ultimately absorbed by) the bulk -- see
    :meth:`field_profile` for the depth distribution. Presentation shape.
    """
    _, _, _, s_inc, fields = self._solve(polarization)
    return present(_poynting_z(fields[-1]) / s_inc)

transmission_coefficients

transmission_coefficients()

Amplitude transmission coefficients t_pp, t_ss, t_ps, t_sp.

These are the field-amplitude analogues of r_pp … -- exactly the coefficients pyGTM exposes from calculate_r_t and the quantities the field reconstruction is built on. They are the (s, p) forward transmission matrix T = M2⁻¹ (M2 = the forward sub-block of Γ), so c_exit = T · [a_s, a_p]. Naming follows the package convention t_{in→out} (matching r_*): e.g. t_ps is the s wave transmitted by incident p. They equal Structure.calculate_transmissivity exactly; transmittance then converts |amplitude|² to power with the correct Poynting/impedance weighting. Presentation shape.

Source code in hyperbolic_optics/fields.py

def transmission_coefficients(self) -> dict[str, np.ndarray]:
    """Amplitude transmission coefficients ``t_pp, t_ss, t_ps, t_sp``.

    These are the field-amplitude analogues of ``r_pp …`` -- exactly the
    coefficients pyGTM exposes from ``calculate_r_t`` and the quantities the
    field reconstruction is built on. They are the ``(s, p)`` forward
    transmission matrix ``T = M2⁻¹`` (``M2`` = the forward sub-block of ``Γ``),
    so ``c_exit = T · [a_s, a_p]``. Naming follows the package convention
    ``t_{in→out}`` (matching ``r_*``): e.g. ``t_ps`` is the s wave transmitted
    by incident p. They equal ``Structure.calculate_transmissivity`` exactly;
    ``transmittance`` then converts |amplitude|² to *power* with the correct
    Poynting/impedance weighting. Presentation shape.
    """
    g = self.gamma
    row0 = np.stack([g[..., _S_FWD, _S_FWD], g[..., _S_FWD, _P_FWD]], axis=-1)
    row1 = np.stack([g[..., _P_FWD, _S_FWD], g[..., _P_FWD, _P_FWD]], axis=-1)
    m2 = np.stack([row0, row1], axis=-2)  # [A, B, F, 2, 2]; rows/cols = (s, p)
    t = np.linalg.inv(m2)  # T[out, in]: row = transmitted (s, p), col = incident (s, p)
    return {
        "t_ss": present(t[..., 0, 0]),
        "t_ps": present(t[..., 0, 1]),  # in p -> out s
        "t_sp": present(t[..., 1, 0]),  # in s -> out p
        "t_pp": present(t[..., 1, 1]),
    }

layer_absorption

layer_absorption(polarization='p')

Layer-resolved absorptance for each finite interior layer.

Returns:

• list[dict[str, Any]] –

  A list (top→bottom) of dicts {"index", "type", "absorptance"}, where

• list[dict[str, Any]] –

  absorptance is the fraction of incident power dissipated in that

• list[dict[str, Any]] –

  layer, in presentation shape. Empty for a prism + semi-infinite-exit

• list[dict[str, Any]] –

  stack (no interior layers -- all absorbed power shows up as T into

• list[dict[str, Any]] –

  the bulk; use :meth:field_profile to resolve it with depth).

Source code in hyperbolic_optics/fields.py

def layer_absorption(
    self, polarization: str | tuple[complex, complex] = "p"
) -> list[dict[str, Any]]:
    """Layer-resolved absorptance for each finite interior layer.

    Returns:
        A list (top→bottom) of dicts ``{"index", "type", "absorptance"}``, where
        ``absorptance`` is the fraction of incident power dissipated in that
        layer, in presentation shape. Empty for a prism + semi-infinite-exit
        stack (no interior layers -- all absorbed power shows up as ``T`` into
        the bulk; use :meth:`field_profile` to resolve it with depth).
    """
    _, _, _, s_inc, fields = self._solve(polarization)
    return [
        {"index": i, "type": layer_type, "absorptance": present(absorbed)}
        for i, layer_type, absorbed in self._layer_absorption_canonical(s_inc, fields)
    ]

summary

summary(polarization='p')

One-call R, T, per-layer absorption, total, and conservation residual.

conservation_residual = max|R + T + ΣAᵢ − 1| over the batch; it should be ~machine-epsilon for a correct, energy-conserving result.

Source code in hyperbolic_optics/fields.py

def summary(self, polarization: str | tuple[complex, complex] = "p") -> dict[str, Any]:
    """One-call ``R``, ``T``, per-layer absorption, total, and conservation residual.

    ``conservation_residual = max|R + T + ΣAᵢ − 1|`` over the batch; it should
    be ~machine-epsilon for a correct, energy-conserving result.
    """
    _, _, _, s_inc, fields = self._solve(polarization)
    r = 1.0 - _poynting_z(fields[1]) / s_inc
    t = _poynting_z(fields[-1]) / s_inc
    per_layer = self._layer_absorption_canonical(s_inc, fields)
    total_abs = sum((a for _, _, a in per_layer), start=np.zeros_like(r))
    residual = np.abs(r + t + total_abs - 1.0)
    return {
        "R": present(r),
        "T": present(t),
        "layers": [
            {"index": i, "type": ty, "absorptance": present(a)} for i, ty, a in per_layer
        ],
        "total_absorption": present(total_abs),
        "conservation_residual": float(np.max(residual)),
    }

check_conservation

check_conservation(polarization='p')

Return max|R + T + ΣAᵢ − 1| over the batch (≈ 0 when correct).

Source code in hyperbolic_optics/fields.py

def check_conservation(self, polarization: str | tuple[complex, complex] = "p") -> float:
    """Return ``max|R + T + ΣAᵢ − 1|`` over the batch (≈ 0 when correct)."""
    return self.summary(polarization)["conservation_residual"]

polarization_resolved

polarization_resolved(polarization='p')

Experimental. Split R and T into co- and cross-polarized channels.

For a pure "p" or "s" incidence, the reflected/transmitted power is decomposed into the co-polarized channel (same polarization as incident) and the cross-polarized channel (the polarization-converted light — the power analogue of r_ps / t_ps). The conversion fractions cross / (co + cross) quantify how much light changed polarization.

Reflection is split exactly (R_co = |r_co|², R_cross = |r_cross|², reflection being back into the same prism). Transmission is split via the per-mode flux of the two forward exit waves.

Experimental / caveat: the transmitted s/p split is rigorous only for an isotropic exit medium, where the two exit eigenmodes are clean s and p and their fluxes add (T_co + T_cross = T). For an anisotropic exit the eigenmodes are elliptical, so the split is an eigenmode-resolved approximation and may not sum exactly to the total transmittance.

Parameters:

• polarization (str, default: 'p' ) –

  "p" or "s" (co/cross is undefined for a mixed Jones state).

Returns:

• dict[str, ndarray] –

  Dict (presentation shape) with R_co, R_cross, T_co, T_cross and

• dict[str, ndarray] –

  conversion_reflection, conversion_transmission.

Raises:

• ValueError –

  If polarization is not "p" or "s".

Source code in hyperbolic_optics/fields.py

def polarization_resolved(self, polarization: str = "p") -> dict[str, np.ndarray]:
    """**Experimental.** Split R and T into co- and cross-polarized channels.

    For a pure ``"p"`` or ``"s"`` incidence, the reflected/transmitted power is
    decomposed into the *co-polarized* channel (same polarization as incident)
    and the *cross-polarized* channel (the polarization-converted light — the
    power analogue of ``r_ps`` / ``t_ps``). The conversion fractions
    ``cross / (co + cross)`` quantify how much light changed polarization.

    Reflection is split exactly (``R_co = |r_co|²``, ``R_cross = |r_cross|²``,
    reflection being back into the same prism). Transmission is split via the
    per-mode flux of the two forward exit waves.

    Experimental / caveat: the transmitted s/p split is rigorous only for an
    **isotropic exit** medium, where the two exit eigenmodes are clean s and p
    and their fluxes add (``T_co + T_cross = T``). For an anisotropic exit the
    eigenmodes are elliptical, so the split is an eigenmode-resolved
    approximation and may not sum exactly to the total transmittance.

    Args:
        polarization: ``"p"`` or ``"s"`` (co/cross is undefined for a mixed
            Jones state).

    Returns:
        Dict (presentation shape) with ``R_co, R_cross, T_co, T_cross`` and
        ``conversion_reflection``, ``conversion_transmission``.

    Raises:
        ValueError: If ``polarization`` is not ``"p"`` or ``"s"``.
    """
    key = polarization.lower() if isinstance(polarization, str) else None
    if key not in ("p", "s"):
        raise ValueError("polarization_resolved requires pure 'p' or 's' incidence.")
    _, _, c_exit, s_inc, _ = self._solve(polarization)

    # Reflected amplitudes are the backward slots of c_inc = Gamma @ c_exit.
    c_inc = _matvec(self.gamma, c_exit)
    refl_s, refl_p = c_inc[..., _S_BWD], c_inc[..., _P_BWD]

    # Per-mode transmitted flux (exit modes: index 0 = s-like, 1 = p-like).
    w_full, _, amps = self._exit_transmitted_modes(c_exit)
    tangential = w_full[..., [0, 1, 3, 4], :]  # [..., (Ex,Ey,Hx,Hy), 2]

    def mode_transmittance(mode: int) -> np.ndarray:
        field = amps[..., mode][..., np.newaxis] * tangential[..., :, mode]
        return _poynting_z(field) / s_inc

    t_s, t_p = mode_transmittance(0), mode_transmittance(1)

    if key == "p":
        r_co, r_cross, t_co, t_cross = refl_p, refl_s, t_p, t_s
    else:
        r_co, r_cross, t_co, t_cross = refl_s, refl_p, t_s, t_p

    r_co_pow = np.abs(r_co) ** 2
    r_cross_pow = np.abs(r_cross) ** 2

    def _fraction(cross: np.ndarray, co: np.ndarray) -> np.ndarray:
        total = co + cross
        return np.where(np.abs(total) > 1e-30, cross / total, 0.0)

    return {
        "R_co": present(r_co_pow),
        "R_cross": present(r_cross_pow),
        "T_co": present(t_co),
        "T_cross": present(t_cross),
        "conversion_reflection": present(_fraction(r_cross_pow, r_co_pow)),
        "conversion_transmission": present(_fraction(t_cross, t_co)),
    }

field_profile

field_profile(polarization='p', n_points=200, semi_inf_thickness=None)

Reconstruct E(z), H(z) and S_z(z) through the whole stack.

Parameters:

• polarization (str | tuple[complex, complex], default: 'p' ) –

  "p", "s" or a complex (a_s, a_p) Jones pair.

• n_points (int, default: 200 ) –

  Depth samples per layer (interior layers and the exit window).

• semi_inf_thickness (float | None, default: None ) –

  Display window for the semi-infinite exit, in µm. None auto-selects ~5 decay lengths of the slowest transmitted mode (a few wavelengths if lossless).

Returns:

• dict[str, ndarray] –

  Dict with z (µm, measured from the first interface) and the squeezed

• dict[str, ndarray] –

  field arrays Ex, Ey, Ez, Hx, Hy, Hz (complex) plus Sz (real,

• dict[str, ndarray] –

  normalized to incident flux). Also absorption_cumulative =

• dict[str, ndarray] –

  (Sz[0] − Sz)/1 running from the surface and layer_boundaries (µm).

Note

Intended for Simple / single-point use (the depth axis multiplies the batch size); see the class docstring.

Raises:

• ValueError –

  If a layer thickness is being swept (canonical T axis size > 1) — the depth grid of the swept layer would itself change along T, so a field-vs-depth profile is not well defined.

Source code in hyperbolic_optics/fields.py

def field_profile(
    self,
    polarization: str | tuple[complex, complex] = "p",
    n_points: int = 200,
    semi_inf_thickness: float | None = None,
) -> dict[str, np.ndarray]:
    """Reconstruct ``E(z), H(z)`` and ``S_z(z)`` through the whole stack.

    Args:
        polarization: ``"p"``, ``"s"`` or a complex ``(a_s, a_p)`` Jones pair.
        n_points: Depth samples per layer (interior layers and the exit window).
        semi_inf_thickness: Display window for the semi-infinite exit, in **µm**.
            ``None`` auto-selects ~5 decay lengths of the slowest transmitted
            mode (a few wavelengths if lossless).

    Returns:
        Dict with ``z`` (µm, measured from the first interface) and the squeezed
        field arrays ``Ex, Ey, Ez, Hx, Hy, Hz`` (complex) plus ``Sz`` (real,
        normalized to incident flux). Also ``absorption_cumulative`` =
        ``(Sz[0] − Sz)/1`` running from the surface and ``layer_boundaries`` (µm).

    Note:
        Intended for ``Simple`` / single-point use (the depth axis multiplies
        the batch size); see the class docstring.

    Raises:
        ValueError: If a layer thickness is being swept (canonical ``T`` axis
            size > 1) — the depth grid of the swept layer would itself change
            along ``T``, so a field-vs-depth profile is not well defined.
    """
    if self.gamma.shape[T] > 1:
        raise ValueError(
            "field_profile is undefined while sweeping a layer thickness "
            f"(T={self.gamma.shape[T]}); the depth grid changes along T. Profile a "
            "single thickness, or use the power quantities (reflectance/"
            "transmittance/layer_absorption), which support the T axis."
        )
    _, _, c_exit, s_inc, fields = self._solve(polarization)
    cm_to_um = 1.0e4

    z_segments: list[np.ndarray] = []
    comp_segments: list[np.ndarray] = []  # each [..., Nz, 6]
    boundaries = [0.0]
    z_offset = 0.0  # cm

    # Interior finite layers: amplitudes from V^{-1} G_i over all four modes.
    for i in range(1, self.n_layers - 1):
        layer = self.layers[i]
        # T == 1 here (guarded above), so a swept layer's [1,1,1,1] thickness
        # still reduces to a single scalar depth.
        thickness = float(np.asarray(layer.thickness).reshape(-1)[0])
        v_tan, kz = layer.profile.tangential_modes()
        w_full, _ = layer.profile.full_modes()
        amps = np.linalg.solve(v_tan, fields[i][..., np.newaxis])[..., 0]  # [..., 4]
        depth = np.linspace(0.0, thickness, n_points)
        comp_segments.append(self._reconstruct(w_full, kz, amps, depth))
        z_segments.append((z_offset + depth) * cm_to_um)
        z_offset += thickness
        boundaries.append(z_offset * cm_to_um)

    # Semi-infinite exit: two transmitted modes over a display window.
    w_full, kz_exit, amps = self._exit_transmitted_modes(c_exit)
    if semi_inf_thickness is None:
        exit_len = self._auto_semi_infinite_thickness(kz_exit)
    else:
        exit_len = float(semi_inf_thickness) / cm_to_um
    depth = np.linspace(0.0, exit_len, n_points)
    comp_segments.append(self._reconstruct(w_full, kz_exit, amps, depth))
    z_segments.append((z_offset + depth) * cm_to_um)
    boundaries.append((z_offset + exit_len) * cm_to_um)

    components = np.concatenate(comp_segments, axis=-2)  # [..., Nz_total, 6]
    z = np.concatenate(z_segments, axis=-1)
    sz = _poynting_z(components[..., [0, 1, 3, 4]]) / s_inc[..., np.newaxis]

    labels = ("Ex", "Ey", "Ez", "Hx", "Hy", "Hz")
    result: dict[str, np.ndarray] = {"z": z, "layer_boundaries": np.asarray(boundaries)}
    for j, label in enumerate(labels):
        result[label] = np.squeeze(components[..., j])
    sz_sq = np.squeeze(sz)
    result["Sz"] = sz_sq
    result["absorption_cumulative"] = np.take(sz_sq, 0, axis=-1)[..., np.newaxis] - sz_sq
    return result

stokes_from_field_profile

stokes_from_field_profile(polarization='p', n_points=200, semi_inf_thickness=None)

Polarization state of the transverse field vs depth (Stokes + ellipse).

Reconstructs the field with :meth:field_profile and forms the Stokes parameters of the transverse pair (Ex, Ey) at every depth, so you can watch the polarization ellipse evolve as the wave crosses a birefringent layer. Ex is the in-plane (p-like) component, Ey the out-of-plane (s-like) one; the convention is e^{-iωt}.

Parameters:

• polarization (str | tuple[complex, complex], default: 'p' ) –

  "p", "s" or a complex (a_s, a_p) Jones pair.

• n_points (int, default: 200 ) –

  Depth samples per layer (passed to :meth:field_profile).

• semi_inf_thickness (float | None, default: None ) –

  Semi-infinite exit display window in µm (see :meth:field_profile).

Returns:

• dict[str, ndarray] –

  Dict with z (µm) and layer_boundaries plus the Stokes profiles

• dict[str, ndarray] –

  S0 (intensity) and S1, S2, S3, and the ellipse parameters

• dict[str, ndarray] –

  azimuth (ψ, rad) and ellipticity (χ, rad). A single coherent

• dict[str, ndarray] –

  field is fully polarized, so the degree of polarization is 1 by

• dict[str, ndarray] –

  construction and is not returned.

Note

Like :meth:field_profile, this is undefined while sweeping a layer thickness (raises if the canonical T axis size > 1).

Source code in hyperbolic_optics/fields.py

def stokes_from_field_profile(
    self,
    polarization: str | tuple[complex, complex] = "p",
    n_points: int = 200,
    semi_inf_thickness: float | None = None,
) -> dict[str, np.ndarray]:
    """Polarization state of the transverse field vs depth (Stokes + ellipse).

    Reconstructs the field with :meth:`field_profile` and forms the Stokes
    parameters of the transverse pair ``(Ex, Ey)`` at every depth, so you can
    watch the polarization ellipse evolve as the wave crosses a birefringent
    layer. ``Ex`` is the in-plane (p-like) component, ``Ey`` the out-of-plane
    (s-like) one; the convention is ``e^{-iωt}``.

    Args:
        polarization: ``"p"``, ``"s"`` or a complex ``(a_s, a_p)`` Jones pair.
        n_points: Depth samples per layer (passed to :meth:`field_profile`).
        semi_inf_thickness: Semi-infinite exit display window in µm (see
            :meth:`field_profile`).

    Returns:
        Dict with ``z`` (µm) and ``layer_boundaries`` plus the Stokes profiles
        ``S0`` (intensity) and ``S1, S2, S3``, and the ellipse parameters
        ``azimuth`` (ψ, rad) and ``ellipticity`` (χ, rad). A single coherent
        field is fully polarized, so the degree of polarization is 1 by
        construction and is not returned.

    Note:
        Like :meth:`field_profile`, this is undefined while sweeping a layer
        thickness (raises if the canonical ``T`` axis size > 1).
    """
    profile = self.field_profile(polarization, n_points, semi_inf_thickness)
    ex, ey = profile["Ex"], profile["Ey"]
    ex_sq, ey_sq = np.abs(ex) ** 2, np.abs(ey) ** 2
    s0 = ex_sq + ey_sq
    s1 = ex_sq - ey_sq
    s2 = 2.0 * np.real(ex * np.conj(ey))
    s3 = -2.0 * np.imag(ex * np.conj(ey))
    return {
        "z": profile["z"],
        "layer_boundaries": profile["layer_boundaries"],
        "S0": s0,
        "S1": s1,
        "S2": s2,
        "S3": s3,
        "azimuth": 0.5 * np.arctan2(s2, s1),
        "ellipticity": 0.5 * np.arctan2(s3, np.sqrt(s1**2 + s2**2)),
    }