def cma_coordinates(omega, n_e, B):
    r"""CMA-diagram coordinates of a cold electron plasma.

    $$X = \frac{\omega_{pe}^2}{\omega^2},\qquad Y = \frac{|\Omega_e|}{\omega}$$

    Parameters
    ----------
    omega : float or np.ndarray
        Wave angular frequency, positive [rad/s].
    n_e : float or np.ndarray
        Electron density, non-negative [m^-3].
    B : float or np.ndarray
        Magnetic-field magnitude [T].

    Returns
    -------
    tuple of float or np.ndarray
        ``(X, Y)``, both non-negative [-].

    Raises
    ------
    ValueError
        A non-positive ``omega`` or a negative density.

    Convention
    ----------
    $Y$ is the positive magnitude $|\Omega_e|/\omega$, the usual CMA axis;
    ``stix_parameters`` uses the signed $\Omega_e$. For electrons only,
    $P = 1 - X$, $R = 1 - X/(1 - Y)$, $L = 1 - X/(1 + Y)$ and
    $S = 1 - X/(1 - Y^2)$.

    Physical interpretation
    -----------------------
    A point $(X, Y)$ is a local plasma at one frequency: moving along a
    profile at fixed $\omega$ traces a curve in the CMA plane, and each
    boundary it crosses is a cutoff or resonance of some mode.

    References
    ----------
    .. [1] T. H. Stix, *Waves in Plasmas*, AIP (1992), Sec. 2-2.
    """
    omega = np.asarray(omega, dtype=float)
    if np.any(~(omega > 0.0)):
        raise ValueError("omega must be positive")
    X = np.asarray(plasma_frequency(n_e, -QE, ME), dtype=float) ** 2 / omega**2
    Y = np.abs(np.asarray(gyrofrequency(-QE, ME, B), dtype=float)) / omega
    return _out(X), _out(Y)

def dielectric_tensor(R, L, P):
    r"""Cold-plasma dielectric tensor in the frame with $\mathbf B_0 \parallel \hat z$.

    $$\boldsymbol\epsilon = \begin{pmatrix} S & -iD & 0\\ iD & S & 0\\ 0 & 0 & P\end{pmatrix}$$

    Parameters
    ----------
    R : float or np.ndarray
        Stix $R$ [-].
    L : float or np.ndarray
        Stix $L$ [-].
    P : float or np.ndarray
        Stix $P$ [-].

    Returns
    -------
    np.ndarray
        Complex relative permittivity, shape ``broadcast + (3, 3)`` [-].

    Convention
    ----------
    Fields $\propto e^{i(\mathbf k\cdot\mathbf x - \omega t)}$, $\hat z$ along
    $\mathbf B_0$; with the opposite time convention the off-diagonal signs
    flip. $S = (R + L)/2$, $D = (R - L)/2$ as ``stix_parameters`` returns them.

    Physical interpretation
    -----------------------
    Hermitian for real Stix parameters: a collisionless cold plasma is
    lossless. The antisymmetric $\pm iD$ is the gyrotropy that makes the
    circular polarisations $R$ and $L$ the eigenmodes along $\mathbf B_0$.

    References
    ----------
    .. [1] T. H. Stix, *Waves in Plasmas*, AIP (1992), Sec. 1-2.
    """
    R, L, P = np.broadcast_arrays(*(np.asarray(v, dtype=float) for v in (R, L, P)))
    S, D = (R + L) / 2.0, (R - L) / 2.0
    eps = np.zeros(R.shape + (3, 3), dtype=complex)
    eps[..., 0, 0] = eps[..., 1, 1] = S
    eps[..., 0, 1] = -1j * D
    eps[..., 1, 0] = 1j * D
    eps[..., 2, 2] = P
    return eps

def perpendicular_refractive_index_squared(R, L, P):
    r"""Ordinary and extraordinary modes at perpendicular propagation.

    $$n_O^2 = P,\qquad n_X^2 = \frac{RL}{S}$$

    Parameters
    ----------
    R : float or np.ndarray
        Stix $R$ [-].
    L : float or np.ndarray
        Stix $L$ [-].
    P : float or np.ndarray
        Stix $P$ [-].

    Returns
    -------
    tuple of float or np.ndarray
        ``(n2_O, n2_X)``; $n_X^2$ is infinite at $S = 0$, the hybrid
        resonances [-].

    Convention
    ----------
    $\theta = \pi/2$: the O mode has $\mathbf E \parallel \mathbf B_0$ and
    feels only $P$; the X mode has $\mathbf E \perp \mathbf B_0$. These are
    the two roots of ``cold_plasma_refractive_index_squared`` at
    $\theta = \pi/2$, named by polarisation.

    Physical interpretation
    -----------------------
    O-mode cutoff: $P = 0$, $\omega = \omega_{pe}$ (with cold ions). X-mode
    cutoffs: $R = 0$ and $L = 0$; its resonance $S = 0$ is the upper (and
    lower) hybrid layer, reached from the high-field side or by tunnelling.

    References
    ----------
    .. [1] T. H. Stix, *Waves in Plasmas*, AIP (1992), Sec. 1-4.
    """
    R, L, P = (np.asarray(v, dtype=float) for v in (R, L, P))
    S = (R + L) / 2.0
    with np.errstate(divide="ignore", invalid="ignore"):
        n2_X = np.where(S != 0.0, R * L / S, np.inf)
    return _out(P), _out(n2_X)

def plasma_frequency(n, q, m):
    r"""Plasma frequency of one species.

    $$\omega_{ps} = \sqrt{\frac{n_sq_s^2}{\epsilon_0m_s}}$$

    Parameters
    ----------
    n : float or np.ndarray
        Density, non-negative [m^-3].
    q : float or np.ndarray
        Charge, either sign [C].
    m : float or np.ndarray
        Mass, positive [kg].

    Returns
    -------
    float or np.ndarray
        $\omega_{ps}$, non-negative [rad/s].

    Raises
    ------
    ValueError
        A negative density, a non-positive mass, or a non-finite input.

    Convention
    ----------
    An angular frequency: divide by $2\pi$ for hertz ($f_{pe} \approx
    8.98\sqrt{n_e}$ Hz with $n_e$ in m$^{-3}$). The charge enters squared,
    so its sign does not matter here; it does in ``stix_parameters``.

    Physical interpretation
    -----------------------
    The frequency at which the species oscillates about a neutralising
    background when displaced; waves below $\omega_{pe}$ without a magnetic
    field cannot propagate.

    References
    ----------
    .. [1] T. H. Stix, *Waves in Plasmas*, AIP (1992), Sec. 1-2.
    """
    n = np.asarray(n, dtype=float)
    q = np.asarray(q, dtype=float)
    m = np.asarray(m, dtype=float)
    if not (np.all(np.isfinite(n)) and np.all(np.isfinite(q)) and np.all(np.isfinite(m))):
        raise ValueError("n, q and m must be finite")
    if np.any(n < 0.0) or np.any(m <= 0.0):
        raise ValueError("n must be non-negative and m positive")
    return _out(np.sqrt(n * q**2 / (EPS0 * m)))

def propagation_regime(n2, *, atol=0.0):
    r"""Classify a refractive index squared: propagating, evanescent, cutoff or resonance.

    $$n^2 > 0: \text{propagating},\quad n^2 < 0: \text{evanescent},\quad
      n^2 = 0: \text{cutoff},\quad |n^2| \to \infty: \text{resonance}$$

    Parameters
    ----------
    n2 : float or np.ndarray
        Refractive index squared [-].
    atol : float
        Half-width of the band $|n^2| \le$ ``atol`` counted as a cutoff;
        zero by default [-].

    Returns
    -------
    str or np.ndarray
        ``"propagating"``, ``"evanescent"``, ``"cutoff"`` or ``"resonance"``
        per element [-].

    Raises
    ------
    ValueError
        A negative ``atol``, or a NaN in ``n2``.

    Convention
    ----------
    A resonance is $\pm\infty$ (as ``cold_plasma_refractive_index_squared``
    returns at $A = 0$); a large finite $n^2$ near it is still
    "propagating" or "evanescent". Cutoff and resonance are layers of zero
    width, so on a sampled profile they appear only through ``atol`` or as
    sign changes between samples.

    References
    ----------
    .. [1] T. H. Stix, *Waves in Plasmas*, AIP (1992), Sec. 1-3.
    """
    if not atol >= 0.0:
        raise ValueError("atol must be non-negative")
    n2 = np.asarray(n2, dtype=float)
    if np.any(np.isnan(n2)):
        raise ValueError("n2 must not be NaN")
    result = np.where(np.isinf(n2), "resonance",
                      np.where(np.abs(n2) <= atol, "cutoff", np.where(n2 > 0.0, "propagating", "evanescent")))
    return str(result) if result.ndim == 0 else result