Wakefields

Resistive Wall Wakefield Classes

ResistiveWallWakefield dataclass

ResistiveWallWakefield(radius: float, conductivity: float, relaxation_time: float, geometry: Geometry = ROUND)

Bases: ResistiveWallWakefieldBase, ImpedanceWakefield

              flowchart TD
              beamphysics.wakefields.ResistiveWallWakefield[ResistiveWallWakefield]
              beamphysics.wakefields.resistive_wall.base.ResistiveWallWakefieldBase[ResistiveWallWakefieldBase]
              beamphysics.wakefields.impedance.ImpedanceWakefield[ImpedanceWakefield]
              beamphysics.wakefields.base.WakefieldBase[WakefieldBase]

                              beamphysics.wakefields.resistive_wall.base.ResistiveWallWakefieldBase --> beamphysics.wakefields.ResistiveWallWakefield
                                beamphysics.wakefields.base.WakefieldBase --> beamphysics.wakefields.resistive_wall.base.ResistiveWallWakefieldBase
                

                beamphysics.wakefields.impedance.ImpedanceWakefield --> beamphysics.wakefields.ResistiveWallWakefield
                                beamphysics.wakefields.base.WakefieldBase --> beamphysics.wakefields.impedance.ImpedanceWakefield
                



              click beamphysics.wakefields.ResistiveWallWakefield href "" "beamphysics.wakefields.ResistiveWallWakefield"
              click beamphysics.wakefields.resistive_wall.base.ResistiveWallWakefieldBase href "" "beamphysics.wakefields.resistive_wall.base.ResistiveWallWakefieldBase"
              click beamphysics.wakefields.impedance.ImpedanceWakefield href "" "beamphysics.wakefields.impedance.ImpedanceWakefield"
              click beamphysics.wakefields.base.WakefieldBase href "" "beamphysics.wakefields.base.WakefieldBase"
            

Accurate impedance-based resistive wall wakefield model.

This class computes the short-range resistive wall wakefield by numerically evaluating the full longitudinal impedance Z(k) and transforming to the wakefield W(z) via FFT. This approach captures all physical effects including the anomalous skin effect at high frequencies.

Physics Background

When a relativistic charged particle travels through a conducting beam pipe, electromagnetic fields penetrate into the conductor due to the finite conductivity. The skin depth δ decreases with frequency:

\[\delta = \sqrt{\frac{2}{\omega \mu_0 \sigma}}\]

At very high frequencies (relevant for short-range wakes), the electron mean free path becomes comparable to the skin depth, and the simple Ohmic model breaks down. The Drude model accounts for this through a frequency- dependent conductivity with relaxation time τ:

\[\sigma(k) = \frac{\sigma_0}{1 - ikc\tau}\]

The longitudinal impedance per unit length for a round pipe is:

\[Z(k) = \frac{Z_0}{2\pi a} \cdot \frac{\zeta(k)}{1 - i k a \zeta(k) / 2}\]

where ζ(k) is the surface impedance and a is the pipe radius.

The wakefield is obtained via cosine transform (for z ≤ 0):

\[W(z) = \frac{2c}{\pi} \int_0^{\infty} \text{Re}[Z(k)] \cos(kz) \, dk\]

Characteristic Scales

• s₀: Characteristic length scale:

$\(s_0 = \left( \frac{2a^2}{Z_0 \sigma_0} \right)^{1/3}\)$

For copper with a = 2.5 mm, s₀ ≈ 8 µm.

• W₀: Characteristic wake amplitude at z = 0:
• Round: \(W_0 = c Z_0 / (\pi a^2)\)
• Flat: \(W_0 = c Z_0 \pi / (16 a^2)\)

Geometry

• Round (circular pipe): Closed-form impedance expression.
• Flat (parallel plates): Requires numerical integration over transverse modes; more computationally expensive.

Parameters:

• radius (float) –

  Radius of the beam pipe [m]. For flat geometry, this is the half-gap.

• conductivity (float) –

  Electrical DC conductivity of the wall material [S/m].

• relaxation_time (float) –

  Drude-model relaxation time τ of the conductor [s]. Typical values: Cu ≈ 27 fs, Al ≈ 8 fs.

• geometry (str, default: ROUND ) –

  Geometry of the beam pipe: 'round' or 'flat'. Default is 'round'.

Attributes:

• s0 (float) –

  Characteristic length scale [m]

• W0 (float) –

  Characteristic wake amplitude [V/C/m]

Notes

• This model is ~10-20× slower than ResistiveWallPseudomode but more accurate, especially for the first oscillation peak.
• The FFT-based convolution uses O(N log N) complexity for density convolution, compared to O(N) for the pseudomode model.
• For flat geometry, impedance calculation requires numerical integration and is significantly slower than round geometry.

See Also

ResistiveWallPseudomode : Fast approximate model using damped sinusoid

References

.. [1] K. Bane and G. Stupakov, "Resistive wall wakefield in the LCLS undulator beam pipe," SLAC-PUB-10707 (2004). https://www.slac.stanford.edu/cgi-wrap/getdoc/slac-pub-10707.pdf

.. [2] A. Chao, "Physics of Collective Beam Instabilities in High Energy Accelerators," Wiley, 1993, Chapter 2.

Examples:

::

wake = ResistiveWallWakefield.from_material(
    "copper-slac-pub-10707", radius=2.5e-3, geometry="round"
)
wake.wake(-10e-6)  # Wake at 10 µm behind source
wake.impedance(1e5)  # Impedance at k = 100/mm

Methods:

• __call__ –

  Evaluate the wakefield at position z (convenience method).

• convolve_density –

  Compute integrated wakefield from a charge density distribution.

• from_material –

  Create a wakefield from a known material preset.

• impedance –

  Evaluate the impedance at wavenumber k.

• material_from_properties –

  Attempt to identify the material based on conductivity and relaxation time.

• particle_kicks –

  Compute wakefield-induced energy kicks per unit length.

• plot –

  Plot the resistive wall wakefield W(z).

• plot_impedance –

  Plot the resistive wall impedance Z(k).

• wake –

  Evaluate the wakefield at position z.

W0 property

W0

Characteristic wake amplitude W₀ at z=0 [V/C/m].

From SLAC-PUB-10707: - Round: W₀ = c·Z₀ / (π·a²) - Flat: W₀ = c·Z₀·π / (16·a²) = W₀_round · (π²/16)

The dimensionless scaled wake is Ŵ(z/s₀) = W(z) / W₀.

k_max property writable

k_max: float

Maximum wavenumber for integration.

s0 property

s0

Characteristic length scale s₀ of the resistive wall wakefield [m].

From SLAC-PUB-10707 Eq. (5):

\[s_0 = \left( \frac{2 a^2}{Z_0 \sigma_0} \right)^{1/3}\]

where a is the pipe radius (round) or half-gap (flat), Z₀ is the impedance of free space, and σ₀ is the DC conductivity.

__call__

__call__(z: float | ndarray) -> float | ndarray

Evaluate the wakefield at position z (convenience method).

Source code in beamphysics/wakefields/resistive_wall/base.py

def __call__(self, z: float | np.ndarray) -> float | np.ndarray:
    """Evaluate the wakefield at position z (convenience method)."""
    return self.wake(z)

convolve_density

convolve_density(density: ndarray, dz: float, offset: float = 0, include_self_kick: bool = True, plot: bool = False, ax=None) -> ndarray

Compute integrated wakefield from a charge density distribution.

This computes the causal wake potential:

\[V(z_i) = \sum_{j>i} Q_j \cdot W(z_j - z_i) + \frac{1}{2} Q_i \cdot W(0)\]

where only particles ahead (larger z index) contribute to the wake felt by each particle. This is mathematically a correlation, not a convolution, and is computed efficiently using FFT.

Parameters:

• density (ndarray) –

  Charge density array [C/m]. Index 0 is the tail (back), index n-1 is the head (front).

• dz (float) –

  Grid spacing [m].

• offset (float, default: 0 ) –

  Offset for the z coordinate [m]. Default is 0.

• include_self_kick (bool, default: True ) –

  Whether to include the half self-kick term. Default is True.

• plot (bool, default: False ) –

  If True, plot the density profile and wake potential. Default is False.

• ax (Axes, default: None ) –

  Axes to plot on. If None and plot=True, creates a new figure.

Returns:

• integrated_wake ( ndarray ) –

  Integrated longitudinal wakefield [V].

Source code in beamphysics/wakefields/base.py

def convolve_density(
    self,
    density: np.ndarray,
    dz: float,
    offset: float = 0,
    include_self_kick: bool = True,
    plot: bool = False,
    ax=None,
) -> np.ndarray:
    """
    Compute integrated wakefield from a charge density distribution.

    This computes the causal wake potential:

    $$V(z_i) = \\sum_{j>i} Q_j \\cdot W(z_j - z_i) + \\frac{1}{2} Q_i \\cdot W(0)$$

    where only particles ahead (larger z index) contribute to the wake
    felt by each particle. This is mathematically a correlation, not
    a convolution, and is computed efficiently using FFT.

    Parameters
    ----------
    density : np.ndarray
        Charge density array [C/m]. Index 0 is the tail (back),
        index n-1 is the head (front).
    dz : float
        Grid spacing [m].
    offset : float, optional
        Offset for the z coordinate [m]. Default is 0.
    include_self_kick : bool, optional
        Whether to include the half self-kick term. Default is True.
    plot : bool, optional
        If True, plot the density profile and wake potential. Default is False.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None and plot=True, creates a new figure.

    Returns
    -------
    integrated_wake : np.ndarray
        Integrated longitudinal wakefield [V].
    """
    from scipy.signal import correlate

    n = len(density)
    charge = density * dz  # Convert to charge per bin [C]

    # Build causal wake array: W[k] = wake at lag k*dz
    # W[0] = W(0), W[1] = W(-dz), W[2] = W(-2*dz), ...
    z_wake = -np.arange(n) * dz + offset
    W_causal = self.wake(z_wake)

    # Correlation: result[i] = sum_k charge[i+k] * W[k]
    # This gives the wake from particles ahead (larger indices)
    corr_result = correlate(charge, W_causal, mode="full")

    # Extract the portion corresponding to non-negative lags
    # correlate output has length 2*n-1, with zero lag at index n-1
    integrated_wake = corr_result[n - 1 : 2 * n - 1]

    # Correlation includes full W[0] contribution from self
    # Remove it, then add half if self-kick is requested
    W0 = W_causal[0]
    integrated_wake = integrated_wake - charge * W0
    if include_self_kick:
        integrated_wake = integrated_wake + 0.5 * W0 * charge

    if plot:
        self._plot_convolve_density(density, dz, integrated_wake, ax=ax)

    return integrated_wake

from_material classmethod

from_material(material: str, radius: float, geometry: Geometry | str = ROUND)

Create a wakefield from a known material preset.

Parameters:

• material (str) –

  Material name. Must be in list(cls.MATERIALS)

• radius (float) –

  Pipe radius [m]

• geometry (str, default: ROUND ) –

  Geometry type: 'round' or 'flat'

Source code in beamphysics/wakefields/resistive_wall/base.py

@classmethod
def from_material(
    cls,
    material: str,
    radius: float,
    geometry: Geometry | str = Geometry.ROUND,
):
    """
    Create a wakefield from a known material preset.

    Parameters
    ----------
    material : str
        Material name. Must be in list(cls.MATERIALS)
    radius : float
        Pipe radius [m]
    geometry : str
        Geometry type: 'round' or 'flat'
    """
    if material not in cls.MATERIALS:
        raise ValueError(
            f"Unknown material {material!r}. Available: {list(cls.MATERIALS)}"
        )

    props = cls.MATERIALS[material]
    return cls(
        radius=radius,
        conductivity=props["conductivity"],
        relaxation_time=props["relaxation_time"],
        geometry=geometry,
    )

impedance

impedance(k)

Evaluate the impedance at wavenumber k.

Parameters:

• k (float or ndarray) –

  Wavenumber [1/m].

Returns:

• Z ( complex or ndarray ) –

  Impedance [Ohm/m].

Source code in beamphysics/wakefields/resistive_wall/impedance.py

def impedance(self, k):
    """
    Evaluate the impedance at wavenumber k.

    Parameters
    ----------
    k : float or np.ndarray
        Wavenumber [1/m].

    Returns
    -------
    Z : complex or np.ndarray
        Impedance [Ohm/m].
    """
    if self.geometry == Geometry.ROUND:
        return longitudinal_impedance_round(
            k, self.radius, self.conductivity, self._ctau
        )
    else:  # FLAT
        return longitudinal_impedance_flat(
            k, self.radius, self.conductivity, self._ctau
        )

material_from_properties

material_from_properties(tol=0.01)

Attempt to identify the material based on conductivity and relaxation time.

Parameters:

• tol (float, default: 0.01 ) –

  Relative tolerance for matching, default is 1% (0.01)

Returns:

• material ( str or None ) –

  Name of the matched material, or None if no match is found.

Source code in beamphysics/wakefields/resistive_wall/base.py

def material_from_properties(self, tol=0.01):
    """
    Attempt to identify the material based on conductivity and relaxation time.

    Parameters
    ----------
    tol : float, optional
        Relative tolerance for matching, default is 1% (0.01)

    Returns
    -------
    material : str or None
        Name of the matched material, or None if no match is found.
    """
    from math import isclose

    for name, props in self.MATERIALS.items():
        cond_match = isclose(self.conductivity, props["conductivity"], rel_tol=tol)
        tau_match = isclose(
            self.relaxation_time, props["relaxation_time"], rel_tol=tol
        )
        if cond_match and tau_match:
            return name
    return None

particle_kicks

particle_kicks(z: ndarray, weight: ndarray, include_self_kick: bool = True, n_bins: int = None) -> ndarray

Compute wakefield-induced energy kicks per unit length.

Uses FFT-based density convolution with interpolation back to particle positions. This is O(N + M log M) where N is the number of particles and M is the number of grid points.

Parameters:

• z (ndarray) –

  Particle positions [m].

• weight (ndarray) –

  Particle charges [C].

• include_self_kick (bool, default: True ) –

  Whether to include the self-kick term. Default is True.

• n_bins (int, default: None ) –

  Number of bins for the density grid. Default is max(100, N//10).

Returns:

• kicks ( ndarray ) –

  Energy kick per unit length at each particle [eV/m].

Source code in beamphysics/wakefields/impedance.py

def particle_kicks(
    self,
    z: np.ndarray,
    weight: np.ndarray,
    include_self_kick: bool = True,
    n_bins: int = None,
) -> np.ndarray:
    """
    Compute wakefield-induced energy kicks per unit length.

    Uses FFT-based density convolution with interpolation back to
    particle positions. This is O(N + M log M) where N is the number
    of particles and M is the number of grid points.

    Parameters
    ----------
    z : np.ndarray
        Particle positions [m].
    weight : np.ndarray
        Particle charges [C].
    include_self_kick : bool, optional
        Whether to include the self-kick term. Default is True.
    n_bins : int, optional
        Number of bins for the density grid. Default is max(100, N//10).

    Returns
    -------
    kicks : np.ndarray
        Energy kick per unit length at each particle [eV/m].
    """
    from scipy.interpolate import interp1d

    z = np.asarray(z)
    weight = np.asarray(weight)
    n = len(z)

    if n_bins is None:
        n_bins = max(100, n // 10)

    # Bin particles into density distribution
    z_min, z_max = z.min(), z.max()
    z_range = z_max - z_min
    if z_range == 0:
        z_range = 1e-6  # Avoid division by zero for single particle

    # Add padding to avoid edge effects
    padding = 0.1 * z_range
    z_min_padded = z_min - padding
    z_max_padded = z_max + padding

    dz = (z_max_padded - z_min_padded) / n_bins
    z_grid = np.linspace(z_min_padded + dz / 2, z_max_padded - dz / 2, n_bins)

    # Create density histogram (charge per bin)
    density, bin_edges = np.histogram(
        z, bins=n_bins, range=(z_min_padded, z_max_padded), weights=weight
    )
    # Convert to charge density [C/m]
    density = density / dz

    # Compute wake potential via FFT convolution
    wake_potential = self.convolve_density(
        density, dz, include_self_kick=include_self_kick
    )

    # Interpolate wake potential to particle positions
    interp = interp1d(
        z_grid, wake_potential, kind="linear", bounds_error=False, fill_value=0.0
    )
    kicks = -interp(z)  # Negative sign: energy loss for trailing particles

    return kicks

plot

plot(zmax=None, zmin=0, n=200, normalized=False, ax=None)

Plot the resistive wall wakefield W(z).

Parameters:

• zmax (float, default: None ) –

  Maximum trailing distance [m]. Defaults to 100 * s0.

• zmin (float, default: 0 ) –

  Minimum trailing distance [m]. Default is 0.

• n (int, default: 200 ) –

  Number of points. Default is 200.

• normalized (bool, default: False ) –

  If True, plot dimensionless Ŵ(z/s₀) = W(z)/W₀ vs z/s₀. Default is False.

• ax (Axes, default: None ) –

  Axes to plot on. If None, creates a new figure.

Source code in beamphysics/wakefields/resistive_wall/impedance.py

def plot(self, zmax=None, zmin=0, n=200, normalized=False, ax=None):
    """
    Plot the resistive wall wakefield W(z).

    Parameters
    ----------
    zmax : float, optional
        Maximum trailing distance [m]. Defaults to 100 * s0.
    zmin : float, optional
        Minimum trailing distance [m]. Default is 0.
    n : int, optional
        Number of points. Default is 200.
    normalized : bool, optional
        If True, plot dimensionless Ŵ(z/s₀) = W(z)/W₀ vs z/s₀.
        Default is False.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if zmax is None:
        zmax = 20 * self.s0

    zlist = np.linspace(zmin, zmax, n)
    Wz = self.wake(-zlist)

    if ax is None:
        _, ax = plt.subplots()

    if normalized:
        ax.plot(zlist / self.s0, Wz / self.W0)
        ax.set_xlabel(r"$|z|/s_0$")
        ax.set_ylabel(r"$W(z)/W_0$")
    else:
        ax.plot(zlist * 1e6, Wz * 1e-12)
        ax.set_xlabel(r"Distance behind source $|z|$ (µm)")
        ax.set_ylabel(r"$W(z)$ (V/pC/m)")

plot_impedance

plot_impedance(kmax: float = None, kmin: float = 0, n: int = 500, ax=None)

Plot the resistive wall impedance Z(k).

Parameters:

• kmax (float, default: None ) –

  Maximum wavenumber [1/m]. Defaults to 10/s0.

• kmin (float, default: 0 ) –

  Minimum wavenumber [1/m]. Default is 0.

• n (int, default: 500 ) –

  Number of points. Default is 500.

• ax (Axes, default: None ) –

  Axes to plot on. If None, creates a new figure.

Source code in beamphysics/wakefields/resistive_wall/impedance.py

def plot_impedance(
    self, kmax: float = None, kmin: float = 0, n: int = 500, ax=None
):
    """
    Plot the resistive wall impedance Z(k).

    Parameters
    ----------
    kmax : float, optional
        Maximum wavenumber [1/m]. Defaults to 10/s0.
    kmin : float, optional
        Minimum wavenumber [1/m]. Default is 0.
    n : int, optional
        Number of points. Default is 500.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if kmax is None:
        kmax = 10 / self.s0
    super().plot_impedance(kmax=kmax, kmin=kmin, n=n, ax=ax)

wake

wake(z: ndarray | float, n_fft: int = 4096) -> ndarray | float

Evaluate the wakefield at position z.

If a wakefield function was provided, uses it directly. Otherwise, computes via cosine transform of Re[Z(k)]: - Scalars: numerical quadrature (accurate) - Arrays: FFT with interpolation (fast)

Parameters:

• z (float or ndarray) –

  Longitudinal position [m].

• n_fft (int, default: 4096 ) –

  Number of FFT points for array evaluation. Default is 4096.

Returns:

• W ( float or ndarray ) –

  Wakefield value [V/C/m]. Returns 0 for z > 0.

Source code in beamphysics/wakefields/impedance.py

def wake(self, z: np.ndarray | float, n_fft: int = 4096) -> np.ndarray | float:
    """
    Evaluate the wakefield at position z.

    If a wakefield function was provided, uses it directly.
    Otherwise, computes via cosine transform of Re[Z(k)]:
    - Scalars: numerical quadrature (accurate)
    - Arrays: FFT with interpolation (fast)

    Parameters
    ----------
    z : float or np.ndarray
        Longitudinal position [m].
    n_fft : int, optional
        Number of FFT points for array evaluation. Default is 4096.

    Returns
    -------
    W : float or np.ndarray
        Wakefield value [V/C/m]. Returns 0 for z > 0.
    """
    if self._wakefield_func is not None:
        return self._wakefield_func(z)

    z = np.asarray(z)
    scalar_input = z.ndim == 0
    z = np.atleast_1d(z)

    if scalar_input:
        # Single point: use quadrature for accuracy
        from scipy.integrate import quad

        z_val = float(z[0])
        if z_val > 0:
            return 0.0

        def integrand(k):
            return np.real(self._impedance_func(k)) * np.cos(k * z_val)

        result, _ = quad(integrand, 0, self._k_max, limit=200)
        return (2 * c_light / np.pi) * result
    else:
        # Array: use FFT for speed
        from scipy.fft import irfft
        from scipy.interpolate import interp1d

        k_max = self._k_max
        dk = k_max / (n_fft - 1)
        k_grid = np.linspace(0, k_max, n_fft)

        Zk_grid = self._impedance_func(k_grid)
        ReZ = np.real(Zk_grid)

        n_full = 2 * (n_fft - 1)
        dz = 2 * np.pi / (n_full * dk)
        z_grid = np.arange(n_full) * dz

        W_grid = irfft(ReZ, n=n_full) * n_full * dk * (c_light / np.pi)

        interp = interp1d(
            z_grid, W_grid, kind="cubic", bounds_error=False, fill_value=0.0
        )

        result = interp(-z)
        result = np.where(z > 0, 0.0, result)

        return result

ResistiveWallPseudomode dataclass

ResistiveWallPseudomode(radius: float, conductivity: float, relaxation_time: float, geometry: Geometry = ROUND)

Bases: ResistiveWallWakefieldBase, PseudomodeWakefield

              flowchart TD
              beamphysics.wakefields.ResistiveWallPseudomode[ResistiveWallPseudomode]
              beamphysics.wakefields.resistive_wall.base.ResistiveWallWakefieldBase[ResistiveWallWakefieldBase]
              beamphysics.wakefields.pseudomode.PseudomodeWakefield[PseudomodeWakefield]
              beamphysics.wakefields.base.WakefieldBase[WakefieldBase]

                              beamphysics.wakefields.resistive_wall.base.ResistiveWallWakefieldBase --> beamphysics.wakefields.ResistiveWallPseudomode
                                beamphysics.wakefields.base.WakefieldBase --> beamphysics.wakefields.resistive_wall.base.ResistiveWallWakefieldBase
                

                beamphysics.wakefields.pseudomode.PseudomodeWakefield --> beamphysics.wakefields.ResistiveWallPseudomode
                                beamphysics.wakefields.base.WakefieldBase --> beamphysics.wakefields.pseudomode.PseudomodeWakefield
                



              click beamphysics.wakefields.ResistiveWallPseudomode href "" "beamphysics.wakefields.ResistiveWallPseudomode"
              click beamphysics.wakefields.resistive_wall.base.ResistiveWallWakefieldBase href "" "beamphysics.wakefields.resistive_wall.base.ResistiveWallWakefieldBase"
              click beamphysics.wakefields.pseudomode.PseudomodeWakefield href "" "beamphysics.wakefields.pseudomode.PseudomodeWakefield"
              click beamphysics.wakefields.base.WakefieldBase href "" "beamphysics.wakefields.base.WakefieldBase"
            

Fast pseudomode-based resistive wall wakefield model.

This class models the short-range resistive wall wakefield using a damped sinusoidal "pseudomode" approximation. When a relativistic charged particle travels through a conducting beam pipe, it induces image currents in the pipe walls. Due to the finite conductivity of the wall material, these currents penetrate into the conductor and dissipate energy, creating a wakefield that acts on trailing particles.

Physics Background

The resistive wall impedance arises from the skin effect in the conducting walls. At high frequencies (short distances), the AC conductivity of metals deviates from DC behavior due to the Drude relaxation time τ:

\[\sigma(\omega) = \frac{\sigma_0}{1 - i\omega\tau}\]

This frequency dependence causes the wakefield to oscillate and decay exponentially behind the source particle. The wakefield can be well approximated by a single damped sinusoid (pseudomode):

\[W(z) = A \cdot e^{d \cdot z} \cdot \sin(k_r z + \phi) \quad (z \le 0)\]

where the parameters \(k_r\) (oscillation wavenumber) and \(Q_r\) (quality factor, related to decay rate \(d = k_r / 2Q_r\)) depend on the dimensionless relaxation parameter \(\Gamma = c\tau / s_0\).

Characteristic Scales

• s₀: Characteristic length scale where the wake transitions from the \(1/\sqrt{z}\) DC behavior to oscillatory AC behavior:

$\(s_0 = \left( \frac{2a^2}{Z_0 \sigma_0} \right)^{1/3}\)$

• Γ: Dimensionless relaxation parameter \(\Gamma = c\tau / s_0\). For copper, Γ ≈ 0.8; for aluminum, Γ ≈ 0.2.

• W₀: Characteristic wake amplitude at z = 0:

• Round: \(W_0 = c Z_0 / (\pi a^2)\)
• Flat: \(W_0 = c Z_0 \pi / (16 a^2)\)

Algorithm

The pseudomode form enables an O(N) algorithm for computing particle kicks, compared to O(N²) or O(N log N) for general wakefields. This makes it suitable for multi-pass tracking simulations.

Parameters:

• radius (float) –

  Radius of the beam pipe [m]. For flat geometry, this is half the gap.

• conductivity (float) –

  Electrical DC conductivity of the wall material [S/m].

• relaxation_time (float) –

  Drude-model relaxation time τ of the conductor [s]. Typical values: Cu ≈ 27 fs, Al ≈ 8 fs.

• geometry (str, default: ROUND ) –

  Geometry of the beam pipe: 'round' or 'flat'. Default is 'round'.

Attributes:

• s0 (float) –

  Characteristic length scale [m]

• Gamma (float) –

  Dimensionless relaxation parameter Γ = cτ/s₀

• kr (float) –

  Resonant wavenumber [1/m]

• Qr (float) –

  Quality factor of the pseudomode

Notes

• The pseudomode approximation has ~10-20% error compared to the full impedance model, primarily in the first oscillation peak.
• The polynomial fits for k_r and Q_r are valid for Γ ≲ 3. A warning is issued for larger values.
• Materials with known conductivity and τ values are available via from_material().

See Also

ResistiveWallWakefield : Accurate impedance-based model (slower)

References

.. [1] K. Bane and G. Stupakov, "Resistive wall wakefield in the LCLS undulator beam pipe," SLAC-PUB-10707 (2004). https://www.slac.stanford.edu/cgi-wrap/getdoc/slac-pub-10707.pdf

.. [2] D. Sagan, "Bmad Manual," Section 24.6: Short-Range Wakefields. https://www.classe.cornell.edu/bmad/manual.html

.. [3] A. Chao, "Physics of Collective Beam Instabilities in High Energy Accelerators," Wiley, 1993, Chapter 2.

Examples:

::

wake = ResistiveWallPseudomode.from_material(
    "copper-slac-pub-10707", radius=2.5e-3, geometry="round"
)
wake.wake(-10e-6)  # Wake at 10 µm behind source

Methods:

• __call__ –

  Evaluate the wakefield at position z (convenience method).

• convolve_density –

  Compute integrated wakefield from a charge density distribution.

• from_material –

  Create a wakefield from a known material preset.

• impedance –

  Compute the impedance Z(k) analytically.

• material_from_properties –

  Attempt to identify the material based on conductivity and relaxation time.

• particle_kicks –

  Compute short-range wakefield energy kicks per unit length.

• plot –

  Plot the resistive wall wakefield W(z).

• plot_impedance –

  Plot the resistive wall impedance Z(k).

• to_bmad –

  Export wakefield in Bmad format.

• wake –

  Evaluate the pseudomode wakefield at position z.

Gamma property

Gamma

Dimensionless relaxation time Γ = c * τ / s₀.

Qr property

Qr

Dimensionless quality factor Q_r of the wakefield pseudomode.

W0 property

W0

Characteristic wake amplitude W₀ at z=0 [V/C/m].

From SLAC-PUB-10707: - Round: W₀ = c·Z₀ / (π·a²) - Flat: W₀ = c·Z₀·π / (16·a²) = W₀_round · (π²/16)

The dimensionless scaled wake is Ŵ(z/s₀) = W(z) / W₀.

kr property

kr

Real-valued wave number k_r of the wakefield pseudomode [1/m].

modes property

modes: list

List of Pseudomode objects.

n_modes property

n_modes: int

Number of pseudomodes.

s0 property

s0

Characteristic length scale s₀ of the resistive wall wakefield [m].

From SLAC-PUB-10707 Eq. (5):

\[s_0 = \left( \frac{2 a^2}{Z_0 \sigma_0} \right)^{1/3}\]

where a is the pipe radius (round) or half-gap (flat), Z₀ is the impedance of free space, and σ₀ is the DC conductivity.

__call__

__call__(z: float | ndarray) -> float | ndarray

Evaluate the wakefield at position z (convenience method).

Source code in beamphysics/wakefields/resistive_wall/base.py

def __call__(self, z: float | np.ndarray) -> float | np.ndarray:
    """Evaluate the wakefield at position z (convenience method)."""
    return self.wake(z)

convolve_density

convolve_density(density: ndarray, dz: float, offset: float = 0, include_self_kick: bool = True, plot: bool = False, ax=None) -> ndarray

Compute integrated wakefield from a charge density distribution.

This computes the causal wake potential:

\[V(z_i) = \sum_{j>i} Q_j \cdot W(z_j - z_i) + \frac{1}{2} Q_i \cdot W(0)\]

where only particles ahead (larger z index) contribute to the wake felt by each particle. This is mathematically a correlation, not a convolution, and is computed efficiently using FFT.

Parameters:

• density (ndarray) –

  Charge density array [C/m]. Index 0 is the tail (back), index n-1 is the head (front).

• dz (float) –

  Grid spacing [m].

• offset (float, default: 0 ) –

  Offset for the z coordinate [m]. Default is 0.

• include_self_kick (bool, default: True ) –

  Whether to include the half self-kick term. Default is True.

• plot (bool, default: False ) –

  If True, plot the density profile and wake potential. Default is False.

• ax (Axes, default: None ) –

  Axes to plot on. If None and plot=True, creates a new figure.

Returns:

• integrated_wake ( ndarray ) –

  Integrated longitudinal wakefield [V].

Source code in beamphysics/wakefields/base.py

def convolve_density(
    self,
    density: np.ndarray,
    dz: float,
    offset: float = 0,
    include_self_kick: bool = True,
    plot: bool = False,
    ax=None,
) -> np.ndarray:
    """
    Compute integrated wakefield from a charge density distribution.

    This computes the causal wake potential:

    $$V(z_i) = \\sum_{j>i} Q_j \\cdot W(z_j - z_i) + \\frac{1}{2} Q_i \\cdot W(0)$$

    where only particles ahead (larger z index) contribute to the wake
    felt by each particle. This is mathematically a correlation, not
    a convolution, and is computed efficiently using FFT.

    Parameters
    ----------
    density : np.ndarray
        Charge density array [C/m]. Index 0 is the tail (back),
        index n-1 is the head (front).
    dz : float
        Grid spacing [m].
    offset : float, optional
        Offset for the z coordinate [m]. Default is 0.
    include_self_kick : bool, optional
        Whether to include the half self-kick term. Default is True.
    plot : bool, optional
        If True, plot the density profile and wake potential. Default is False.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None and plot=True, creates a new figure.

    Returns
    -------
    integrated_wake : np.ndarray
        Integrated longitudinal wakefield [V].
    """
    from scipy.signal import correlate

    n = len(density)
    charge = density * dz  # Convert to charge per bin [C]

    # Build causal wake array: W[k] = wake at lag k*dz
    # W[0] = W(0), W[1] = W(-dz), W[2] = W(-2*dz), ...
    z_wake = -np.arange(n) * dz + offset
    W_causal = self.wake(z_wake)

    # Correlation: result[i] = sum_k charge[i+k] * W[k]
    # This gives the wake from particles ahead (larger indices)
    corr_result = correlate(charge, W_causal, mode="full")

    # Extract the portion corresponding to non-negative lags
    # correlate output has length 2*n-1, with zero lag at index n-1
    integrated_wake = corr_result[n - 1 : 2 * n - 1]

    # Correlation includes full W[0] contribution from self
    # Remove it, then add half if self-kick is requested
    W0 = W_causal[0]
    integrated_wake = integrated_wake - charge * W0
    if include_self_kick:
        integrated_wake = integrated_wake + 0.5 * W0 * charge

    if plot:
        self._plot_convolve_density(density, dz, integrated_wake, ax=ax)

    return integrated_wake

from_material classmethod

from_material(material: str, radius: float, geometry: Geometry | str = ROUND)

Create a wakefield from a known material preset.

Parameters:

• material (str) –

  Material name. Must be in list(cls.MATERIALS)

• radius (float) –

  Pipe radius [m]

• geometry (str, default: ROUND ) –

  Geometry type: 'round' or 'flat'

Source code in beamphysics/wakefields/resistive_wall/base.py

@classmethod
def from_material(
    cls,
    material: str,
    radius: float,
    geometry: Geometry | str = Geometry.ROUND,
):
    """
    Create a wakefield from a known material preset.

    Parameters
    ----------
    material : str
        Material name. Must be in list(cls.MATERIALS)
    radius : float
        Pipe radius [m]
    geometry : str
        Geometry type: 'round' or 'flat'
    """
    if material not in cls.MATERIALS:
        raise ValueError(
            f"Unknown material {material!r}. Available: {list(cls.MATERIALS)}"
        )

    props = cls.MATERIALS[material]
    return cls(
        radius=radius,
        conductivity=props["conductivity"],
        relaxation_time=props["relaxation_time"],
        geometry=geometry,
    )

impedance

impedance(k: ndarray | float) -> ndarray | complex

Compute the impedance Z(k) analytically.

Sums the contribution from each pseudomode using the analytic Fourier transform.

Parameters:

• k (float or ndarray) –

  Wavenumber [1/m].

Returns:

• Z ( complex or ndarray ) –

  Impedance [Ohm/m].

Source code in beamphysics/wakefields/pseudomode.py

def impedance(self, k: np.ndarray | float) -> np.ndarray | complex:
    """
    Compute the impedance Z(k) analytically.

    Sums the contribution from each pseudomode using the analytic
    Fourier transform.

    Parameters
    ----------
    k : float or np.ndarray
        Wavenumber [1/m].

    Returns
    -------
    Z : complex or np.ndarray
        Impedance [Ohm/m].
    """
    k = np.asarray(k)
    scalar_input = k.ndim == 0
    k = np.atleast_1d(k)

    Z = np.zeros(len(k), dtype=complex)
    for mode in self._modes:
        Z += mode.impedance(k)

    if scalar_input:
        return complex(Z[0])
    return Z

material_from_properties

material_from_properties(tol=0.01)

Attempt to identify the material based on conductivity and relaxation time.

Parameters:

• tol (float, default: 0.01 ) –

  Relative tolerance for matching, default is 1% (0.01)

Returns:

• material ( str or None ) –

  Name of the matched material, or None if no match is found.

Source code in beamphysics/wakefields/resistive_wall/base.py

def material_from_properties(self, tol=0.01):
    """
    Attempt to identify the material based on conductivity and relaxation time.

    Parameters
    ----------
    tol : float, optional
        Relative tolerance for matching, default is 1% (0.01)

    Returns
    -------
    material : str or None
        Name of the matched material, or None if no match is found.
    """
    from math import isclose

    for name, props in self.MATERIALS.items():
        cond_match = isclose(self.conductivity, props["conductivity"], rel_tol=tol)
        tau_match = isclose(
            self.relaxation_time, props["relaxation_time"], rel_tol=tol
        )
        if cond_match and tau_match:
            return name
    return None

particle_kicks

particle_kicks(z: ndarray, weight: ndarray, include_self_kick: bool = True) -> ndarray

Compute short-range wakefield energy kicks per unit length.

Uses an O(N) single-pass algorithm for each mode by exploiting the exponential form of the pseudomode wakefield.

Parameters:

• z (ndarray) –

  Particle positions [m].

• weight (ndarray) –

  Particle charges [C].

• include_self_kick (bool, default: True ) –

  If True, applies the self-kick. Default is True.

Returns:

• kicks ( ndarray ) –

  Wake-induced energy kick per unit length [eV/m].

Source code in beamphysics/wakefields/pseudomode.py

def particle_kicks(
    self,
    z: np.ndarray,
    weight: np.ndarray,
    include_self_kick: bool = True,
) -> np.ndarray:
    """
    Compute short-range wakefield energy kicks per unit length.

    Uses an O(N) single-pass algorithm for each mode by exploiting
    the exponential form of the pseudomode wakefield.

    Parameters
    ----------
    z : np.ndarray
        Particle positions [m].
    weight : np.ndarray
        Particle charges [C].
    include_self_kick : bool, optional
        If True, applies the self-kick. Default is True.

    Returns
    -------
    kicks : np.ndarray
        Wake-induced energy kick per unit length [eV/m].
    """
    z = np.asarray(z)
    weight = np.asarray(weight)

    if z.shape != weight.shape:
        raise ValueError(
            f"Mismatched shapes: z.shape={z.shape}, weight.shape={weight.shape}"
        )
    if z.ndim != 1:
        raise ValueError("z and weight must be 1D arrays")

    # Sort particles from tail to head
    ix = z.argsort()
    z_sorted = z[ix].copy()
    z_sorted -= z_sorted.max()  # Offset for numerical stability
    weight_sorted = weight[ix]

    n = len(z)
    delta_E = np.zeros(n)

    # Process each mode with O(N) algorithm
    for mode in self._modes:
        s = mode.d + 1j * mode.k
        c = mode.A * np.exp(1j * mode.phi)

        # O(N) accumulator for this mode
        b = 0.0 + 0.0j

        # Iterate from tail to head
        for i in range(n - 1, -1, -1):
            zi = z_sorted[i]
            qi = weight_sorted[i]

            # Kick from trailing particles
            delta_E[i] -= np.imag(c * np.exp(s * zi) * b)

            # Add this particle to accumulator
            b += qi * np.exp(-s * zi)

        if include_self_kick:
            delta_E -= 0.5 * mode.A * weight_sorted * np.sin(mode.phi)

    # Restore original order
    kicks = np.empty_like(delta_E)
    kicks[ix] = delta_E

    return kicks

plot

plot(zmax=None, zmin=0, n=200, normalized=False, ax=None)

Plot the resistive wall wakefield W(z).

Parameters:

• zmax (float, default: None ) –

  Maximum trailing distance [m]. Defaults to 10 decay lengths.

• zmin (float, default: 0 ) –

  Minimum trailing distance [m]. Default is 0.

• n (int, default: 200 ) –

  Number of points. Default is 200.

• normalized (bool, default: False ) –

  If True, plot dimensionless Ŵ(z/s₀) = W(z)/W₀ vs z/s₀. Default is False.

• ax (Axes, default: None ) –

  Axes to plot on. If None, creates a new figure.

Source code in beamphysics/wakefields/resistive_wall/pseudomode.py

def plot(self, zmax=None, zmin=0, n=200, normalized=False, ax=None):
    """
    Plot the resistive wall wakefield W(z).

    Parameters
    ----------
    zmax : float, optional
        Maximum trailing distance [m]. Defaults to 10 decay lengths.
    zmin : float, optional
        Minimum trailing distance [m]. Default is 0.
    n : int, optional
        Number of points. Default is 200.
    normalized : bool, optional
        If True, plot dimensionless Ŵ(z/s₀) = W(z)/W₀ vs z/s₀.
        Default is False.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if zmax is None:
        zmax = 1 / (self.kr / (2 * self.Qr)) * 10

    zlist = np.linspace(zmin, zmax, n)
    Wz = self.wake(-zlist)

    if ax is None:
        _, ax = plt.subplots()

    if normalized:
        ax.plot(zlist / self.s0, Wz / self.W0)
        ax.set_xlabel(r"$|z|/s_0$")
        ax.set_ylabel(r"$W(z)/W_0$")
    else:
        ax.plot(zlist * 1e6, Wz * 1e-12)
        ax.set_xlabel(r"Distance behind source $|z|$ (µm)")
        ax.set_ylabel(r"$W(z)$ (V/pC/m)")

plot_impedance

plot_impedance(kmax: float = None, kmin: float = 0, n: int = 500, ax=None)

Plot the resistive wall impedance Z(k).

Parameters:

• kmax (float, default: None ) –

  Maximum wavenumber [1/m]. Defaults to 10/s0.

• kmin (float, default: 0 ) –

  Minimum wavenumber [1/m]. Default is 0.

• n (int, default: 500 ) –

  Number of points. Default is 500.

• ax (Axes, default: None ) –

  Axes to plot on. If None, creates a new figure.

Source code in beamphysics/wakefields/resistive_wall/pseudomode.py

def plot_impedance(
    self, kmax: float = None, kmin: float = 0, n: int = 500, ax=None
):
    """
    Plot the resistive wall impedance Z(k).

    Parameters
    ----------
    kmax : float, optional
        Maximum wavenumber [1/m]. Defaults to 10/s0.
    kmin : float, optional
        Minimum wavenumber [1/m]. Default is 0.
    n : int, optional
        Number of points. Default is 500.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if kmax is None:
        kmax = 10 / self.s0
    super().plot_impedance(kmax=kmax, kmin=kmin, n=n, ax=ax)

to_bmad

to_bmad(file=None, z_max=100, amp_scale=1, scale_with_length=True, z_scale=1)

Export wakefield in Bmad format.

Parameters:

• file (str, default: None ) –

  Output file path. If None, returns string only.

• z_max (float, default: 100 ) –

  Trailing z distance for Bmad.

• amp_scale (float, default: 1 ) –

  Amplitude scaling factor.

• scale_with_length (bool, default: True ) –

  Whether to scale with length.

• z_scale (float, default: 1 ) –

  Z scaling factor.

Returns:

• str –

  Bmad-formatted wakefield string.

Source code in beamphysics/wakefields/resistive_wall/pseudomode.py

    def to_bmad(
        self, file=None, z_max=100, amp_scale=1, scale_with_length=True, z_scale=1
    ):
        """
        Export wakefield in Bmad format.

        Parameters
        ----------
        file : str, optional
            Output file path. If None, returns string only.
        z_max : float
            Trailing z distance for Bmad.
        amp_scale : float
            Amplitude scaling factor.
        scale_with_length : bool
            Whether to scale with length.
        z_scale : float
            Z scaling factor.

        Returns
        -------
        str
            Bmad-formatted wakefield string.
        """
        s = f"""! AC Resistive wall wakefield
! Adapted from SLAC-PUB-10707
!    Material        : {self.material_from_properties()}
!    Conductivity    : {self.conductivity} S/m
!    Relaxation time : {self.relaxation_time} s
!    Geometry        : {self.geometry.value}
"""

        if self.geometry == Geometry.ROUND:
            s += f"!    Radius          : {self.radius} m\n"
        elif self.geometry == Geometry.FLAT:
            s += f"!    full gap        : {2*self.radius} m\n"

        s += f"!    s₀              : {self.s0}  m\n"
        s += f"!    Γ               : {self.Gamma} \n"
        s += "! sr_wake =  \n"

        s += f"{{{z_scale=}, {amp_scale=}, {scale_with_length=}, {z_max=},\n"
        s += PseudomodeWakefield.to_bmad(self) + "}\n"

        if file is not None:
            with open(file, "w") as f:
                f.write(s)

        return s

wake

wake(z: ndarray | float) -> ndarray | float

Evaluate the pseudomode wakefield at position z.

Parameters:

• z (float or ndarray) –

  Longitudinal position [m]. Negative z is behind the source.

Returns:

• W ( float or ndarray ) –

  Wakefield value [V/C/m]. Returns 0 for z > 0.

Source code in beamphysics/wakefields/pseudomode.py

def wake(self, z: np.ndarray | float) -> np.ndarray | float:
    """
    Evaluate the pseudomode wakefield at position z.

    Parameters
    ----------
    z : float or np.ndarray
        Longitudinal position [m]. Negative z is behind the source.

    Returns
    -------
    W : float or np.ndarray
        Wakefield value [V/C/m]. Returns 0 for z > 0.
    """
    z = np.asarray(z)
    scalar_input = z.ndim == 0
    z = np.atleast_1d(z)

    out = np.zeros_like(z, dtype=float)
    mask = z <= 0
    z_masked = z[mask]

    # Sum contributions from all modes
    for mode in self._modes:
        out[mask] += mode(z_masked)

    if scalar_input:
        return float(out[0])
    return out

Base Classes

WakefieldBase

Bases: ABC

              flowchart TD
              beamphysics.wakefields.WakefieldBase[WakefieldBase]

              

              click beamphysics.wakefields.WakefieldBase href "" "beamphysics.wakefields.WakefieldBase"
            

Abstract base class for longitudinal wakefield models.

All wakefield implementations should inherit from this class and implement the required abstract methods.

Sign Convention

The wakefield W(z) is defined for a trailing particle at position z relative to the source particle at z=0:

• z < 0: Behind the source (trailing particle feels the wake)
• z > 0: Ahead of the source (causality requires W=0)

The wake is positive for an energy-losing interaction (the trailing particle loses energy).

Methods:

• wake –

  Evaluate the wakefield at position z [V/C/m]

• impedance –

  Evaluate the impedance at wavenumber k [Ohm/m]

• convolve_density –

  Compute integrated wake from a charge density distribution [V/m]

• particle_kicks –

  Compute momentum kicks for a distribution of particles [eV/m]

• plot –

  Plot the wakefield over a range of z values

• plot_impedance –

  Plot the impedance over a range of wavenumbers

See Also

ParticleGroup.apply_wakefield : Apply wakefield kicks to a ParticleGroup

Subclassed by:

• API Wakefields
  •  ImpedanceWakefield
  •  PseudomodeWakefield
  •  TabularWakefield

Used by:

• API Particles  ParticleGroup
  •  apply_wakefield
  •  wakefield_plot

__call__

__call__(z: ndarray | float) -> ndarray | float

Evaluate the wakefield at position z (convenience method).

This is equivalent to calling wake(z).

Source code in beamphysics/wakefields/base.py

def __call__(self, z: np.ndarray | float) -> np.ndarray | float:
    """
    Evaluate the wakefield at position z (convenience method).

    This is equivalent to calling wake(z).
    """
    return self.wake(z)

convolve_density

convolve_density(density: ndarray, dz: float, offset: float = 0, include_self_kick: bool = True, plot: bool = False, ax=None) -> ndarray

Compute integrated wakefield from a charge density distribution.

This computes the causal wake potential:

\[V(z_i) = \sum_{j>i} Q_j \cdot W(z_j - z_i) + \frac{1}{2} Q_i \cdot W(0)\]

where only particles ahead (larger z index) contribute to the wake felt by each particle. This is mathematically a correlation, not a convolution, and is computed efficiently using FFT.

Parameters:

• density (ndarray) –

  Charge density array [C/m]. Index 0 is the tail (back), index n-1 is the head (front).

• dz (float) –

  Grid spacing [m].

• offset (float, default: 0 ) –

  Offset for the z coordinate [m]. Default is 0.

• include_self_kick (bool, default: True ) –

  Whether to include the half self-kick term. Default is True.

• plot (bool, default: False ) –

  If True, plot the density profile and wake potential. Default is False.

• ax (Axes, default: None ) –

  Axes to plot on. If None and plot=True, creates a new figure.

Returns:

• integrated_wake ( ndarray ) –

  Integrated longitudinal wakefield [V].

Source code in beamphysics/wakefields/base.py

def convolve_density(
    self,
    density: np.ndarray,
    dz: float,
    offset: float = 0,
    include_self_kick: bool = True,
    plot: bool = False,
    ax=None,
) -> np.ndarray:
    """
    Compute integrated wakefield from a charge density distribution.

    This computes the causal wake potential:

    $$V(z_i) = \\sum_{j>i} Q_j \\cdot W(z_j - z_i) + \\frac{1}{2} Q_i \\cdot W(0)$$

    where only particles ahead (larger z index) contribute to the wake
    felt by each particle. This is mathematically a correlation, not
    a convolution, and is computed efficiently using FFT.

    Parameters
    ----------
    density : np.ndarray
        Charge density array [C/m]. Index 0 is the tail (back),
        index n-1 is the head (front).
    dz : float
        Grid spacing [m].
    offset : float, optional
        Offset for the z coordinate [m]. Default is 0.
    include_self_kick : bool, optional
        Whether to include the half self-kick term. Default is True.
    plot : bool, optional
        If True, plot the density profile and wake potential. Default is False.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None and plot=True, creates a new figure.

    Returns
    -------
    integrated_wake : np.ndarray
        Integrated longitudinal wakefield [V].
    """
    from scipy.signal import correlate

    n = len(density)
    charge = density * dz  # Convert to charge per bin [C]

    # Build causal wake array: W[k] = wake at lag k*dz
    # W[0] = W(0), W[1] = W(-dz), W[2] = W(-2*dz), ...
    z_wake = -np.arange(n) * dz + offset
    W_causal = self.wake(z_wake)

    # Correlation: result[i] = sum_k charge[i+k] * W[k]
    # This gives the wake from particles ahead (larger indices)
    corr_result = correlate(charge, W_causal, mode="full")

    # Extract the portion corresponding to non-negative lags
    # correlate output has length 2*n-1, with zero lag at index n-1
    integrated_wake = corr_result[n - 1 : 2 * n - 1]

    # Correlation includes full W[0] contribution from self
    # Remove it, then add half if self-kick is requested
    W0 = W_causal[0]
    integrated_wake = integrated_wake - charge * W0
    if include_self_kick:
        integrated_wake = integrated_wake + 0.5 * W0 * charge

    if plot:
        self._plot_convolve_density(density, dz, integrated_wake, ax=ax)

    return integrated_wake

impedance abstractmethod

impedance(k: ndarray | float) -> ndarray | complex

Evaluate the impedance at wavenumber k.

Parameters:

• k (float or ndarray) –

  Wavenumber [1/m].

Returns:

• Z ( complex or ndarray ) –

  Impedance [Ohm/m].

Source code in beamphysics/wakefields/base.py

@abstractmethod
def impedance(self, k: np.ndarray | float) -> np.ndarray | complex:
    """
    Evaluate the impedance at wavenumber k.

    Parameters
    ----------
    k : float or np.ndarray
        Wavenumber [1/m].

    Returns
    -------
    Z : complex or np.ndarray
        Impedance [Ohm/m].
    """
    pass

particle_kicks

particle_kicks(z: ndarray, weight: ndarray, include_self_kick: bool = True) -> ndarray

Compute wakefield-induced energy kicks per unit length.

This is a default O(N²) implementation. Subclasses may override with more efficient algorithms.

Parameters:

• z (ndarray) –

  Particle positions [m].

• weight (ndarray) –

  Particle charges [C].

• include_self_kick (bool, default: True ) –

  Whether to include the self-kick term. Default is True.

Returns:

• kicks ( ndarray ) –

  Energy kick per unit length at each particle [eV/m].

Source code in beamphysics/wakefields/base.py

def particle_kicks(
    self,
    z: np.ndarray,
    weight: np.ndarray,
    include_self_kick: bool = True,
) -> np.ndarray:
    """
    Compute wakefield-induced energy kicks per unit length.

    This is a default O(N²) implementation. Subclasses may override
    with more efficient algorithms.

    Parameters
    ----------
    z : np.ndarray
        Particle positions [m].
    weight : np.ndarray
        Particle charges [C].
    include_self_kick : bool, optional
        Whether to include the self-kick term. Default is True.

    Returns
    -------
    kicks : np.ndarray
        Energy kick per unit length at each particle [eV/m].
    """
    z = np.asarray(z)
    weight = np.asarray(weight)
    n = len(z)
    kicks = np.zeros(n)

    # O(N²) algorithm: sum contribution from all particles
    for i in range(n):
        for j in range(n):
            if i == j:
                if include_self_kick:
                    # Self-kick uses W(0⁻) / 2
                    kicks[i] -= 0.5 * weight[j] * self._self_kick_value()
            else:
                dz = z[i] - z[j]
                if dz < 0:  # j is ahead of i, so i feels wake from j
                    kicks[i] -= weight[j] * self.wake(dz)

    return kicks

plot

plot(zmax: float = None, zmin: float = 0, n: int = 200, ax=None)

Plot the wakefield over a range of z values.

Parameters:

• zmax (float, default: None ) –

  Maximum trailing distance [m]. If None, uses a sensible default.

• zmin (float, default: 0 ) –

  Minimum trailing distance [m]. Default is 0.

• n (int, default: 200 ) –

  Number of points. Default is 200.

• ax (Axes, default: None ) –

  Axes to plot on. If None, creates a new figure.

Source code in beamphysics/wakefields/base.py

def plot(self, zmax: float = None, zmin: float = 0, n: int = 200, ax=None):
    """
    Plot the wakefield over a range of z values.

    Parameters
    ----------
    zmax : float, optional
        Maximum trailing distance [m]. If None, uses a sensible default.
    zmin : float, optional
        Minimum trailing distance [m]. Default is 0.
    n : int, optional
        Number of points. Default is 200.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if zmax is None:
        zmax = 1e-3  # 1 mm default

    zlist = np.linspace(zmin, zmax, n)
    Wz = self.wake(-zlist)  # Negative z for trailing particles

    if ax is None:
        _, ax = plt.subplots()
    ax.plot(zlist * 1e6, Wz * 1e-12)
    ax.set_xlabel(r"Distance behind source $|z|$ (µm)")
    ax.set_ylabel(r"$W(z)$ (V/pC/m)")

plot_impedance

plot_impedance(kmax: float = None, kmin: float = 0, n: int = 200, ax=None)

Plot the impedance over a range of wavenumbers.

Parameters:

• kmax (float, default: None ) –

  Maximum wavenumber [1/m]. If None, uses a sensible default.

• kmin (float, default: 0 ) –

  Minimum wavenumber [1/m]. Default is 0.

• n (int, default: 200 ) –

  Number of points. Default is 200.

• ax (Axes, default: None ) –

  Axes to plot on. If None, creates a new figure.

Source code in beamphysics/wakefields/base.py

def plot_impedance(
    self, kmax: float = None, kmin: float = 0, n: int = 200, ax=None
):
    """
    Plot the impedance over a range of wavenumbers.

    Parameters
    ----------
    kmax : float, optional
        Maximum wavenumber [1/m]. If None, uses a sensible default.
    kmin : float, optional
        Minimum wavenumber [1/m]. Default is 0.
    n : int, optional
        Number of points. Default is 200.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if kmax is None:
        kmax = 1e6  # 1/µm default

    k = np.linspace(kmin, kmax, n)
    Z = self.impedance(k)

    if ax is None:
        _, ax = plt.subplots()
    ax.plot(k * 1e-3, np.real(Z), label=r"Re[$Z$]")
    ax.plot(k * 1e-3, np.imag(Z), label=r"Im[$Z$]")
    ax.set_xlabel(r"$k$ (1/mm)")
    ax.set_ylabel(r"$Z(k)$ (Ω/m)")
    ax.legend()

wake abstractmethod

wake(z: ndarray | float) -> ndarray | float

Evaluate the wakefield at longitudinal position z.

Parameters:

• z (float or ndarray) –

  Longitudinal position [m]. Negative z is behind the source.

Returns:

• W ( float or ndarray ) –

  Wakefield value [V/C/m]. Returns 0 for z > 0 (causality).

Source code in beamphysics/wakefields/base.py

@abstractmethod
def wake(self, z: np.ndarray | float) -> np.ndarray | float:
    """
    Evaluate the wakefield at longitudinal position z.

    Parameters
    ----------
    z : float or np.ndarray
        Longitudinal position [m]. Negative z is behind the source.

    Returns
    -------
    W : float or np.ndarray
        Wakefield value [V/C/m]. Returns 0 for z > 0 (causality).
    """
    pass

PseudomodeWakefield

PseudomodeWakefield(modes: list)

Bases: WakefieldBase

              flowchart TD
              beamphysics.wakefields.PseudomodeWakefield[PseudomodeWakefield]
              beamphysics.wakefields.base.WakefieldBase[WakefieldBase]

                              beamphysics.wakefields.base.WakefieldBase --> beamphysics.wakefields.PseudomodeWakefield
                


              click beamphysics.wakefields.PseudomodeWakefield href "" "beamphysics.wakefields.PseudomodeWakefield"
              click beamphysics.wakefields.base.WakefieldBase href "" "beamphysics.wakefields.base.WakefieldBase"
            

Wakefield represented as a sum of damped sinusoidal pseudomodes.

Models the longitudinal wakefield as:

\[W(z) = \sum_i A_i \cdot e^{d_i \cdot z} \cdot \sin(k_i \cdot z + \phi_i)\]

This form is used to approximate short-range wakefields such as the resistive wall wake.

Parameters:

• modes (list of Pseudomode or list of dict) –

  List of pseudomodes. Each can be a Pseudomode object or a dict with keys 'A', 'd', 'k', 'phi'.

Notes

• Wakefields are defined for z ≤ 0 (trailing the source particle).
• For z > 0, returns 0 (causality).

Examples:

Single mode::

pm = PseudomodeWakefield([Pseudomode(A=1e15, d=1e4, k=1e5, phi=np.pi/2)])
pm.wake(-10e-6)  # Wake at 10 µm behind source

Multiple modes::

modes = [
    Pseudomode(A=1e15, d=1e4, k=1e5, phi=np.pi/2),
    Pseudomode(A=5e14, d=2e4, k=2e5, phi=np.pi/4),
]
pm = PseudomodeWakefield(modes)

Subclassed by:

• API Wakefields  ResistiveWallPseudomode

Methods:

• __call__ –

  Evaluate the wakefield at position z (convenience method).

• convolve_density –

  Compute integrated wakefield from a charge density distribution.

• impedance –

  Compute the impedance Z(k) analytically.

• particle_kicks –

  Compute short-range wakefield energy kicks per unit length.

• plot –

  Plot the pseudomode wakefield.

• plot_impedance –

  Plot the impedance over a range of wavenumbers.

• to_bmad –

  Format pseudomode parameters as Bmad-compatible strings.

• wake –

  Evaluate the pseudomode wakefield at position z.

Attributes:

• modes (list) –

  List of Pseudomode objects.

• n_modes (int) –

  Number of pseudomodes.

Source code in beamphysics/wakefields/pseudomode.py

def __init__(self, modes: list) -> None:
    self._modes = []
    for mode in modes:
        if isinstance(mode, Pseudomode):
            self._modes.append(mode)
        elif isinstance(mode, dict):
            self._modes.append(Pseudomode(**mode))
        else:
            raise TypeError(f"Mode must be Pseudomode or dict, got {type(mode)}")

modes property

modes: list

List of Pseudomode objects.

n_modes property

n_modes: int

Number of pseudomodes.

__call__

__call__(z: ndarray | float) -> ndarray | float

Evaluate the wakefield at position z (convenience method).

This is equivalent to calling wake(z).

Source code in beamphysics/wakefields/base.py

def __call__(self, z: np.ndarray | float) -> np.ndarray | float:
    """
    Evaluate the wakefield at position z (convenience method).

    This is equivalent to calling wake(z).
    """
    return self.wake(z)

convolve_density

convolve_density(density: ndarray, dz: float, offset: float = 0, include_self_kick: bool = True, plot: bool = False, ax=None) -> ndarray

Compute integrated wakefield from a charge density distribution.

This computes the causal wake potential:

\[V(z_i) = \sum_{j>i} Q_j \cdot W(z_j - z_i) + \frac{1}{2} Q_i \cdot W(0)\]

where only particles ahead (larger z index) contribute to the wake felt by each particle. This is mathematically a correlation, not a convolution, and is computed efficiently using FFT.

Parameters:

• density (ndarray) –

  Charge density array [C/m]. Index 0 is the tail (back), index n-1 is the head (front).

• dz (float) –

  Grid spacing [m].

• offset (float, default: 0 ) –

  Offset for the z coordinate [m]. Default is 0.

• include_self_kick (bool, default: True ) –

  Whether to include the half self-kick term. Default is True.

• plot (bool, default: False ) –

  If True, plot the density profile and wake potential. Default is False.

• ax (Axes, default: None ) –

  Axes to plot on. If None and plot=True, creates a new figure.

Returns:

• integrated_wake ( ndarray ) –

  Integrated longitudinal wakefield [V].

Source code in beamphysics/wakefields/base.py

def convolve_density(
    self,
    density: np.ndarray,
    dz: float,
    offset: float = 0,
    include_self_kick: bool = True,
    plot: bool = False,
    ax=None,
) -> np.ndarray:
    """
    Compute integrated wakefield from a charge density distribution.

    This computes the causal wake potential:

    $$V(z_i) = \\sum_{j>i} Q_j \\cdot W(z_j - z_i) + \\frac{1}{2} Q_i \\cdot W(0)$$

    where only particles ahead (larger z index) contribute to the wake
    felt by each particle. This is mathematically a correlation, not
    a convolution, and is computed efficiently using FFT.

    Parameters
    ----------
    density : np.ndarray
        Charge density array [C/m]. Index 0 is the tail (back),
        index n-1 is the head (front).
    dz : float
        Grid spacing [m].
    offset : float, optional
        Offset for the z coordinate [m]. Default is 0.
    include_self_kick : bool, optional
        Whether to include the half self-kick term. Default is True.
    plot : bool, optional
        If True, plot the density profile and wake potential. Default is False.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None and plot=True, creates a new figure.

    Returns
    -------
    integrated_wake : np.ndarray
        Integrated longitudinal wakefield [V].
    """
    from scipy.signal import correlate

    n = len(density)
    charge = density * dz  # Convert to charge per bin [C]

    # Build causal wake array: W[k] = wake at lag k*dz
    # W[0] = W(0), W[1] = W(-dz), W[2] = W(-2*dz), ...
    z_wake = -np.arange(n) * dz + offset
    W_causal = self.wake(z_wake)

    # Correlation: result[i] = sum_k charge[i+k] * W[k]
    # This gives the wake from particles ahead (larger indices)
    corr_result = correlate(charge, W_causal, mode="full")

    # Extract the portion corresponding to non-negative lags
    # correlate output has length 2*n-1, with zero lag at index n-1
    integrated_wake = corr_result[n - 1 : 2 * n - 1]

    # Correlation includes full W[0] contribution from self
    # Remove it, then add half if self-kick is requested
    W0 = W_causal[0]
    integrated_wake = integrated_wake - charge * W0
    if include_self_kick:
        integrated_wake = integrated_wake + 0.5 * W0 * charge

    if plot:
        self._plot_convolve_density(density, dz, integrated_wake, ax=ax)

    return integrated_wake

impedance

impedance(k: ndarray | float) -> ndarray | complex

Compute the impedance Z(k) analytically.

Sums the contribution from each pseudomode using the analytic Fourier transform.

Parameters:

• k (float or ndarray) –

  Wavenumber [1/m].

Returns:

• Z ( complex or ndarray ) –

  Impedance [Ohm/m].

Source code in beamphysics/wakefields/pseudomode.py

def impedance(self, k: np.ndarray | float) -> np.ndarray | complex:
    """
    Compute the impedance Z(k) analytically.

    Sums the contribution from each pseudomode using the analytic
    Fourier transform.

    Parameters
    ----------
    k : float or np.ndarray
        Wavenumber [1/m].

    Returns
    -------
    Z : complex or np.ndarray
        Impedance [Ohm/m].
    """
    k = np.asarray(k)
    scalar_input = k.ndim == 0
    k = np.atleast_1d(k)

    Z = np.zeros(len(k), dtype=complex)
    for mode in self._modes:
        Z += mode.impedance(k)

    if scalar_input:
        return complex(Z[0])
    return Z

particle_kicks

particle_kicks(z: ndarray, weight: ndarray, include_self_kick: bool = True) -> ndarray

Compute short-range wakefield energy kicks per unit length.

Uses an O(N) single-pass algorithm for each mode by exploiting the exponential form of the pseudomode wakefield.

Parameters:

• z (ndarray) –

  Particle positions [m].

• weight (ndarray) –

  Particle charges [C].

• include_self_kick (bool, default: True ) –

  If True, applies the self-kick. Default is True.

Returns:

• kicks ( ndarray ) –

  Wake-induced energy kick per unit length [eV/m].

Source code in beamphysics/wakefields/pseudomode.py

def particle_kicks(
    self,
    z: np.ndarray,
    weight: np.ndarray,
    include_self_kick: bool = True,
) -> np.ndarray:
    """
    Compute short-range wakefield energy kicks per unit length.

    Uses an O(N) single-pass algorithm for each mode by exploiting
    the exponential form of the pseudomode wakefield.

    Parameters
    ----------
    z : np.ndarray
        Particle positions [m].
    weight : np.ndarray
        Particle charges [C].
    include_self_kick : bool, optional
        If True, applies the self-kick. Default is True.

    Returns
    -------
    kicks : np.ndarray
        Wake-induced energy kick per unit length [eV/m].
    """
    z = np.asarray(z)
    weight = np.asarray(weight)

    if z.shape != weight.shape:
        raise ValueError(
            f"Mismatched shapes: z.shape={z.shape}, weight.shape={weight.shape}"
        )
    if z.ndim != 1:
        raise ValueError("z and weight must be 1D arrays")

    # Sort particles from tail to head
    ix = z.argsort()
    z_sorted = z[ix].copy()
    z_sorted -= z_sorted.max()  # Offset for numerical stability
    weight_sorted = weight[ix]

    n = len(z)
    delta_E = np.zeros(n)

    # Process each mode with O(N) algorithm
    for mode in self._modes:
        s = mode.d + 1j * mode.k
        c = mode.A * np.exp(1j * mode.phi)

        # O(N) accumulator for this mode
        b = 0.0 + 0.0j

        # Iterate from tail to head
        for i in range(n - 1, -1, -1):
            zi = z_sorted[i]
            qi = weight_sorted[i]

            # Kick from trailing particles
            delta_E[i] -= np.imag(c * np.exp(s * zi) * b)

            # Add this particle to accumulator
            b += qi * np.exp(-s * zi)

        if include_self_kick:
            delta_E -= 0.5 * mode.A * weight_sorted * np.sin(mode.phi)

    # Restore original order
    kicks = np.empty_like(delta_E)
    kicks[ix] = delta_E

    return kicks

plot

plot(zmax: float = 0.001, zmin: float = 0, n: int = 200, ax=None)

Plot the pseudomode wakefield.

Source code in beamphysics/wakefields/pseudomode.py

def plot(self, zmax: float = 0.001, zmin: float = 0, n: int = 200, ax=None):
    """Plot the pseudomode wakefield."""
    zlist = np.linspace(zmin, zmax, n)
    Wz = self.wake(-zlist)

    if ax is None:
        _, ax = plt.subplots()
    ax.plot(zlist * 1e6, Wz * 1e-12)
    ax.set_xlabel(r"Distance behind source $|z|$ (µm)")
    ax.set_ylabel(r"$W(z)$ (V/pC/m)")
    title = f"Pseudomode Wakefield ({self.n_modes} mode"
    if self.n_modes > 1:
        title += "s"
    title += ")"
    ax.set_title(title)

plot_impedance

plot_impedance(kmax: float = None, kmin: float = 0, n: int = 200, ax=None)

Plot the impedance over a range of wavenumbers.

Parameters:

• kmax (float, default: None ) –

  Maximum wavenumber [1/m]. If None, uses a sensible default.

• kmin (float, default: 0 ) –

  Minimum wavenumber [1/m]. Default is 0.

• n (int, default: 200 ) –

  Number of points. Default is 200.

• ax (Axes, default: None ) –

  Axes to plot on. If None, creates a new figure.

Source code in beamphysics/wakefields/base.py

def plot_impedance(
    self, kmax: float = None, kmin: float = 0, n: int = 200, ax=None
):
    """
    Plot the impedance over a range of wavenumbers.

    Parameters
    ----------
    kmax : float, optional
        Maximum wavenumber [1/m]. If None, uses a sensible default.
    kmin : float, optional
        Minimum wavenumber [1/m]. Default is 0.
    n : int, optional
        Number of points. Default is 200.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if kmax is None:
        kmax = 1e6  # 1/µm default

    k = np.linspace(kmin, kmax, n)
    Z = self.impedance(k)

    if ax is None:
        _, ax = plt.subplots()
    ax.plot(k * 1e-3, np.real(Z), label=r"Re[$Z$]")
    ax.plot(k * 1e-3, np.imag(Z), label=r"Im[$Z$]")
    ax.set_xlabel(r"$k$ (1/mm)")
    ax.set_ylabel(r"$Z(k)$ (Ω/m)")
    ax.legend()

to_bmad

to_bmad(type: str = 'longitudinal', transverse_dependence: str = 'none') -> str

Format pseudomode parameters as Bmad-compatible strings.

Parameters:

• type (str, default: 'longitudinal' ) –

  Wake type. Default is "longitudinal".

• transverse_dependence (str, default: 'none' ) –

  Transverse dependence. Default is "none".

Returns:

• str –

  Bmad-formatted string (one line per mode).

Source code in beamphysics/wakefields/pseudomode.py

def to_bmad(
    self,
    type: str = "longitudinal",
    transverse_dependence: str = "none",
) -> str:
    """
    Format pseudomode parameters as Bmad-compatible strings.

    Parameters
    ----------
    type : str, optional
        Wake type. Default is "longitudinal".
    transverse_dependence : str, optional
        Transverse dependence. Default is "none".

    Returns
    -------
    str
        Bmad-formatted string (one line per mode).
    """
    lines = [mode.to_bmad(type, transverse_dependence) for mode in self._modes]
    return "\n".join(lines)

wake

wake(z: ndarray | float) -> ndarray | float

Evaluate the pseudomode wakefield at position z.

Parameters:

• z (float or ndarray) –

  Longitudinal position [m]. Negative z is behind the source.

Returns:

• W ( float or ndarray ) –

  Wakefield value [V/C/m]. Returns 0 for z > 0.

Source code in beamphysics/wakefields/pseudomode.py

def wake(self, z: np.ndarray | float) -> np.ndarray | float:
    """
    Evaluate the pseudomode wakefield at position z.

    Parameters
    ----------
    z : float or np.ndarray
        Longitudinal position [m]. Negative z is behind the source.

    Returns
    -------
    W : float or np.ndarray
        Wakefield value [V/C/m]. Returns 0 for z > 0.
    """
    z = np.asarray(z)
    scalar_input = z.ndim == 0
    z = np.atleast_1d(z)

    out = np.zeros_like(z, dtype=float)
    mask = z <= 0
    z_masked = z[mask]

    # Sum contributions from all modes
    for mode in self._modes:
        out[mask] += mode(z_masked)

    if scalar_input:
        return float(out[0])
    return out

ImpedanceWakefield

ImpedanceWakefield(impedance_func: Callable[[ndarray], ndarray], wakefield_func: Callable[[ndarray], ndarray] = None, k_max: float = 10000000.0)

Bases: WakefieldBase

              flowchart TD
              beamphysics.wakefields.ImpedanceWakefield[ImpedanceWakefield]
              beamphysics.wakefields.base.WakefieldBase[WakefieldBase]

                              beamphysics.wakefields.base.WakefieldBase --> beamphysics.wakefields.ImpedanceWakefield
                


              click beamphysics.wakefields.ImpedanceWakefield href "" "beamphysics.wakefields.ImpedanceWakefield"
              click beamphysics.wakefields.base.WakefieldBase href "" "beamphysics.wakefields.base.WakefieldBase"
            

Wakefield defined through its longitudinal impedance Z(k).

Uses FFT-based convolution in frequency domain for efficient wake potential calculations.

Parameters:

• impedance_func (callable) –

  Function that returns complex impedance Z(k) [Ohm/m] for wave number k [1/m]. Should handle arrays.

• wakefield_func (callable, default: None ) –

  Function that returns wakefield W(z) [V/C/m] for position z [m]. If not provided, uses numerical cosine transform of impedance.

• k_max (float, default: 10000000.0 ) –

  Maximum wavenumber for integration [1/m]. Default is 1e7.

Notes

The wakefield W(z) is defined for z ≤ 0 (trailing particles), with W(z) = 0 for z > 0 (causality). The impedance Z(k) is related to the wakefield by:

\[Z(k) = \frac{1}{c} \int_{-\infty}^{0} W(z) e^{-ikz} dz\]

And the inverse:

\[W(z) = \frac{2c}{\pi} \int_0^{\infty} \text{Re}[Z(k)] \cos(kz) dk \quad (z \le 0)\]

Examples:

::

def my_impedance(k):
    return 100 / (1 + 1j * k * 1e-6)  # Simple resonator
wake = ImpedanceWakefield(my_impedance)

Subclassed by:

• API Wakefields  ResistiveWallWakefield

Methods:

• __call__ –

  Evaluate the wakefield at position z (convenience method).

• convolve_density –

  Compute integrated wakefield from a charge density distribution.

• impedance –

  Evaluate the impedance at wavenumber k.

• particle_kicks –

  Compute wakefield-induced energy kicks per unit length.

• plot –

  Plot the wakefield over a range of z values.

• plot_impedance –

  Plot the impedance over a range of wavenumbers.

• wake –

  Evaluate the wakefield at position z.

Attributes:

• k_max (float) –

  Maximum wavenumber for integration.

Source code in beamphysics/wakefields/impedance.py

def __init__(
    self,
    impedance_func: Callable[[np.ndarray], np.ndarray],
    wakefield_func: Callable[[np.ndarray], np.ndarray] = None,
    k_max: float = 1e7,
) -> None:
    self._impedance_func = impedance_func
    self._wakefield_func = wakefield_func
    self._k_max = k_max

k_max property writable

k_max: float

Maximum wavenumber for integration.

__call__

__call__(z: ndarray | float) -> ndarray | float

Evaluate the wakefield at position z (convenience method).

This is equivalent to calling wake(z).

Source code in beamphysics/wakefields/base.py

def __call__(self, z: np.ndarray | float) -> np.ndarray | float:
    """
    Evaluate the wakefield at position z (convenience method).

    This is equivalent to calling wake(z).
    """
    return self.wake(z)

convolve_density

convolve_density(density: ndarray, dz: float, offset: float = 0, include_self_kick: bool = True, plot: bool = False, ax=None) -> ndarray

Compute integrated wakefield from a charge density distribution.

This computes the causal wake potential:

\[V(z_i) = \sum_{j>i} Q_j \cdot W(z_j - z_i) + \frac{1}{2} Q_i \cdot W(0)\]

where only particles ahead (larger z index) contribute to the wake felt by each particle. This is mathematically a correlation, not a convolution, and is computed efficiently using FFT.

Parameters:

• density (ndarray) –

  Charge density array [C/m]. Index 0 is the tail (back), index n-1 is the head (front).

• dz (float) –

  Grid spacing [m].

• offset (float, default: 0 ) –

  Offset for the z coordinate [m]. Default is 0.

• include_self_kick (bool, default: True ) –

  Whether to include the half self-kick term. Default is True.

• plot (bool, default: False ) –

  If True, plot the density profile and wake potential. Default is False.

• ax (Axes, default: None ) –

  Axes to plot on. If None and plot=True, creates a new figure.

Returns:

• integrated_wake ( ndarray ) –

  Integrated longitudinal wakefield [V].

Source code in beamphysics/wakefields/base.py

def convolve_density(
    self,
    density: np.ndarray,
    dz: float,
    offset: float = 0,
    include_self_kick: bool = True,
    plot: bool = False,
    ax=None,
) -> np.ndarray:
    """
    Compute integrated wakefield from a charge density distribution.

    This computes the causal wake potential:

    $$V(z_i) = \\sum_{j>i} Q_j \\cdot W(z_j - z_i) + \\frac{1}{2} Q_i \\cdot W(0)$$

    where only particles ahead (larger z index) contribute to the wake
    felt by each particle. This is mathematically a correlation, not
    a convolution, and is computed efficiently using FFT.

    Parameters
    ----------
    density : np.ndarray
        Charge density array [C/m]. Index 0 is the tail (back),
        index n-1 is the head (front).
    dz : float
        Grid spacing [m].
    offset : float, optional
        Offset for the z coordinate [m]. Default is 0.
    include_self_kick : bool, optional
        Whether to include the half self-kick term. Default is True.
    plot : bool, optional
        If True, plot the density profile and wake potential. Default is False.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None and plot=True, creates a new figure.

    Returns
    -------
    integrated_wake : np.ndarray
        Integrated longitudinal wakefield [V].
    """
    from scipy.signal import correlate

    n = len(density)
    charge = density * dz  # Convert to charge per bin [C]

    # Build causal wake array: W[k] = wake at lag k*dz
    # W[0] = W(0), W[1] = W(-dz), W[2] = W(-2*dz), ...
    z_wake = -np.arange(n) * dz + offset
    W_causal = self.wake(z_wake)

    # Correlation: result[i] = sum_k charge[i+k] * W[k]
    # This gives the wake from particles ahead (larger indices)
    corr_result = correlate(charge, W_causal, mode="full")

    # Extract the portion corresponding to non-negative lags
    # correlate output has length 2*n-1, with zero lag at index n-1
    integrated_wake = corr_result[n - 1 : 2 * n - 1]

    # Correlation includes full W[0] contribution from self
    # Remove it, then add half if self-kick is requested
    W0 = W_causal[0]
    integrated_wake = integrated_wake - charge * W0
    if include_self_kick:
        integrated_wake = integrated_wake + 0.5 * W0 * charge

    if plot:
        self._plot_convolve_density(density, dz, integrated_wake, ax=ax)

    return integrated_wake

impedance

impedance(k: ndarray | float) -> ndarray | complex

Evaluate the impedance at wavenumber k.

Parameters:

• k (float or ndarray) –

  Wavenumber [1/m].

Returns:

• Z ( complex or ndarray ) –

  Impedance [Ohm/m].

Source code in beamphysics/wakefields/impedance.py

def impedance(self, k: np.ndarray | float) -> np.ndarray | complex:
    """
    Evaluate the impedance at wavenumber k.

    Parameters
    ----------
    k : float or np.ndarray
        Wavenumber [1/m].

    Returns
    -------
    Z : complex or np.ndarray
        Impedance [Ohm/m].
    """
    return self._impedance_func(k)

particle_kicks

particle_kicks(z: ndarray, weight: ndarray, include_self_kick: bool = True, n_bins: int = None) -> ndarray

Compute wakefield-induced energy kicks per unit length.

Uses FFT-based density convolution with interpolation back to particle positions. This is O(N + M log M) where N is the number of particles and M is the number of grid points.

Parameters:

• z (ndarray) –

  Particle positions [m].

• weight (ndarray) –

  Particle charges [C].

• include_self_kick (bool, default: True ) –

  Whether to include the self-kick term. Default is True.

• n_bins (int, default: None ) –

  Number of bins for the density grid. Default is max(100, N//10).

Returns:

• kicks ( ndarray ) –

  Energy kick per unit length at each particle [eV/m].

Source code in beamphysics/wakefields/impedance.py

def particle_kicks(
    self,
    z: np.ndarray,
    weight: np.ndarray,
    include_self_kick: bool = True,
    n_bins: int = None,
) -> np.ndarray:
    """
    Compute wakefield-induced energy kicks per unit length.

    Uses FFT-based density convolution with interpolation back to
    particle positions. This is O(N + M log M) where N is the number
    of particles and M is the number of grid points.

    Parameters
    ----------
    z : np.ndarray
        Particle positions [m].
    weight : np.ndarray
        Particle charges [C].
    include_self_kick : bool, optional
        Whether to include the self-kick term. Default is True.
    n_bins : int, optional
        Number of bins for the density grid. Default is max(100, N//10).

    Returns
    -------
    kicks : np.ndarray
        Energy kick per unit length at each particle [eV/m].
    """
    from scipy.interpolate import interp1d

    z = np.asarray(z)
    weight = np.asarray(weight)
    n = len(z)

    if n_bins is None:
        n_bins = max(100, n // 10)

    # Bin particles into density distribution
    z_min, z_max = z.min(), z.max()
    z_range = z_max - z_min
    if z_range == 0:
        z_range = 1e-6  # Avoid division by zero for single particle

    # Add padding to avoid edge effects
    padding = 0.1 * z_range
    z_min_padded = z_min - padding
    z_max_padded = z_max + padding

    dz = (z_max_padded - z_min_padded) / n_bins
    z_grid = np.linspace(z_min_padded + dz / 2, z_max_padded - dz / 2, n_bins)

    # Create density histogram (charge per bin)
    density, bin_edges = np.histogram(
        z, bins=n_bins, range=(z_min_padded, z_max_padded), weights=weight
    )
    # Convert to charge density [C/m]
    density = density / dz

    # Compute wake potential via FFT convolution
    wake_potential = self.convolve_density(
        density, dz, include_self_kick=include_self_kick
    )

    # Interpolate wake potential to particle positions
    interp = interp1d(
        z_grid, wake_potential, kind="linear", bounds_error=False, fill_value=0.0
    )
    kicks = -interp(z)  # Negative sign: energy loss for trailing particles

    return kicks

plot

plot(zmax: float = None, zmin: float = 0, n: int = 200, ax=None)

Plot the wakefield over a range of z values.

Parameters:

• zmax (float, default: None ) –

  Maximum trailing distance [m]. If None, uses a sensible default.

• zmin (float, default: 0 ) –

  Minimum trailing distance [m]. Default is 0.

• n (int, default: 200 ) –

  Number of points. Default is 200.

• ax (Axes, default: None ) –

  Axes to plot on. If None, creates a new figure.

Source code in beamphysics/wakefields/base.py

def plot(self, zmax: float = None, zmin: float = 0, n: int = 200, ax=None):
    """
    Plot the wakefield over a range of z values.

    Parameters
    ----------
    zmax : float, optional
        Maximum trailing distance [m]. If None, uses a sensible default.
    zmin : float, optional
        Minimum trailing distance [m]. Default is 0.
    n : int, optional
        Number of points. Default is 200.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if zmax is None:
        zmax = 1e-3  # 1 mm default

    zlist = np.linspace(zmin, zmax, n)
    Wz = self.wake(-zlist)  # Negative z for trailing particles

    if ax is None:
        _, ax = plt.subplots()
    ax.plot(zlist * 1e6, Wz * 1e-12)
    ax.set_xlabel(r"Distance behind source $|z|$ (µm)")
    ax.set_ylabel(r"$W(z)$ (V/pC/m)")

plot_impedance

plot_impedance(kmax: float = None, kmin: float = 0, n: int = 200, ax=None)

Plot the impedance over a range of wavenumbers.

Parameters:

• kmax (float, default: None ) –

  Maximum wavenumber [1/m]. If None, uses a sensible default.

• kmin (float, default: 0 ) –

  Minimum wavenumber [1/m]. Default is 0.

• n (int, default: 200 ) –

  Number of points. Default is 200.

• ax (Axes, default: None ) –

  Axes to plot on. If None, creates a new figure.

Source code in beamphysics/wakefields/base.py

def plot_impedance(
    self, kmax: float = None, kmin: float = 0, n: int = 200, ax=None
):
    """
    Plot the impedance over a range of wavenumbers.

    Parameters
    ----------
    kmax : float, optional
        Maximum wavenumber [1/m]. If None, uses a sensible default.
    kmin : float, optional
        Minimum wavenumber [1/m]. Default is 0.
    n : int, optional
        Number of points. Default is 200.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if kmax is None:
        kmax = 1e6  # 1/µm default

    k = np.linspace(kmin, kmax, n)
    Z = self.impedance(k)

    if ax is None:
        _, ax = plt.subplots()
    ax.plot(k * 1e-3, np.real(Z), label=r"Re[$Z$]")
    ax.plot(k * 1e-3, np.imag(Z), label=r"Im[$Z$]")
    ax.set_xlabel(r"$k$ (1/mm)")
    ax.set_ylabel(r"$Z(k)$ (Ω/m)")
    ax.legend()

wake

wake(z: ndarray | float, n_fft: int = 4096) -> ndarray | float

Evaluate the wakefield at position z.

If a wakefield function was provided, uses it directly. Otherwise, computes via cosine transform of Re[Z(k)]: - Scalars: numerical quadrature (accurate) - Arrays: FFT with interpolation (fast)

Parameters:

• z (float or ndarray) –

  Longitudinal position [m].

• n_fft (int, default: 4096 ) –

  Number of FFT points for array evaluation. Default is 4096.

Returns:

• W ( float or ndarray ) –

  Wakefield value [V/C/m]. Returns 0 for z > 0.

Source code in beamphysics/wakefields/impedance.py

def wake(self, z: np.ndarray | float, n_fft: int = 4096) -> np.ndarray | float:
    """
    Evaluate the wakefield at position z.

    If a wakefield function was provided, uses it directly.
    Otherwise, computes via cosine transform of Re[Z(k)]:
    - Scalars: numerical quadrature (accurate)
    - Arrays: FFT with interpolation (fast)

    Parameters
    ----------
    z : float or np.ndarray
        Longitudinal position [m].
    n_fft : int, optional
        Number of FFT points for array evaluation. Default is 4096.

    Returns
    -------
    W : float or np.ndarray
        Wakefield value [V/C/m]. Returns 0 for z > 0.
    """
    if self._wakefield_func is not None:
        return self._wakefield_func(z)

    z = np.asarray(z)
    scalar_input = z.ndim == 0
    z = np.atleast_1d(z)

    if scalar_input:
        # Single point: use quadrature for accuracy
        from scipy.integrate import quad

        z_val = float(z[0])
        if z_val > 0:
            return 0.0

        def integrand(k):
            return np.real(self._impedance_func(k)) * np.cos(k * z_val)

        result, _ = quad(integrand, 0, self._k_max, limit=200)
        return (2 * c_light / np.pi) * result
    else:
        # Array: use FFT for speed
        from scipy.fft import irfft
        from scipy.interpolate import interp1d

        k_max = self._k_max
        dk = k_max / (n_fft - 1)
        k_grid = np.linspace(0, k_max, n_fft)

        Zk_grid = self._impedance_func(k_grid)
        ReZ = np.real(Zk_grid)

        n_full = 2 * (n_fft - 1)
        dz = 2 * np.pi / (n_full * dk)
        z_grid = np.arange(n_full) * dz

        W_grid = irfft(ReZ, n=n_full) * n_full * dk * (c_light / np.pi)

        interp = interp1d(
            z_grid, W_grid, kind="cubic", bounds_error=False, fill_value=0.0
        )

        result = interp(-z)
        result = np.where(z > 0, 0.0, result)

        return result

TabularWakefield

TabularWakefield(z: ndarray, W: ndarray, fill_value: float = 0.0, kind: str = 'cubic')

Bases: WakefieldBase

              flowchart TD
              beamphysics.wakefields.TabularWakefield[TabularWakefield]
              beamphysics.wakefields.base.WakefieldBase[WakefieldBase]

                              beamphysics.wakefields.base.WakefieldBase --> beamphysics.wakefields.TabularWakefield
                


              click beamphysics.wakefields.TabularWakefield href "" "beamphysics.wakefields.TabularWakefield"
              click beamphysics.wakefields.base.WakefieldBase href "" "beamphysics.wakefields.base.WakefieldBase"
            

Wakefield defined by user-supplied tabular data with interpolation.

Uses cubic spline interpolation to evaluate the wakefield at arbitrary positions between the supplied data points.

Parameters:

• z (ndarray) –

  Longitudinal positions [m]. Should be negative (behind source) and sorted in ascending order.

• W (ndarray) –

  Wakefield values [V/C/m] at each z position.

• fill_value (float, default: 0.0 ) –

  Value to return outside the interpolation range. Default is 0.

• kind (str, default: 'cubic' ) –

  Interpolation method. Default is 'cubic'.

Examples:

::

z_data = -np.linspace(1e-6, 1e-3, 100)
W_data = 1e15 * np.exp(z_data / 100e-6) * np.sin(1e5 * z_data)
wake = TabularWakefield(z_data, W_data)
wake.wake(-50e-6)  # Interpolated wake at 50 µm behind source

Methods:

• __call__ –

  Evaluate the wakefield at position z (convenience method).

• convolve_density –

  Compute integrated wakefield from a charge density distribution.

• impedance –

  Compute the impedance Z(k) via numerical FFT.

• particle_kicks –

  Compute wakefield-induced energy kicks per unit length.

• plot –

  Plot the wakefield over a range of z values.

• plot_impedance –

  Plot the impedance over a range of wavenumbers.

• wake –

  Evaluate the wakefield at position z using interpolation.

Attributes:

• W_data (ndarray) –

  Return the W data points.

• z_data (ndarray) –

  Return the z data points.

Source code in beamphysics/wakefields/tabular.py

def __init__(
    self,
    z: np.ndarray,
    W: np.ndarray,
    fill_value: float = 0.0,
    kind: str = "cubic",
) -> None:
    z = np.asarray(z)
    W = np.asarray(W)

    if z.shape != W.shape:
        raise ValueError(f"Shape mismatch: z.shape={z.shape}, W.shape={W.shape}")
    if len(z) < 4:
        raise ValueError("Need at least 4 points for cubic interpolation")

    # Store data
    self._z = z
    self._W = W
    self._fill_value = fill_value

    # Create interpolator
    self._interp = interp1d(
        z,
        W,
        kind=kind,
        bounds_error=False,
        fill_value=fill_value,
    )

W_data property

W_data: ndarray

Return the W data points.

z_data property

z_data: ndarray

Return the z data points.

__call__

__call__(z: ndarray | float) -> ndarray | float

Evaluate the wakefield at position z (convenience method).

This is equivalent to calling wake(z).

Source code in beamphysics/wakefields/base.py

def __call__(self, z: np.ndarray | float) -> np.ndarray | float:
    """
    Evaluate the wakefield at position z (convenience method).

    This is equivalent to calling wake(z).
    """
    return self.wake(z)

convolve_density

convolve_density(density: ndarray, dz: float, offset: float = 0, include_self_kick: bool = True, plot: bool = False, ax=None) -> ndarray

Compute integrated wakefield from a charge density distribution.

This computes the causal wake potential:

\[V(z_i) = \sum_{j>i} Q_j \cdot W(z_j - z_i) + \frac{1}{2} Q_i \cdot W(0)\]

where only particles ahead (larger z index) contribute to the wake felt by each particle. This is mathematically a correlation, not a convolution, and is computed efficiently using FFT.

Parameters:

• density (ndarray) –

  Charge density array [C/m]. Index 0 is the tail (back), index n-1 is the head (front).

• dz (float) –

  Grid spacing [m].

• offset (float, default: 0 ) –

  Offset for the z coordinate [m]. Default is 0.

• include_self_kick (bool, default: True ) –

  Whether to include the half self-kick term. Default is True.

• plot (bool, default: False ) –

  If True, plot the density profile and wake potential. Default is False.

• ax (Axes, default: None ) –

  Axes to plot on. If None and plot=True, creates a new figure.

Returns:

• integrated_wake ( ndarray ) –

  Integrated longitudinal wakefield [V].

Source code in beamphysics/wakefields/base.py

def convolve_density(
    self,
    density: np.ndarray,
    dz: float,
    offset: float = 0,
    include_self_kick: bool = True,
    plot: bool = False,
    ax=None,
) -> np.ndarray:
    """
    Compute integrated wakefield from a charge density distribution.

    This computes the causal wake potential:

    $$V(z_i) = \\sum_{j>i} Q_j \\cdot W(z_j - z_i) + \\frac{1}{2} Q_i \\cdot W(0)$$

    where only particles ahead (larger z index) contribute to the wake
    felt by each particle. This is mathematically a correlation, not
    a convolution, and is computed efficiently using FFT.

    Parameters
    ----------
    density : np.ndarray
        Charge density array [C/m]. Index 0 is the tail (back),
        index n-1 is the head (front).
    dz : float
        Grid spacing [m].
    offset : float, optional
        Offset for the z coordinate [m]. Default is 0.
    include_self_kick : bool, optional
        Whether to include the half self-kick term. Default is True.
    plot : bool, optional
        If True, plot the density profile and wake potential. Default is False.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None and plot=True, creates a new figure.

    Returns
    -------
    integrated_wake : np.ndarray
        Integrated longitudinal wakefield [V].
    """
    from scipy.signal import correlate

    n = len(density)
    charge = density * dz  # Convert to charge per bin [C]

    # Build causal wake array: W[k] = wake at lag k*dz
    # W[0] = W(0), W[1] = W(-dz), W[2] = W(-2*dz), ...
    z_wake = -np.arange(n) * dz + offset
    W_causal = self.wake(z_wake)

    # Correlation: result[i] = sum_k charge[i+k] * W[k]
    # This gives the wake from particles ahead (larger indices)
    corr_result = correlate(charge, W_causal, mode="full")

    # Extract the portion corresponding to non-negative lags
    # correlate output has length 2*n-1, with zero lag at index n-1
    integrated_wake = corr_result[n - 1 : 2 * n - 1]

    # Correlation includes full W[0] contribution from self
    # Remove it, then add half if self-kick is requested
    W0 = W_causal[0]
    integrated_wake = integrated_wake - charge * W0
    if include_self_kick:
        integrated_wake = integrated_wake + 0.5 * W0 * charge

    if plot:
        self._plot_convolve_density(density, dz, integrated_wake, ax=ax)

    return integrated_wake

impedance

impedance(k: ndarray | float) -> ndarray | complex

Compute the impedance Z(k) via numerical FFT.

Uses FFT to compute the Fourier transform of the tabular wake data.

Parameters:

• k (float or ndarray) –

  Wavenumber [1/m].

Returns:

• Z ( complex or ndarray ) –

  Impedance [Ohm/m].

Source code in beamphysics/wakefields/tabular.py

def impedance(self, k: np.ndarray | float) -> np.ndarray | complex:
    """
    Compute the impedance Z(k) via numerical FFT.

    Uses FFT to compute the Fourier transform of the tabular wake data.

    Parameters
    ----------
    k : float or np.ndarray
        Wavenumber [1/m].

    Returns
    -------
    Z : complex or np.ndarray
        Impedance [Ohm/m].
    """
    k = np.asarray(k)
    scalar_input = k.ndim == 0
    k = np.atleast_1d(k)

    # Use the stored data for FFT
    z_data = self._z
    W_data = self._W

    # Sort by z (ascending)
    sort_idx = np.argsort(z_data)
    z_sorted = z_data[sort_idx]
    W_sorted = W_data[sort_idx]

    # Z(k) = (1/c) * integral of W(z) * exp(-ikz) dz
    # For each k, compute numerical integral
    Z = np.zeros(len(k), dtype=complex)
    for i, ki in enumerate(k):
        integrand = W_sorted * np.exp(-1j * ki * z_sorted)
        Z[i] = np.trapezoid(integrand, z_sorted) / c_light

    if scalar_input:
        return complex(Z[0])
    return Z

particle_kicks

particle_kicks(z: ndarray, weight: ndarray, include_self_kick: bool = True) -> ndarray

Compute wakefield-induced energy kicks per unit length.

This is a default O(N²) implementation. Subclasses may override with more efficient algorithms.

Parameters:

• z (ndarray) –

  Particle positions [m].

• weight (ndarray) –

  Particle charges [C].

• include_self_kick (bool, default: True ) –

  Whether to include the self-kick term. Default is True.

Returns:

• kicks ( ndarray ) –

  Energy kick per unit length at each particle [eV/m].

Source code in beamphysics/wakefields/base.py

def particle_kicks(
    self,
    z: np.ndarray,
    weight: np.ndarray,
    include_self_kick: bool = True,
) -> np.ndarray:
    """
    Compute wakefield-induced energy kicks per unit length.

    This is a default O(N²) implementation. Subclasses may override
    with more efficient algorithms.

    Parameters
    ----------
    z : np.ndarray
        Particle positions [m].
    weight : np.ndarray
        Particle charges [C].
    include_self_kick : bool, optional
        Whether to include the self-kick term. Default is True.

    Returns
    -------
    kicks : np.ndarray
        Energy kick per unit length at each particle [eV/m].
    """
    z = np.asarray(z)
    weight = np.asarray(weight)
    n = len(z)
    kicks = np.zeros(n)

    # O(N²) algorithm: sum contribution from all particles
    for i in range(n):
        for j in range(n):
            if i == j:
                if include_self_kick:
                    # Self-kick uses W(0⁻) / 2
                    kicks[i] -= 0.5 * weight[j] * self._self_kick_value()
            else:
                dz = z[i] - z[j]
                if dz < 0:  # j is ahead of i, so i feels wake from j
                    kicks[i] -= weight[j] * self.wake(dz)

    return kicks

plot

plot(zmax: float = None, zmin: float = 0, n: int = 200, ax=None)

Plot the wakefield over a range of z values.

Parameters:

• zmax (float, default: None ) –

  Maximum trailing distance [m]. If None, uses a sensible default.

• zmin (float, default: 0 ) –

  Minimum trailing distance [m]. Default is 0.

• n (int, default: 200 ) –

  Number of points. Default is 200.

• ax (Axes, default: None ) –

  Axes to plot on. If None, creates a new figure.

Source code in beamphysics/wakefields/base.py

def plot(self, zmax: float = None, zmin: float = 0, n: int = 200, ax=None):
    """
    Plot the wakefield over a range of z values.

    Parameters
    ----------
    zmax : float, optional
        Maximum trailing distance [m]. If None, uses a sensible default.
    zmin : float, optional
        Minimum trailing distance [m]. Default is 0.
    n : int, optional
        Number of points. Default is 200.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if zmax is None:
        zmax = 1e-3  # 1 mm default

    zlist = np.linspace(zmin, zmax, n)
    Wz = self.wake(-zlist)  # Negative z for trailing particles

    if ax is None:
        _, ax = plt.subplots()
    ax.plot(zlist * 1e6, Wz * 1e-12)
    ax.set_xlabel(r"Distance behind source $|z|$ (µm)")
    ax.set_ylabel(r"$W(z)$ (V/pC/m)")

plot_impedance

plot_impedance(kmax: float = None, kmin: float = 0, n: int = 200, ax=None)

Plot the impedance over a range of wavenumbers.

Parameters:

• kmax (float, default: None ) –

  Maximum wavenumber [1/m]. If None, uses a sensible default.

• kmin (float, default: 0 ) –

  Minimum wavenumber [1/m]. Default is 0.

• n (int, default: 200 ) –

  Number of points. Default is 200.

• ax (Axes, default: None ) –

  Axes to plot on. If None, creates a new figure.

Source code in beamphysics/wakefields/base.py

def plot_impedance(
    self, kmax: float = None, kmin: float = 0, n: int = 200, ax=None
):
    """
    Plot the impedance over a range of wavenumbers.

    Parameters
    ----------
    kmax : float, optional
        Maximum wavenumber [1/m]. If None, uses a sensible default.
    kmin : float, optional
        Minimum wavenumber [1/m]. Default is 0.
    n : int, optional
        Number of points. Default is 200.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if kmax is None:
        kmax = 1e6  # 1/µm default

    k = np.linspace(kmin, kmax, n)
    Z = self.impedance(k)

    if ax is None:
        _, ax = plt.subplots()
    ax.plot(k * 1e-3, np.real(Z), label=r"Re[$Z$]")
    ax.plot(k * 1e-3, np.imag(Z), label=r"Im[$Z$]")
    ax.set_xlabel(r"$k$ (1/mm)")
    ax.set_ylabel(r"$Z(k)$ (Ω/m)")
    ax.legend()

wake

wake(z: ndarray | float) -> ndarray | float

Evaluate the wakefield at position z using interpolation.

Parameters:

• z (float or ndarray) –

  Longitudinal position [m].

Returns:

• W ( float or ndarray ) –

  Interpolated wakefield value [V/C/m]. Returns fill_value outside the data range, and 0 for z > 0 (causality).

Source code in beamphysics/wakefields/tabular.py

def wake(self, z: np.ndarray | float) -> np.ndarray | float:
    """
    Evaluate the wakefield at position z using interpolation.

    Parameters
    ----------
    z : float or np.ndarray
        Longitudinal position [m].

    Returns
    -------
    W : float or np.ndarray
        Interpolated wakefield value [V/C/m]. Returns fill_value
        outside the data range, and 0 for z > 0 (causality).
    """
    z = np.asarray(z)
    scalar_input = z.ndim == 0
    z = np.atleast_1d(z)

    # Apply causality
    result = np.where(z > 0, 0.0, self._interp(z))

    if scalar_input:
        return float(result[0])
    return result

Pseudomode dataclass

Pseudomode(A: float, d: float, k: float, phi: float)

Single pseudomode parameters for a damped sinusoidal wakefield component.

Represents one term in the sum: W(z) = A·exp(d·z)·sin(k·z + φ)

Parameters:

• A (float) –

  Amplitude coefficient [V/C/m].

• d (float) –

  Exponential decay rate [1/m].

• k (float) –

  Oscillation wavenumber [1/m].

• phi (float) –

  Phase offset [rad].

Methods:

• __call__ –

  Evaluate this mode at position z (array input assumed).

• impedance –

  Compute the impedance Z(k) for this pseudomode analytically.

• to_bmad –

  Format as Bmad-compatible string.

__call__

__call__(z: ndarray) -> ndarray

Evaluate this mode at position z (array input assumed).

Source code in beamphysics/wakefields/pseudomode.py

def __call__(self, z: np.ndarray) -> np.ndarray:
    """Evaluate this mode at position z (array input assumed)."""
    return self.A * np.exp(self.d * z) * np.sin(self.k * z + self.phi)

impedance

impedance(k: ndarray | float) -> ndarray | complex

Compute the impedance Z(k) for this pseudomode analytically.

For a pseudomode wakefield defined for z ≤ 0:

\[W(z) = A \cdot e^{dz} \cdot \sin(k_0 z + \phi)\]

The impedance is defined as:

\[Z(k) = \frac{1}{c} \int_{-\infty}^{0} W(z) \cdot e^{-ikz} \, dz\]

Derivation:

Substituting the wakefield:

\[Z(k) = \frac{A}{c} \int_{-\infty}^{0} e^{dz} \sin(k_0 z + \phi) e^{-ikz} \, dz\]

Using \(\sin(\theta) = \frac{e^{i\theta} - e^{-i\theta}}{2i}\):

\[Z(k) = \frac{A}{2ic} \left[ e^{i\phi} \int_{-\infty}^{0} e^{(d + i(k_0 - k))z} dz - e^{-i\phi} \int_{-\infty}^{0} e^{(d - i(k_0 + k))z} dz \right]\]

Since \(d > 0\), both integrals converge:

\[\int_{-\infty}^{0} e^{az} dz = \frac{1}{a} \quad \text{for } \text{Re}(a) > 0\]

Evaluating:

\[Z(k) = \frac{A}{2ic} \left[ \frac{e^{i\phi}}{d + i(k_0 - k)} - \frac{e^{-i\phi}}{d - i(k_0 + k)} \right]\]

Mathematica verification:

(* Define the wakefield and compute impedance numerically *)
W[z_, A_, d_, k0_, phi_] := A Exp[d z] Sin[k0 z + phi]

(* Analytical result *)
Zanalytic[k_, A_, d_, k0_, phi_] :=
    A/(2 I c) (Exp[I phi]/(d + I (k0 - k)) - Exp[-I phi]/(d - I (k0 + k)))

(* Numerical integration (should match) *)
Znumeric[k_, A_, d_, k0_, phi_] :=
    1/c NIntegrate[W[z, A, d, k0, phi] Exp[-I k z], {z, -Infinity, 0}]

(* Test with sample values *)
c = 299792458;
{A, d, k0, phi} = {1*^15, 1*^4, 1*^5, Pi/4};
ktest = 5*^4;
{Zanalytic[ktest, A, d, k0, phi], Znumeric[ktest, A, d, k0, phi]}
(* Both should give the same complex number *)

Parameters:

• k (float or ndarray) –

  Wavenumber [1/m].

Returns:

• Z ( complex or ndarray ) –

  Impedance [Ohm/m].

Source code in beamphysics/wakefields/pseudomode.py

def impedance(self, k: np.ndarray | float) -> np.ndarray | complex:
    """
    Compute the impedance Z(k) for this pseudomode analytically.

    For a pseudomode wakefield defined for z ≤ 0:

    $$W(z) = A \\cdot e^{dz} \\cdot \\sin(k_0 z + \\phi)$$

    The impedance is defined as:

    $$Z(k) = \\frac{1}{c} \\int_{-\\infty}^{0} W(z) \\cdot e^{-ikz} \\, dz$$

    **Derivation:**

    Substituting the wakefield:

    $$Z(k) = \\frac{A}{c} \\int_{-\\infty}^{0} e^{dz} \\sin(k_0 z + \\phi) e^{-ikz} \\, dz$$

    Using $\\sin(\\theta) = \\frac{e^{i\\theta} - e^{-i\\theta}}{2i}$:

    $$Z(k) = \\frac{A}{2ic} \\left[
        e^{i\\phi} \\int_{-\\infty}^{0} e^{(d + i(k_0 - k))z} dz
        - e^{-i\\phi} \\int_{-\\infty}^{0} e^{(d - i(k_0 + k))z} dz
    \\right]$$

    Since $d > 0$, both integrals converge:

    $$\\int_{-\\infty}^{0} e^{az} dz = \\frac{1}{a} \\quad \\text{for } \\text{Re}(a) > 0$$

    Evaluating:

    $$Z(k) = \\frac{A}{2ic} \\left[
        \\frac{e^{i\\phi}}{d + i(k_0 - k)}
        - \\frac{e^{-i\\phi}}{d - i(k_0 + k)}
    \\right]$$

    **Mathematica verification:**

    ```mathematica
    (* Define the wakefield and compute impedance numerically *)
    W[z_, A_, d_, k0_, phi_] := A Exp[d z] Sin[k0 z + phi]

    (* Analytical result *)
    Zanalytic[k_, A_, d_, k0_, phi_] :=
        A/(2 I c) (Exp[I phi]/(d + I (k0 - k)) - Exp[-I phi]/(d - I (k0 + k)))

    (* Numerical integration (should match) *)
    Znumeric[k_, A_, d_, k0_, phi_] :=
        1/c NIntegrate[W[z, A, d, k0, phi] Exp[-I k z], {z, -Infinity, 0}]

    (* Test with sample values *)
    c = 299792458;
    {A, d, k0, phi} = {1*^15, 1*^4, 1*^5, Pi/4};
    ktest = 5*^4;
    {Zanalytic[ktest, A, d, k0, phi], Znumeric[ktest, A, d, k0, phi]}
    (* Both should give the same complex number *)
    ```

    Parameters
    ----------
    k : float or np.ndarray
        Wavenumber [1/m].

    Returns
    -------
    Z : complex or np.ndarray
        Impedance [Ohm/m].
    """
    k = np.asarray(k)
    scalar_input = k.ndim == 0
    k = np.atleast_1d(k)

    d = self.d
    k0 = self.k
    A = self.A
    phi = self.phi

    # Z(k) = A/(2ic) * [e^{iφ}/(d + i(k₀-k)) - e^{-iφ}/(d - i(k₀+k))]
    denom1 = d + 1j * (k0 - k)
    denom2 = d - 1j * (k0 + k)

    Z = (A / (2j * c_light)) * (
        np.exp(1j * phi) / denom1 - np.exp(-1j * phi) / denom2
    )

    if scalar_input:
        return complex(Z[0])
    return Z

to_bmad

to_bmad(type: str = 'longitudinal', transverse_dependence: str = 'none') -> str

Format as Bmad-compatible string.

Source code in beamphysics/wakefields/pseudomode.py

def to_bmad(
    self,
    type: str = "longitudinal",
    transverse_dependence: str = "none",
) -> str:
    """Format as Bmad-compatible string."""
    return (
        f"{type} = {{{self.A}, {self.d}, {self.k}, "
        f"{self.phi / (2 * np.pi)}, {transverse_dependence}}}"
    )

Low-level Functions

longitudinal_impedance_round

longitudinal_impedance_round(k: float | ndarray, a: float, sigma0: float, ctau: float) -> complex | ndarray

Compute longitudinal impedance Z(k) for round (circular pipe) geometry.

Uses closed-form expression (SLAC-PUB-10707 Eq. 2).

Parameters:

• k (float or ndarray) –

  Longitudinal wave number [1/m]

• a (float) –

  Pipe radius [m]

• sigma0 (float) –

  DC conductivity [S/m]

• ctau (float) –

  Relaxation distance c·τ [m]

Returns:

• Zk ( complex or np.ndarray of complex ) –

  Longitudinal impedance [Ohm/m]

Source code in beamphysics/wakefields/resistive_wall/base.py

def longitudinal_impedance_round(
    k: float | np.ndarray,
    a: float,
    sigma0: float,
    ctau: float,
) -> complex | np.ndarray:
    """
    Compute longitudinal impedance Z(k) for round (circular pipe) geometry.

    Uses closed-form expression (SLAC-PUB-10707 Eq. 2).

    Parameters
    ----------
    k : float or np.ndarray
        Longitudinal wave number [1/m]
    a : float
        Pipe radius [m]
    sigma0 : float
        DC conductivity [S/m]
    ctau : float
        Relaxation distance c·τ [m]

    Returns
    -------
    Zk : complex or np.ndarray of complex
        Longitudinal impedance [Ohm/m]
    """
    k = np.asarray(k)
    scalar_input = k.ndim == 0
    k = np.atleast_1d(k)

    prefactor = Z0 / (2 * np.pi * a)
    zeta = surface_impedance(k, sigma0, ctau)

    with np.errstate(divide="ignore", invalid="ignore"):
        Zk = prefactor / (1.0 / zeta - 1j * k * a / 2)
        Zk = np.where(k == 0, 0.0 + 0.0j, Zk)

    if scalar_input:
        return complex(Zk[0])
    return Zk

longitudinal_impedance_flat

longitudinal_impedance_flat(k: float | ndarray, a: float, sigma0: float, ctau: float) -> complex | ndarray

Compute longitudinal impedance Z(k) for flat (parallel plate) geometry.

Uses numerical integration over transverse modes (SLAC-PUB-10707 Eq. 52).

Parameters:

• k (float or ndarray) –

  Longitudinal wave number [1/m]

• a (float) –

  Half-gap height between parallel plates [m]

• sigma0 (float) –

  DC conductivity [S/m]

• ctau (float) –

  Relaxation distance c·τ [m]

Returns:

• Zk ( complex or np.ndarray of complex ) –

  Longitudinal impedance [Ohm/m]

Source code in beamphysics/wakefields/resistive_wall/base.py

def longitudinal_impedance_flat(
    k: float | np.ndarray,
    a: float,
    sigma0: float,
    ctau: float,
) -> complex | np.ndarray:
    """
    Compute longitudinal impedance Z(k) for flat (parallel plate) geometry.

    Uses numerical integration over transverse modes (SLAC-PUB-10707 Eq. 52).

    Parameters
    ----------
    k : float or np.ndarray
        Longitudinal wave number [1/m]
    a : float
        Half-gap height between parallel plates [m]
    sigma0 : float
        DC conductivity [S/m]
    ctau : float
        Relaxation distance c·τ [m]

    Returns
    -------
    Zk : complex or np.ndarray of complex
        Longitudinal impedance [Ohm/m]
    """

    @np.vectorize
    def _Zk_scalar(k_val):
        if k_val == 0:
            return 0.0 + 0.0j

        def integrand(x):
            return _impedance_integrand_flat(k_val, x, a, sigma0, ctau)

        return quad_vec(integrand, 0, np.inf)[0]

    return _Zk_scalar(k)

wakefield_from_impedance

wakefield_from_impedance(z: float | ndarray, Zk_func: callable, k_max: float = 10000000.0, epsabs: float = 1e-09, epsrel: float = 1e-06) -> float | ndarray

Compute wakefield W(z) from Re[Z(k)] using a cosine transform (quadrature).

For z ≤ 0 (trailing particles):

\[W(z) = \frac{2c}{\pi} \int_0^{k_{\max}} \text{Re}[Z(k)] \cos(kz) \, dk\]

Returns 0 for z > 0 (causality).

Parameters:

• z (float or ndarray) –

  Longitudinal position [m]. Negative z is behind the source.

• Zk_func (callable) –

  Function returning complex impedance Z(k) [Ohm/m]

• k_max (float, default: 10000000.0 ) –

  Upper limit of k integration [1/m]. Default is 1e7.

• epsabs (float, default: 1e-09 ) –

  Absolute tolerance. Default is 1e-9.

• epsrel (float, default: 1e-06 ) –

  Relative tolerance. Default is 1e-6.

Returns:

• Wz ( float or ndarray ) –

  Wakefield [V/C/m]

Source code in beamphysics/wakefields/resistive_wall/base.py

def wakefield_from_impedance(
    z: float | np.ndarray,
    Zk_func: callable,
    k_max: float = 1e7,
    epsabs: float = 1e-9,
    epsrel: float = 1e-6,
) -> float | np.ndarray:
    """
    Compute wakefield W(z) from Re[Z(k)] using a cosine transform (quadrature).

    For z ≤ 0 (trailing particles):

    $$W(z) = \\frac{2c}{\\pi} \\int_0^{k_{\\max}} \\text{Re}[Z(k)] \\cos(kz) \\, dk$$

    Returns 0 for z > 0 (causality).

    Parameters
    ----------
    z : float or np.ndarray
        Longitudinal position [m]. Negative z is behind the source.
    Zk_func : callable
        Function returning complex impedance Z(k) [Ohm/m]
    k_max : float, optional
        Upper limit of k integration [1/m]. Default is 1e7.
    epsabs : float, optional
        Absolute tolerance. Default is 1e-9.
    epsrel : float, optional
        Relative tolerance. Default is 1e-6.

    Returns
    -------
    Wz : float or np.ndarray
        Wakefield [V/C/m]
    """

    def _wakefield_scalar(z_val):
        if z_val > 0:
            return 0.0

        def integrand(k):
            return np.real(Zk_func(k)) * np.cos(k * z_val)

        result, _ = quad(integrand, 0, k_max, epsabs=epsabs, epsrel=epsrel, limit=200)
        return (2 * c_light / np.pi) * result

    z = np.asarray(z)
    if z.ndim == 0:
        return _wakefield_scalar(float(z))

    return np.array([_wakefield_scalar(zi) for zi in z])

wakefield_from_impedance_fft

wakefield_from_impedance_fft(z: ndarray, Zk_func: callable, k_max: float = 10000000.0, n_fft: int = 8192) -> ndarray

Compute wakefield W(z) from Z(k) using FFT and interpolation.

Much faster than wakefield_from_impedance for arrays.

Parameters:

• z (ndarray) –

  Longitudinal positions [m]. Negative z is behind the source.

• Zk_func (callable) –

  Function returning complex impedance Z(k) [Ohm/m]

• k_max (float, default: 10000000.0 ) –

  Upper limit of k integration [1/m]. Default is 1e7.

• n_fft (int, default: 8192 ) –

  Number of FFT points. Default is 8192.

Returns:

• Wz ( ndarray ) –

  Wakefield [V/C/m]

Source code in beamphysics/wakefields/resistive_wall/base.py

def wakefield_from_impedance_fft(
    z: np.ndarray,
    Zk_func: callable,
    k_max: float = 1e7,
    n_fft: int = 8192,
) -> np.ndarray:
    """
    Compute wakefield W(z) from Z(k) using FFT and interpolation.

    Much faster than `wakefield_from_impedance` for arrays.

    Parameters
    ----------
    z : np.ndarray
        Longitudinal positions [m]. Negative z is behind the source.
    Zk_func : callable
        Function returning complex impedance Z(k) [Ohm/m]
    k_max : float, optional
        Upper limit of k integration [1/m]. Default is 1e7.
    n_fft : int, optional
        Number of FFT points. Default is 8192.

    Returns
    -------
    Wz : np.ndarray
        Wakefield [V/C/m]
    """
    from scipy.fft import irfft
    from scipy.interpolate import interp1d

    z = np.asarray(z)

    dk = k_max / (n_fft - 1)
    k_grid = np.linspace(0, k_max, n_fft)

    Zk_grid = Zk_func(k_grid)
    ReZ = np.real(Zk_grid)

    n_full = 2 * (n_fft - 1)
    dz = 2 * np.pi / (n_full * dk)
    z_grid = np.arange(n_full) * dz

    W_grid = irfft(ReZ, n=n_full) * n_full * dk * (c_light / np.pi)

    interp = interp1d(z_grid, W_grid, kind="cubic", bounds_error=False, fill_value=0.0)

    result = interp(-z)
    result = np.where(z > 0, 0.0, result)

    return result

ac_conductivity

ac_conductivity(k: float | ndarray, sigma0: float, ctau: float) -> complex | ndarray

Frequency-dependent AC conductivity with relaxation time (Drude model).

\[\sigma(k) = \frac{\sigma_0}{1 - i k c \tau}\]

Parameters:

• k (float or ndarray) –

  Longitudinal wave number [1/m]

• sigma0 (float) –

  DC conductivity [S/m]

• ctau (float) –

  Relaxation distance c·τ [m]

Returns:

• sigma ( complex or np.ndarray of complex ) –

  AC conductivity [S/m]

Source code in beamphysics/wakefields/resistive_wall/base.py

def ac_conductivity(
    k: float | np.ndarray,
    sigma0: float,
    ctau: float,
) -> complex | np.ndarray:
    """
    Frequency-dependent AC conductivity with relaxation time (Drude model).

    $$\\sigma(k) = \\frac{\\sigma_0}{1 - i k c \\tau}$$

    Parameters
    ----------
    k : float or np.ndarray
        Longitudinal wave number [1/m]
    sigma0 : float
        DC conductivity [S/m]
    ctau : float
        Relaxation distance c·τ [m]

    Returns
    -------
    sigma : complex or np.ndarray of complex
        AC conductivity [S/m]
    """
    return sigma0 / (1 - 1j * k * ctau)

surface_impedance

surface_impedance(k: float | ndarray, sigma0: float, ctau: float) -> complex | ndarray

Surface impedance for a conducting wall with AC conductivity.

\[\zeta(k) = (1 - i) \sqrt{\frac{k c}{2 \sigma(k) Z_0 c}}\]

Parameters:

• k (float or ndarray) –

  Longitudinal wave number [1/m]

• sigma0 (float) –

  DC conductivity [S/m]

• ctau (float) –

  Relaxation distance c·τ [m]

Returns:

• zeta ( complex or np.ndarray of complex ) –

  Surface impedance [dimensionless]

Source code in beamphysics/wakefields/resistive_wall/base.py

def surface_impedance(
    k: float | np.ndarray,
    sigma0: float,
    ctau: float,
) -> complex | np.ndarray:
    """
    Surface impedance for a conducting wall with AC conductivity.

    $$\\zeta(k) = (1 - i) \\sqrt{\\frac{k c}{2 \\sigma(k) Z_0 c}}$$

    Parameters
    ----------
    k : float or np.ndarray
        Longitudinal wave number [1/m]
    sigma0 : float
        DC conductivity [S/m]
    ctau : float
        Relaxation distance c·τ [m]

    Returns
    -------
    zeta : complex or np.ndarray of complex
        Surface impedance [dimensionless]
    """
    sigma = ac_conductivity(k, sigma0, ctau)
    return (1 - 1j) * np.sqrt(k * c_light / (2 * sigma * Z0 * c_light))

characteristic_length

characteristic_length(a: float, sigma0: float) -> float

Characteristic length scale s₀ for resistive wall wakefield.

From SLAC-PUB-10707 Eq. (5):

\[s_0 = \left( \frac{2 a^2}{Z_0 \sigma_0} \right)^{1/3}\]

Parameters:

• a (float) –

  Half-gap height (flat) or radius (round) [m]

• sigma0 (float) –

  DC conductivity [S/m]

Returns:

• s0 ( float ) –

  Characteristic length [m]

Source code in beamphysics/wakefields/resistive_wall/base.py

def characteristic_length(a: float, sigma0: float) -> float:
    """
    Characteristic length scale s₀ for resistive wall wakefield.

    From SLAC-PUB-10707 Eq. (5):

    $$s_0 = \\left( \\frac{2 a^2}{Z_0 \\sigma_0} \\right)^{1/3}$$

    Parameters
    ----------
    a : float
        Half-gap height (flat) or radius (round) [m]
    sigma0 : float
        DC conductivity [S/m]

    Returns
    -------
    s0 : float
        Characteristic length [m]
    """
    return (2 * a**2 / (Z0 * sigma0)) ** (1 / 3)