def binary_collision_energy_transfer_factor(m_1, m_2):
    r"""Largest fraction of a projectile's energy one elastic collision can give a target atom.

    $$\gamma = \frac{4m_1m_2}{(m_1 + m_2)^2}$$

    Parameters
    ----------
    m_1 : float or np.ndarray
        Projectile mass, positive [kg or u].
    m_2 : float or np.ndarray
        Target-atom mass, positive, same unit [kg or u].

    Returns
    -------
    float or np.ndarray
        $\gamma \in (0, 1]$, one for equal masses [-].

    Raises
    ------
    ValueError
        A mass is not positive.

    Physical interpretation
    -----------------------
    A head-on elastic collision transfers $\gamma E$. A light projectile
    (D on W: $\gamma \approx 0.043$) cannot give a heavy target atom much
    energy in one collision, which is why light ions sputter heavy targets
    only well above the surface binding energy -- the origin of the high
    sputtering threshold of tungsten under hydrogen.

    Assumptions
    -----------
    Two-body elastic collision, target initially at rest; exact kinematics.

    References
    ----------
    .. [1] W. Eckstein, *Computer Simulation of Ion-Solid Interactions*,
           Springer (1991), Ch. 2.
    """
    m_1 = _positive(m_1, "m_1")
    m_2 = _positive(m_2, "m_2")
    return _out(4.0 * m_1 * m_2 / (m_1 + m_2) ** 2)

def mean_reflected_energy_fraction(R_N, R_E):
    r"""Mean energy of a reflected particle, as a fraction of the incident energy.

    $$\frac{\langle E_\mathrm{refl}\rangle}{E_\mathrm{in}} = \frac{R_E}{R_N}$$

    Parameters
    ----------
    R_N : float or np.ndarray
        Particle reflection coefficient, in (0, 1] [-].
    R_E : float or np.ndarray
        Energy reflection coefficient, in [0, R_N] [-].

    Returns
    -------
    float or np.ndarray
        $\langle E_\mathrm{refl}\rangle/E_\mathrm{in}$ [-].

    Raises
    ------
    ValueError
        ``R_N`` outside (0, 1] or ``R_E`` outside [0, ``R_N``].

    Convention
    ----------
    $R_N$ counts reflected particles per incident particle, $R_E$ reflected
    energy per incident energy; they are different quantities with different
    data, and $R_E \le R_N$ because a reflected particle leaves with at most
    its incident energy.

    Physical interpretation
    -----------------------
    Particle balance and energy balance are not the same bookkeeping: a
    surface may return most particles ($R_N$ large) while keeping most of
    their energy ($R_E/R_N$ small), and it is $R_E$, not $R_N$, that enters
    the power the wall absorbs.

    References
    ----------
    .. [1] W. Eckstein, *Computer Simulation of Ion-Solid Interactions*,
           Springer (1991), Ch. 8.
    """
    R_N = np.asarray(R_N, dtype=float)
    R_E = np.asarray(R_E, dtype=float)
    if np.any(~(R_N > 0.0)) or np.any(R_N > 1.0):
        raise ValueError("R_N must lie in (0, 1]")
    if np.any(R_E < 0.0) or np.any(R_E > R_N):
        raise ValueError("R_E must lie in [0, R_N]: a reflected particle cannot leave with more than it brought")
    return _out(R_E / R_N)

def recycling_coefficient(Gamma_reflected, Gamma_reemitted, Gamma_incident):
    r"""Recycling coefficient: particles returned to the plasma per incident particle.

    $$R = \frac{\Gamma_\mathrm{refl} + \Gamma_\mathrm{re\text{-}em}}{\Gamma_\mathrm{in}},\qquad
      1 - R = \frac{\Gamma_\mathrm{retained}}{\Gamma_\mathrm{in}}$$

    Parameters
    ----------
    Gamma_reflected : float or np.ndarray
        Promptly reflected flux (fast atoms), non-negative [m^-2 s^-1].
    Gamma_reemitted : float or np.ndarray
        Flux re-emitted after implantation (thermal atoms and molecules,
        counted as atoms), non-negative [m^-2 s^-1].
    Gamma_incident : float or np.ndarray
        Incident ion plus atom flux, positive [m^-2 s^-1].

    Returns
    -------
    float or np.ndarray
        $R$; below one while the wall retains, above one while it releases a
        previous inventory [-].

    Raises
    ------
    ValueError
        A returned flux is negative or ``Gamma_incident`` is not positive.

    Convention
    ----------
    Counts atoms: a re-emitted D$_2$ molecule is two. Reflection is the
    prompt part only; recycling includes the delayed re-emission; retention
    is what neither returns. $R$ is not bounded by one: an outgassing wall
    returns more than it receives.

    References
    ----------
    .. [1] P. C. Stangeby, *The Plasma Boundary of Magnetic Fusion Devices*,
           IOP (2000), Ch. 3.
    """
    Gamma_reflected = np.asarray(Gamma_reflected, dtype=float)
    Gamma_reemitted = np.asarray(Gamma_reemitted, dtype=float)
    if np.any(~(Gamma_reflected >= 0.0)) or np.any(~(Gamma_reemitted >= 0.0)):
        raise ValueError("returned fluxes must be non-negative and finite")
    Gamma_incident = _positive(Gamma_incident, "Gamma_incident")
    return _out((Gamma_reflected + Gamma_reemitted) / Gamma_incident)

def sputtering_threshold_bohdansky(E_s, m_1, m_2):
    r"""Physical-sputtering threshold energy, Bohdansky's empirical fit.

    $$E_\mathrm{th} = \begin{cases}
      \dfrac{E_s}{\gamma(1 - \gamma)}, & m_1/m_2 \le 0.3 \\[1ex]
      8E_s\left(\dfrac{m_1}{m_2}\right)^{2/5}, & m_1/m_2 > 0.3
      \end{cases},\qquad \gamma = \frac{4m_1m_2}{(m_1 + m_2)^2}$$

    Parameters
    ----------
    E_s : float or np.ndarray
        Surface binding energy of the target (usually its sublimation
        energy), supplied by the caller, positive [eV].
    m_1 : float or np.ndarray
        Projectile mass, positive [kg or u].
    m_2 : float or np.ndarray
        Target-atom mass, positive, same unit [kg or u].

    Returns
    -------
    float or np.ndarray
        Threshold incident energy, normal incidence [eV].

    Raises
    ------
    ValueError
        A mass or ``E_s`` is not positive.

    Convention
    ----------
    A named empirical model (Bohdansky 1984), not a table: $E_s$ is an
    input and no material data are built in. Fits of Eckstein and others
    differ by tens of per cent near threshold; use tabulated thresholds for
    quantitative work. The breakpoint between the light-ion and heavy-ion
    branches is $m_1/m_2 = 0.3$, as every restatement of the fit gives it
    ([2]_ Eq. 12; [3]_); the two branches do not join exactly there: at
    $m_1/m_2 = 0.3$ the light branch gives $4.86E_s$ and the heavy one
    $4.94E_s$, a 1.75 % step. (An earlier revision of this function broke
    at 0.2, which put D on Be, He on N and He on O -- ratios between 0.2
    and 0.3 -- on the heavy branch, 6 % high.)

    Physical interpretation
    -----------------------
    A light projectile cannot sputter by a forward cascade: it is backscattered
    by a deeper target atom, keeping $(1 - \gamma)E$, and on its way out
    strikes a surface atom from below, giving it up to $\gamma(1 - \gamma)E$
    (``binary_collision_energy_transfer_factor`` for $\gamma$). That must
    exceed $E_s$: hence $E_s/\gamma(1 - \gamma)$, far above $E_s$ for D on W.

    References
    ----------
    .. [1] J. Bohdansky, J. Roth and H. L. Bay, J. Appl. Phys. 51 (1980) 2861;
           J. Bohdansky, Nucl. Instrum. Methods B 2 (1984) 587.
    .. [2] T. Nozawa, T. Kozasa and A. Habe, Astrophys. J. 648 (2006) 435,
           Eq. (12): the threshold with its $M_1/M_2 \le 0.3$ breakpoint.
    .. [3] W. Eckstein and R. Preuss, J. Nucl. Mater. 320 (2003) 209 (after
           C. Garcia-Rosales et al., J. Nucl. Mater. 218 (1994) 8): the same
           light/heavy split at $M_1/M_2 = 0.3$.
    """
    E_s = _positive(E_s, "E_s")
    m_1 = _positive(m_1, "m_1")
    m_2 = _positive(m_2, "m_2")
    ratio = m_1 / m_2
    gamma = 4.0 * m_1 * m_2 / (m_1 + m_2) ** 2
    with np.errstate(divide="ignore"):
        light = E_s / (gamma * (1.0 - gamma))  # gamma = 1 (equal masses) only reaches the heavy branch
    return _out(np.where(ratio <= 0.3, light, 8.0 * E_s * ratio**0.4))