def dilution_fraction_from_species(n_e, impurity_densities, impurity_charges, main_ion_charge=1.0):
    r"""Main-ion dilution from an explicit impurity species list.

    $$f_\mathrm{dil} = 1 - \frac{n_\mathrm{main}}{n_e},\qquad
      n_\mathrm{main} = \frac{n_e - \sum_j Z_j n_j}{Z_m}
      \;\xrightarrow{Z_m = 1}\; f_\mathrm{dil} = \sum_j Z_j \frac{n_j}{n_e}$$

    Parameters
    ----------
    n_e : float or array-like
        Electron density, finite and positive [m^-3].
    impurity_densities : array-like
        Density of each impurity species along the last axis, non-negative [m^-3].
    impurity_charges : array-like
        Charge state of each species, positive [-].
    main_ion_charge : float, optional
        Main-ion charge $Z_m$, default 1 [-].

    Returns
    -------
    float or np.ndarray
        $f_\mathrm{dil}$ [-].

    Raises
    ------
    ValueError
        Invalid input, or impurities carrying more charge than there are
        electrons ($\sum_j Z_j n_j > n_e$).

    Convention
    ----------
    Preferred over :func:`dilution_fraction_from_zeff` whenever the species
    are resolved: it needs no effective charge and is exact for any mixture
    and any charge-state distribution (pass each charge state, or each
    element at its $\langle Z\rangle$ -- quasi-neutrality is linear in $Z$).

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 4.25.

    See Also
    --------
    main_ion_density_from_species
    """
    ne = _finite(n_e, "n_e")
    if np.any(ne <= 0.0):
        raise ValueError("n_e must be positive")
    main = np.asarray(main_ion_density_from_species(ne, impurity_densities, impurity_charges, main_ion_charge))
    return _out(1.0 - main / ne)

def dilution_fraction_from_zeff(z_eff, impurity_charge, main_ion_charge=1.0):
    r"""Main-ion dilution from Z_eff and one effective impurity charge.

    $$f_\mathrm{dil} = 1 - \frac{n_\mathrm{main}}{n_e}
      \;\xrightarrow{Z_m = 1}\; \frac{Z_\mathrm{eff} - 1}{Z_I - 1}$$

    Parameters
    ----------
    z_eff : float or array-like
        Plasma effective charge, in $[Z_m, Z_I]$ [-].
    impurity_charge : float or array-like
        Charge $Z_I$ of the one (effective) impurity, greater than $Z_m$ [-].
    main_ion_charge : float, optional
        Main-ion charge $Z_m$, default 1 [-].

    Returns
    -------
    float or np.ndarray
        $f_\mathrm{dil}$ [-].

    Raises
    ------
    ValueError
        As :func:`main_ion_fraction_from_zeff`.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 4.25.

    See Also
    --------
    main_ion_fraction_from_zeff
    dilution_fraction_from_species
    """
    return _out(1.0 - np.asarray(main_ion_fraction_from_zeff(z_eff, impurity_charge, main_ion_charge)))

def expand_effective_impurity(density, charge, weights, charges, mass=None, masses=None):
    r"""Species densities of a known relative composition behind a pseudo-impurity.

    $$n_i = \frac{n_{I,\mathrm{eff}} Z_{I,\mathrm{eff}}}{S_1}\,w_i,
      \qquad\text{requiring}\quad Z_{I,\mathrm{eff}} = \frac{S_2}{S_1}$$

    the inverse of :func:`reduce_impurity_mixture` once the relative
    composition $w_i$ is retained.

    Parameters
    ----------
    density : float or array-like
        Pseudo-impurity density $n_{I,\mathrm{eff}}$, non-negative, in any
        unit [m^-3].
    charge : float or array-like
        Pseudo-impurity charge $Z_{I,\mathrm{eff}}$, which must equal
        $S_2/S_1$ of the composition [-].
    weights : array-like
        Relative particle densities $w_i$ along the last axis, summing to one [-].
    charges : array-like
        Charge state $Z_i$ of each species, positive [-].
    mass : float or array-like, optional
        Pseudo-impurity mass $A_{I,\mathrm{eff}}$; when given with ``masses``
        it must equal $\bar A S_2/S_1^2$ [u].
    masses : array-like, optional
        Mass $A_i$ of each species, positive [u].

    Returns
    -------
    np.ndarray
        Density $n_i$ of each species, species on the last axis, in the unit of
        ``density`` [m^-3].

    Raises
    ------
    ValueError
        Invalid input, a negative density, ``mass`` without ``masses``, or a
        pseudo-impurity whose charge (or mass) is not the one this composition
        reduces to -- then no species densities reproduce both its charge and
        $Z^2$ moments.

    Assumptions
    -----------
    The relative composition is the one the pseudo-impurity was reduced
    from.  A pseudo-impurity alone does not determine a mixture; the
    weights are the extra information the round trip needs.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 4.25.

    See Also
    --------
    reduce_impurity_mixture
    """
    n_eff = _finite(density, "density")
    if np.any(n_eff < 0.0):
        raise ValueError("density must be non-negative")
    z_eff = _finite(charge, "charge")
    moments = impurity_mixture_moments(weights, charges, masses)
    s1 = np.asarray(moments.S1, dtype=float)
    s2 = np.asarray(moments.S2, dtype=float)
    if not np.allclose(z_eff, s2 / s1, rtol=1.0e-6, atol=0.0):
        raise ValueError(
            "charge is not the effective charge S2/S1 of this composition; "
            "no species densities keep both its charge and Z^2 moments"
        )
    if mass is not None and moments.A_bar is None:
        raise ValueError("mass can only be checked against the species masses; pass masses too")
    if mass is not None:
        expected = np.asarray(moments.A_bar) * s2 / s1**2
        if not np.allclose(_finite(mass, "mass"), expected, rtol=1.0e-6, atol=0.0):
            raise ValueError("mass is not the effective mass A_bar S2/S1^2 of this composition")
    w = _weights(weights)
    scale = n_eff * z_eff / s1
    result = np.asarray(scale)[..., None] * w if np.ndim(scale) else scale * w
    return np.asarray(result, dtype=float)

def impurity_mixture_moments(weights, charges, masses=None):
    r"""Charge and mass moments of a normalised impurity mixture.

    $$S_1 = \sum_i w_i Z_i,\qquad S_2 = \sum_i w_i Z_i^2,\qquad
      \bar A = \sum_i w_i A_i,\qquad \sum_i w_i = 1$$

    Parameters
    ----------
    weights : array-like
        Relative particle density $w_i$ of each impurity species along the
        last axis, non-negative and summing to one [-].
    charges : array-like
        Charge state $Z_i$ of each species, positive, broadcastable against
        ``weights`` [-].
    masses : array-like, optional
        Mass $A_i$ of each species, positive; omitted, $\bar A$ is ``None`` [u].

    Returns
    -------
    MixtureMoments
        ``S1`` and ``S2`` [-], and ``A_bar`` or ``None`` [u], each with the
        species axis removed [-].

    Raises
    ------
    ValueError
        Non-finite input, a negative weight, a non-positive charge or mass,
        or weights that do not sum to one.

    Convention
    ----------
    $Z_i$ is the ionic charge of the state the species is in, not the
    element's atomic number: carbon as C$^{4+}$ has $Z = 4$.  A partially
    ionised element contributes its charge-state-averaged $\langle Z\rangle$
    to $S_1$ and $\langle Z^2\rangle$ to $S_2$, not $\langle Z\rangle^2$.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 4.25.

    See Also
    --------
    solve_impurity_mixture_for_target_zeff
    reduce_impurity_mixture
    """
    w = _weights(weights)
    z = _species(charges, "charges", positive=True)
    w, z = np.broadcast_arrays(w, z)
    s1 = np.sum(w * z, axis=-1)
    s2 = np.sum(w * z**2, axis=-1)
    a_bar = None
    if masses is not None:
        a = _species(masses, "masses", positive=True)
        a_bar = _out(np.sum(np.broadcast_to(w, np.broadcast_shapes(w.shape, a.shape)) * a, axis=-1))
    return MixtureMoments(_out(s1), _out(s2), a_bar)

def main_ion_density_from_species(n_e, impurity_densities, impurity_charges, main_ion_charge=1.0):
    r"""Main-ion density that keeps a plasma with these impurities quasi-neutral.

    $$n_\mathrm{main} = \frac{n_e - \sum_j Z_j n_j}{Z_m}$$

    Parameters
    ----------
    n_e : float or array-like
        Electron density, finite and positive [m^-3].
    impurity_densities : array-like
        Density of each impurity species along the last axis, non-negative [m^-3].
    impurity_charges : array-like
        Charge state of each species, positive [-].
    main_ion_charge : float, optional
        Main-ion charge $Z_m$, default 1 [-].

    Returns
    -------
    float or np.ndarray
        $n_\mathrm{main}$ [m^-3].

    Raises
    ------
    ValueError
        Invalid input, or $\sum_j Z_j n_j > n_e$ (a negative main-ion density).

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 4.25.

    See Also
    --------
    dilution_fraction_from_species
    """
    zm = _main_charge(main_ion_charge)
    ne = _finite(n_e, "n_e")
    if np.any(ne <= 0.0):
        raise ValueError("n_e must be positive")
    n = _species(impurity_densities, "impurity_densities", positive=False)
    z = _species(impurity_charges, "impurity_charges", positive=True)
    n, z = np.broadcast_arrays(n, z)
    impurity_electrons = np.sum(n * z, axis=-1)
    main = (ne - impurity_electrons) / zm
    if np.any(main < -_RTOL * ne):
        raise ValueError("the impurities carry more charge than there are electrons")
    return _out(np.clip(main, 0.0, None))

def main_ion_density_from_zeff(n_e, z_eff, impurity_charge, main_ion_charge=1.0):
    r"""Main-ion density from n_e, Z_eff and one effective impurity charge.

    $$n_\mathrm{main} = n_e\,\frac{Z_I - Z_\mathrm{eff}}{Z_m\,(Z_I - Z_m)}$$

    Parameters
    ----------
    n_e : float or array-like
        Electron density, finite and positive [m^-3].
    z_eff : float or array-like
        Plasma effective charge, in $[Z_m, Z_I]$ [-].
    impurity_charge : float or array-like
        Charge $Z_I$ of the one (effective) impurity, greater than $Z_m$ [-].
    main_ion_charge : float, optional
        Main-ion charge $Z_m$, default 1 [-].

    Returns
    -------
    float or np.ndarray
        $n_\mathrm{main}$ [m^-3].

    Raises
    ------
    ValueError
        A non-positive ``n_e``, or as :func:`main_ion_fraction_from_zeff`.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 4.25.

    See Also
    --------
    main_ion_density_from_species
    """
    ne = _finite(n_e, "n_e")
    if np.any(ne <= 0.0):
        raise ValueError("n_e must be positive")
    return _out(ne * np.asarray(main_ion_fraction_from_zeff(z_eff, impurity_charge, main_ion_charge)))

def main_ion_fraction_from_zeff(z_eff, impurity_charge, main_ion_charge=1.0):
    r"""Main-ion fraction of the electron density from Z_eff and one effective impurity charge.

    $$f_\mathrm{main} = \frac{n_\mathrm{main}}{n_e}
      = \frac{Z_I - Z_\mathrm{eff}}{Z_m\,(Z_I - Z_m)}
      \;\xrightarrow{Z_m = 1}\; \frac{Z_I - Z_\mathrm{eff}}{Z_I - 1}$$

    Parameters
    ----------
    z_eff : float or array-like
        Plasma effective charge, in $[Z_m, Z_I]$ [-].
    impurity_charge : float or array-like
        Charge $Z_I$ of the one (effective) impurity, e.g.
        $Z_{I,\mathrm{eff}}(\rho)$ of a reduced mixture, greater than $Z_m$ [-].
    main_ion_charge : float, optional
        Main-ion charge $Z_m$, default 1 [-].

    Returns
    -------
    float or np.ndarray
        $f_\mathrm{main}$ [-].

    Raises
    ------
    ValueError
        Non-finite input, $Z_I \le Z_m$, or $Z_\mathrm{eff}$ outside
        $[Z_m, Z_I]$.

    Physical interpretation
    -----------------------
    A $Z_\mathrm{eff}$ scan is not a dilution scan: the same $Z_\mathrm{eff}$
    removes more main ions when the impurity charge is lower.  At
    $Z_\mathrm{eff} = 2$, carbon ($Z_I = 6$) leaves $f_\mathrm{main} = 0.8$,
    the VEST C/O mixture ($Z_{I,\mathrm{eff}} = 50/7$) leaves $36/43 \approx 0.837$.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 4.25.

    See Also
    --------
    dilution_fraction_from_zeff
    main_ion_density_from_species
    """
    zm = _main_charge(main_ion_charge)
    zi = _impurity_charge(impurity_charge, zm)
    zeff = _finite(z_eff, "z_eff")
    if np.any(zeff < zm - _RTOL) or np.any(zeff > zi + _RTOL * zi):
        raise ValueError("z_eff must lie in [Z_m, Z_I]")
    return _out(np.clip((zi - zeff) / (zm * (zi - zm)), 0.0, 1.0 / zm))

def reduce_impurity_mixture(densities, charges, masses=None):
    r"""The single pseudo-impurity that keeps a mixture's charge, Z^2 and mass moments.

    $$Z_{I,\mathrm{eff}} = \frac{\sum_i n_i Z_i^2}{\sum_i n_i Z_i},\qquad
      n_{I,\mathrm{eff}} = \frac{(\sum_i n_i Z_i)^2}{\sum_i n_i Z_i^2},\qquad
      A_{I,\mathrm{eff}} = \frac{\sum_i n_i A_i}{n_{I,\mathrm{eff}}}$$

    so that $n_{I,\mathrm{eff}} Z_{I,\mathrm{eff}} = \sum n_i Z_i$,
    $n_{I,\mathrm{eff}} Z_{I,\mathrm{eff}}^2 = \sum n_i Z_i^2$ and
    $n_{I,\mathrm{eff}} A_{I,\mathrm{eff}} = \sum n_i A_i$.

    Parameters
    ----------
    densities : array-like
        Density of each impurity species along the last axis, non-negative,
        in any one unit (absolute, or relative to $n_e$) [m^-3].
    charges : array-like
        Charge state $Z_i$ of each species, positive [-].
    masses : array-like, optional
        Mass $A_i$ of each species, positive; omitted, ``mass`` is ``None`` [u].

    Returns
    -------
    EffectiveImpurity
        ``charge`` $Z_{I,\mathrm{eff}}$ [-], ``density`` $n_{I,\mathrm{eff}}$
        in the unit of ``densities`` [m^-3], and ``mass`` $A_{I,\mathrm{eff}}$
        or ``None`` [u].

    Raises
    ------
    ValueError
        Non-finite or negative densities, non-positive charges or masses, or a
        mixture of zero total density.

    Physical interpretation
    -----------------------
    Quasi-neutrality and $Z_\mathrm{eff}$ see only the first two moments, so
    an interface that takes one impurity reproduces both exactly with this
    pseudo-species; the mass moment keeps the impurity mass density.  For the
    VEST C$^{6+}$:O$^{8+}$ = 1:1 mixture at $Z_\mathrm{eff} = 2$ this is
    $Z_{I,\mathrm{eff}} = 50/7$, $n_{I,\mathrm{eff}}/n_e = 49/2150$ and
    $A_{I,\mathrm{eff}} = \bar A S_2/S_1^2 \approx 14.29$.

    Limitations
    -----------
    Higher moments are not kept: a model sensitive to $Z^3$ or to the
    individual masses (impurity transport, collision operators resolving
    each species) sees a different plasma.  Prefer the explicit species list
    wherever the downstream code accepts one.

    References
    ----------
    .. [1] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011),
           Sec. 4.25.

    See Also
    --------
    expand_effective_impurity
    """
    n = _species(densities, "densities", positive=False)
    z = _species(charges, "charges", positive=True)
    n, z = np.broadcast_arrays(n, z)
    charge_density = np.sum(n * z, axis=-1)
    if np.any(charge_density <= 0.0):
        raise ValueError("a mixture of zero total density has no effective impurity")
    z2_density = np.sum(n * z**2, axis=-1)
    effective_charge = z2_density / charge_density
    effective_density = charge_density**2 / z2_density
    mass = None
    if masses is not None:
        a = _species(masses, "masses", positive=True)
        mass_density = np.sum(np.broadcast_to(n, np.broadcast_shapes(n.shape, a.shape)) * a, axis=-1)
        mass = _out(mass_density / effective_density)
    return EffectiveImpurity(_out(effective_charge), _out(effective_density), mass)