def equilibrium_particle_field(equilibrium, *, derivative_step=1e-4):
    """Adapt an existing VAFT equilibrium field interpolator to particle tracing.

    Parameters
    ----------
    equilibrium : EquilibriumData
        Validated equilibrium with explicit resolved COCOS convention [-].
    derivative_step : float, optional
        Cartesian field derivative step, a numerical convenience [m].

    Returns
    -------
    ParticleField
        Cartesian field restricted to the equilibrium map domain [-].

    Processing steps
    ----------------
    Reuse the equilibrium interpolator and its COCOS sign/flux handling. Convert
    its (B_R,B_Z,B_phi) components to Cartesian order and enforce map bounds.

    Applicability
    -------------
    Machine-independent. A resolved static equilibrium snapshot.

    Limitations
    -----------
    The underlying interpolator clips F beyond the profile range; this adapter
    does not validate the equilibrium or invent a vector-potential gauge.

    Provenance
    ----------
    .. [1] vaft.process.equilibrium.make_equilibrium_field_interpolator.
    """
    from vaft.process.equilibrium import make_equilibrium_field_interpolator

    if equilibrium.convention.cocos is None:
        raise ValueError("particle equilibrium adapter requires resolved COCOS")
    native = make_equilibrium_field_interpolator(equilibrium)
    r, z = equilibrium.r, equilibrium.z
    if r is None or z is None:
        raise ValueError("equilibrium map axes are required")

    def B(x):
        R = np.hypot(x[0], x[1])
        br, bz, bt = native(R, x[2])
        return np.array([(br * x[0] - bt * x[1]) / R, (br * x[1] + bt * x[0]) / R, bz])

    return ParticleField(
        B,
        derivative_step=derivative_step,
        contains=lambda x: r[0] < np.hypot(x[0], x[1]) < r[-1] and z[0] < x[2] < z[-1],
        name="VAFT equilibrium particle field",
        provenance={
            "cocos": equilibrium.convention.cocos,
            "field_source": "make_equilibrium_field_interpolator; F edge clipping retained",
        },
    )

def mirror_particle_field(B0=1.0, length=10.0):
    """Construct an exactly solenoidal axisymmetric parabolic mirror field.

    Parameters
    ----------
    B0 : float, optional
        Minimum on-axis field, a benchmark default [T].
    length : float, optional
        Half length of the accepted mirror domain, a benchmark default [m].

    Returns
    -------
    ParticleField
        Field with Bz = B0(1+z^2/length^2) and Br = -r B0 z/length^2 [-].

    Convention
    ----------
    The domain is abs(z) < length and r < length/4. On axis, s=z and the
    nominal mirror ratio at the two ends is two. This is a current-carrying
    analytic benchmark, not an assertion of a vacuum-coil solution.

    Applicability
    -------------
    Machine-independent. Magnetized mirror and harmonic-bounce benchmark.

    Provenance
    ----------
    .. [1] Curl of A=(-y f(z)/2, x f(z)/2, 0), f=B0(1+z^2/L^2),
       an exactly divergence-free benchmark for issue #2005.
    """
    _positive(B0=B0, length=length)
    c = B0 / length**2

    def B(x):
        a, b, z = x
        return np.array([-c * a * z, -c * b * z, B0 + c * z * z])

    def J(x):
        a, b, z = x
        return np.array(
            [[-c * z, 0.0, -c * a], [0.0, -c * z, -c * b], [0.0, 0.0, 2 * c * z]]
        )

    return ParticleField(
        B,
        jacobian=J,
        vector_potential=lambda x: (
            np.array([-x[1], x[0], 0.0]) * (B0 + c * x[2] ** 2) / 2
        ),
        contains=lambda x: abs(x[2]) < length and np.hypot(x[0], x[1]) < length / 4,
        name="solenoidal parabolic magnetic mirror",
        provenance={
            "gauge": "A=(-y f/2,x f/2,0)",
            "length_m": length,
            "B_min_T": B0,
            "B0_T": B0,
            "nominal_end_ratio": 2.0,
            "on_axis_arc_coordinate": "z",
        },
    )

def slab_particle_field(B0=1.0, length=30.0):
    """Construct a divergence-free slab gradient Bz = B0 (1 + x/length).

    Parameters
    ----------
    B0 : float, optional
        Positive reference field, a benchmark default [T].
    length : float, optional
        Positive gradient scale, a benchmark default [m].

    Returns
    -------
    ParticleField
        Slab field restricted to abs(x) < length/2 [-].

    Applicability
    -------------
    Machine-independent. Weak-gradient drift benchmark when rho/length is small.

    Provenance
    ----------
    .. [1] F. F. Chen, Introduction to Plasma Physics and Controlled Fusion,
       3rd ed. (2016), Sec. 2.3.
    """
    _positive(B0=B0, length=length)
    J = np.zeros((3, 3))
    J[2, 0] = B0 / length
    return ParticleField(
        lambda x: np.array([0.0, 0.0, B0 * (1 + x[0] / length)]),
        jacobian=lambda x: J.copy(),
        contains=lambda x: abs(x[0]) < length / 2,
        vector_potential=lambda x: np.array(
            [0.0, B0 * (x[0] + x[0] ** 2 / (2 * length)), 0.0]
        ),
        name="divergence-free slab gradient",
        provenance={"gauge": "Ay=B0(x+x^2/(2L))", "length_m": length, "B0_T": B0},
    )

def uniform_particle_field(B=(0.0, 0.0, 1.0), E=(0.0, 0.0, 0.0)):
    """Construct a uniform magnetic field with a consistent crossed electric field.

    Parameters
    ----------
    B : array_like, optional
        Nonzero Cartesian magnetic field; unit z field is a reference default [T].
    E : array_like, optional
        Uniform electric field; zero is the magnetic-only default [V/m].

    Returns
    -------
    ParticleField
        Uniform field with symmetric gauge A = B cross x / 2 [-].

    Applicability
    -------------
    Machine-independent. Exact uniform-field benchmark.

    Provenance
    ----------
    .. [1] F. F. Chen, Introduction to Plasma Physics and Controlled Fusion,
       3rd ed. (2016), Sec. 2.2.
    """
    B, E = np.asarray(B, float), np.asarray(E, float)
    if (
        B.shape != (3,)
        or E.shape != (3,)
        or not np.all(np.isfinite(np.r_[B, E]))
        or np.linalg.norm(B) == 0
    ):
        raise ValueError("finite 3-vector fields and nonzero B required")
    B, E = B.copy(), E.copy()
    electric = np.any(E != 0)
    return ParticleField(
        lambda x: B.copy(),
        E=(lambda x: E.copy()) if electric else None,
        potential=(lambda x: -float(E @ x)) if electric else None,
        jacobian=lambda x: np.zeros((3, 3)),
        vector_potential=lambda x: np.cross(B, x) / 2,
        derivative_step=1e-4,
        name="uniform Cartesian field",
        provenance={
            "gauge": "A=B cross x/2; Phi=-E dot x", "geometry": "uniform",
            "B_T": B.tolist(), "E_V_per_m": E.tolist(),
        },
    )