# ImpactX

ImpactX is a high-performance beam dynamics simulation code for modeling particle accelerators with collective effects. It is the next generation of the IMPACT-Z code, designed to run efficiently on modern GPUs and CPUs. The code uses the reference trajectory path length `s` as the independent variable of motion and includes effects from externally applied fields (magnets, accelerating cavities) as well as self-fields (space charge, CSR, wakefields). All particle tracking models are symplectic, and space charge is computed by solving the Poisson equation in the beam rest frame.

ImpactX provides both Python and C++ interfaces, supporting particle tracking, envelope (covariance matrix) tracking, and reference orbit tracking modes. The code is part of the Beam, Plasma & Accelerator Simulation Toolkit (BLAST) and features extensive support for accelerator lattice elements, beam distributions, and collective effects including 2D/2.5D/3D space charge, coherent synchrotron radiation (CSR), and incoherent synchrotron radiation (ISR). It outputs data in openPMD format for seamless integration with analysis workflows.

## ImpactX Simulation Class

The central simulation class that manages the entire beam dynamics simulation including particles, lattice, and physics models.

```python
from impactx import ImpactX, distribution, elements

# Initialize the simulation
sim = ImpactX()

# Configure simulation parameters
sim.space_charge = False  # Disable space charge (or "2D", "2p5D", "3D", "Gauss3D", "Gauss2p5D")
sim.slice_step_diagnostics = True  # Enable per-slice diagnostics
sim.particle_shape = 2  # B-spline order: 1 (linear), 2 (quadratic), 3 (cubic)

# Initialize domain decomposition and mesh
sim.init_grids()

# Set up beam reference particle
ref = sim.beam.ref
ref.set_species("electron").set_kin_energy_MeV(2.0e3)  # 2 GeV electron

# Create particle distribution and add particles
distr = distribution.Waterbag(
    lambdaX=3.9984884770e-5,
    lambdaY=3.9984884770e-5,
    lambdaT=1.0e-3,
    lambdaPx=2.6623538760e-5,
    lambdaPy=2.6623538760e-5,
    lambdaPt=2.0e-3,
    muxpx=-0.846574929020762,
    muypy=0.846574929020762,
    mutpt=0.0,
)
sim.add_particles(bunch_charge_C=1.0e-9, distr=distr, npart=10000)

# Define lattice elements
monitor = elements.BeamMonitor("monitor", backend="h5")
sim.lattice.extend([
    monitor,
    elements.Drift(name="drift1", ds=0.25, nslice=25),
    elements.Quad(name="quad1", ds=1.0, k=1.0, nslice=25),
    elements.Drift(name="drift2", ds=0.5, nslice=25),
    elements.Quad(name="quad2", ds=1.0, k=-1.0, nslice=25),
    elements.Drift(name="drift3", ds=0.25, nslice=25),
    monitor,
])

# Run particle tracking simulation
sim.track_particles()

# Clean shutdown
sim.finalize()
```

## Envelope Tracking Mode

Track the 6x6 beam covariance matrix through linearized transport maps instead of individual particles, providing fast beam envelope evolution.

```python
from impactx import ImpactX, distribution, elements

sim = ImpactX()
sim.space_charge = False
sim.slice_step_diagnostics = True
sim.init_grids()

# Reference particle
ref = sim.beam.ref
ref.set_species("electron").set_kin_energy_MeV(2.0e3)

# Initialize envelope from distribution (no particles needed)
distr = distribution.Waterbag(
    lambdaX=3.9984884770e-5,
    lambdaY=3.9984884770e-5,
    lambdaT=1.0e-3,
    lambdaPx=2.6623538760e-5,
    lambdaPy=2.6623538760e-5,
    lambdaPt=2.0e-3,
    muxpx=-0.846574929020762,
    muypy=0.846574929020762,
    mutpt=0.0,
)
sim.init_envelope(ref, distr)

# Define lattice
ns = 25
sim.lattice.extend([
    elements.Drift(name="drift1", ds=0.25, nslice=ns),
    elements.Quad(name="quad1", ds=1.0, k=1.0, nslice=ns),
    elements.Drift(name="drift2", ds=0.5, nslice=ns),
    elements.Quad(name="quad2", ds=1.0, k=-1.0, nslice=ns),
    elements.Drift(name="drift3", ds=0.25, nslice=ns),
])

# Run envelope tracking
sim.track_envelope()
sim.finalize()
```

## Space Charge Simulation

Enable self-consistent space charge calculations with configurable models and mesh parameters.

```python
from impactx import ImpactX, distribution, elements

sim = ImpactX()

# Space charge configuration
sim.max_level = 0  # No mesh refinement
sim.n_cell = [32, 32, 1]  # Grid cells (x, y, z)
sim.blocking_factor_x = [16]
sim.blocking_factor_y = [16]
sim.blocking_factor_z = [1]

sim.particle_shape = 2  # Quadratic B-spline
sim.space_charge = "2p5D"  # Options: "2D", "2p5D", "3D", "Gauss3D", "Gauss2p5D"
sim.poisson_solver = "fft"  # Options: "fft", "multigrid"
sim.dynamic_size = True
sim.prob_relative = [1.1]  # Mesh padding factor

# 2.5D specific parameters
sim.space_charge_num_longitudinal_bins = 100
sim.space_charge_apply_longitudinal_kick = True

sim.slice_step_diagnostics = True
sim.init_grids()

# Reference particle and distribution
ref = sim.beam.ref
ref.set_species("electron").set_kin_energy_MeV(100)

distr = distribution.Gaussian(
    lambdaX=3.9984884770e-5,
    lambdaY=3.9984884770e-5,
    lambdaT=1.0e-3,
    lambdaPx=2.6623538760e-5,
    lambdaPy=2.6623538760e-5,
    lambdaPt=2.0e-3,
    muxpx=-0.846574929020762,
    muypy=0.846574929020762,
    mutpt=0.0,
)
sim.add_particles(bunch_charge_C=1.0e-9, distr=distr, npart=10000)

# Use chromatic elements for 2.5D space charge
ns = 25
sim.lattice.extend([
    elements.ChrDrift(name="drift1", ds=0.25, nslice=ns),
    elements.ChrQuad(name="quad1", ds=1.0, k=1.0, nslice=ns),
    elements.ChrDrift(name="drift2", ds=0.5, nslice=ns),
    elements.ChrQuad(name="quad2", ds=1.0, k=-1.0, nslice=ns),
    elements.ChrDrift(name="drift3", ds=0.25, nslice=ns),
])

sim.track_particles()
sim.finalize()
```

## Particle Distributions

Create initial beam phase space distributions with various models including Gaussian, Waterbag, KV, and thermal distributions.

```python
from impactx import distribution

# Gaussian distribution with optional truncation
gauss = distribution.Gaussian(
    lambdaX=1.0e-4,      # Phase space ellipse x-intercept (m)
    lambdaY=1.0e-4,      # Phase space ellipse y-intercept (m)
    lambdaT=1.0e-3,      # Phase space ellipse t-intercept (m)
    lambdaPx=1.0e-5,     # Phase space ellipse px-intercept (rad)
    lambdaPy=1.0e-5,     # Phase space ellipse py-intercept (rad)
    lambdaPt=1.0e-3,     # Phase space ellipse pt-intercept (rad)
    muxpx=0.0,           # X-Px correlation
    muypy=0.0,           # Y-Py correlation
    mutpt=0.0,           # T-Pt correlation
    meanX=0.0,           # Centroid offset in x (m)
    meanPx=0.0,          # Centroid offset in px (rad)
    cutX=3.0,            # Truncation at 3 sigma in (x,px) plane
    cutY=3.0,            # Truncation at 3 sigma in (y,py) plane
    cutT=0.0,            # No truncation in (t,pt) plane (0 = no cut)
)

# 6D Waterbag distribution (uniform in phase space)
waterbag = distribution.Waterbag(
    lambdaX=1.0e-4, lambdaY=1.0e-4, lambdaT=1.0e-3,
    lambdaPx=1.0e-5, lambdaPy=1.0e-5, lambdaPt=1.0e-3,
)

# K-V distribution (4D transverse uniform + longitudinal Gaussian)
kv = distribution.KVdist(
    lambdaX=1.0e-4, lambdaY=1.0e-4, lambdaT=1.0e-3,
    lambdaPx=1.0e-5, lambdaPy=1.0e-5, lambdaPt=1.0e-3,
)

# 6D Kurth distribution
kurth6d = distribution.Kurth6D(
    lambdaX=1.0e-4, lambdaY=1.0e-4, lambdaT=1.0e-3,
    lambdaPx=1.0e-5, lambdaPy=1.0e-5, lambdaPt=1.0e-3,
)

# Thermal distribution for space charge benchmarks
thermal = distribution.Thermal(
    k=100.0,           # External focusing strength (1/m)
    kT=1.0e-6,         # Temperature of core population
    kT_halo=2.0e-6,    # Temperature of halo population
    normalize=1.0,     # Normalizing constant for core
    normalize_halo=0.1,# Normalizing constant for halo
    halo=0.1,          # Fraction of charge in halo
)
```

## Quadrupole Magnet

Standard quadrupole magnet for transverse focusing with optional chromatic and exact nonlinear models.

```python
from impactx import elements

# Standard linear quadrupole
quad = elements.Quad(
    name="qf1",
    ds=0.5,           # Length (m)
    k=2.5,            # Strength (1/m^2): k>0 horizontal focusing
    dx=0.0,           # Horizontal misalignment error (m)
    dy=0.0,           # Vertical misalignment error (m)
    rotation=0.0,     # Rotation error (degrees)
    aperture_x=0.05,  # Horizontal half-aperture (m)
    aperture_y=0.05,  # Vertical half-aperture (m)
    nslice=10,        # Slices for space charge
)

# Chromatic quadrupole (exact pt dependence)
chrquad = elements.ChrQuad(
    name="qf_chr",
    ds=0.5,
    k=2.5,
    unit=0,           # 0: k in 1/m^2 (MAD-X), 1: k in T/m
    nslice=10,
)

# Exact nonlinear quadrupole (full kinematic nonlinearities)
exact_quad = elements.ExactQuad(
    name="qf_exact",
    ds=0.5,
    k=2.5,
    unit=0,
    int_order=4,      # Symplectic integration order: 2, 4, or 6
    mapsteps=5,       # Integration steps per slice
    nslice=10,
)
```

## Dipole Bending Magnet

Sector bend magnets for beam deflection with various physics models.

```python
from impactx import elements

# Ideal sector bend
sbend = elements.Sbend(
    name="b1",
    ds=1.0,           # Arc length (m)
    rc=5.0,           # Radius of curvature (m), positive = bend right
    aperture_x=0.05,
    aperture_y=0.02,
    nslice=20,
)

# Exact sector bend (nonlinear map)
exact_sbend = elements.ExactSbend(
    name="b1_exact",
    ds=1.0,           # Arc length (m)
    phi=11.46,        # Bend angle (degrees)
    B=0.0,            # Field (T), 0 = derive from ds/phi
    nslice=20,
)

# Combined function bend (dipole + quadrupole)
cfbend = elements.CFbend(
    name="cf1",
    ds=1.0,
    rc=5.0,
    k=0.5,            # Quadrupole strength (1/m^2)
    nslice=20,
)

# Dipole edge focusing
dipedge = elements.DipEdge(
    name="edge1",
    psi=0.1,          # Pole face angle (radians)
    rc=5.0,           # Bend radius (m)
    g=0.05,           # Gap parameter (m)
    K2=0.5,           # FINT parameter (MAD-X convention)
    model="linear",   # "linear" or "nonlinear"
    location="entry", # "entry" or "exit"
)
```

## RF Cavity

Radiofrequency cavity for beam acceleration using on-axis field data or Fourier coefficients.

```python
import numpy as np
from impactx import elements

# Load on-axis field data
data = np.loadtxt("onaxis_data.in")
z = data[:, 0]
ez_onaxis = data[:, 1]

# RF cavity from field data
rf = elements.RFCavity(
    name="cavity1",
    ds=1.31879807,    # Length (m)
    escale=62.0,      # E-field scale: (peak Ez in MV/m) / (rest mass in MeV)
    z=z,              # Longitudinal positions (m)
    field_on_axis=ez_onaxis,  # Normalized Ez values
    ncoef=25,         # Number of Fourier coefficients to compute
    freq=1.3e9,       # RF frequency (Hz)
    phase=85.5,       # RF phase (degrees)
    mapsteps=100,     # Integration steps per slice
    nslice=4,
)

# RF cavity with pre-computed Fourier coefficients
rf_fourier = elements.RFCavity(
    name="cavity2",
    ds=1.0,
    escale=50.0,
    cos_coefficients=[1.0, 0.1, 0.01],  # Cosine Fourier coefficients
    sin_coefficients=[0.0, 0.05, 0.005], # Sine Fourier coefficients
    freq=1.3e9,
    phase=0.0,        # On-crest acceleration
    mapsteps=50,
    nslice=4,
)

# Short RF cavity (thin lens approximation)
short_rf = elements.ShortRF(
    name="buncher",
    V=0.01,           # Normalized voltage: (max energy gain) / (rest mass)
    freq=1.3e9,
    phase=-90.0,      # -90 deg = bunching, 0 = acceleration
)
```

## Solenoid Magnet

Solenoid focusing magnets with hard-edge and soft-edge models.

```python
from impactx import elements

# Ideal hard-edge solenoid
sol = elements.Sol(
    name="sol1",
    ds=0.5,           # Length (m)
    ks=1.5,           # Strength (1/m): Bz(T) / rigidity(T-m)
    aperture_x=0.02,
    aperture_y=0.02,
    nslice=10,
)

# Soft-edge solenoid with realistic fringe fields
soft_sol = elements.SoftSolenoid(
    name="sol_soft",
    ds=0.5,
    bscale=3.0,       # Field scale: Bz(T) / rigidity(T-m) (if unit=0)
    unit=0,           # 0: normalized, 1: SI units (Tesla)
    # cos/sin_coefficients: optional, defaults to thin-shell model
    mapsteps=10,
    nslice=10,
)
```

## Beam Monitor

Write particle data to openPMD files for diagnostics and restart capabilities.

```python
from impactx import elements

# Basic beam monitor
monitor = elements.BeamMonitor(
    name="diag1",
    backend="h5",     # Options: "h5", "bp4", "bp5", "json", "default"
    encoding="g",     # (g)roup, (f)ile, or (v)ariable based
    period_sample_intervals=1,  # Output every N periods (for rings)
)

# Monitor with IOTA nonlinear lens invariants
monitor_nll = elements.BeamMonitor(name="diag_nll", backend="h5")
monitor_nll.nonlinear_lens_invariants = True
monitor_nll.alpha = -1.0  # Twiss alpha at monitor location
monitor_nll.beta = 2.5    # Twiss beta (m) at monitor location
monitor_nll.tn = 0.4      # IOTA NLL strength
monitor_nll.cn = 0.01     # IOTA NLL scale (m^0.5)
```

## Lattice Operations

Build and manipulate accelerator lattices with filtering, selection, and transformation capabilities.

```python
from impactx import elements

# Create and populate lattice
lattice = elements.KnownElementsList()
lattice.extend([
    elements.Drift(name="d1", ds=1.0),
    elements.Quad(name="qf", ds=0.5, k=2.0),
    elements.Drift(name="d2", ds=2.0),
    elements.Quad(name="qd", ds=0.5, k=-2.0),
    elements.Drift(name="d3", ds=1.0),
])

# Filter elements by type or name
drifts = lattice.select(kind="Drift")           # All Drift elements
quads = lattice.select(kind=elements.Quad)      # All Quad elements
all_quads = lattice.select(kind=r".*Quad")      # Regex: Quad, ChrQuad, ExactQuad
named = lattice.select(name="qf")               # By name
combined = lattice.select(kind="Drift").select(name="d1")  # Chained AND filter

# Modify selected elements
drifts[0].ds = 1.5  # Modifies original lattice

# Replace elements
lattice.select(kind=r".*Quad").replace_with_drifts()  # Replace quads with drifts
lattice.select(kind="Quad").replace_each(
    elements.Quad(name="new", ds=0.1, k=1.5),
    keep_name=True,
    keep_ds=False,
)

# Delete selected elements
lattice.select(kind="Drift").select(name="d2").delete()

# Compute transfer map
ref = sim.beam.ref
M = lattice.transfer_map(ref, order="linear")  # 6x6 transfer matrix

# Load from MAD-X file
lattice.load_file("lattice.madx", nslice=10)

# Serialize/deserialize
import json
with open("lattice.json", "w") as f:
    json.dump(lattice.to_dicts(), f, indent=2)

lattice2 = elements.KnownElementsList()
with open("lattice.json") as f:
    lattice2.from_dicts(json.load(f))
```

## Ring Simulation with Multiple Periods

Simulate circular accelerators by repeating the lattice for multiple turns.

```python
from impactx import ImpactX, distribution, elements

sim = ImpactX()
sim.space_charge = False
sim.slice_step_diagnostics = True
sim.init_grids()

# Reference particle (proton beam)
ref = sim.beam.ref
ref.set_species("proton").set_kin_energy_MeV(2.5)

# Distribution
distr = distribution.Waterbag(
    lambdaX=1.588960728035e-3,
    lambdaY=2.496625268437e-3,
    lambdaT=1.0e-3,
    lambdaPx=2.8320397837724e-3,
    lambdaPy=1.802433091137e-3,
    lambdaPt=0.0,
)
sim.add_particles(bunch_charge_C=1.0e-9, distr=distr, npart=10000)

# Build ring lattice (simplified IOTA-like)
ns = 10
rc = 0.82223
sbend = elements.Sbend(name="sbend", ds=0.43, rc=rc)
edge = elements.DipEdge(name="edge", psi=0.0, rc=rc, g=0.058, K2=0.5)
drift = elements.Drift(name="dr", ds=0.5, nslice=ns)
quad_f = elements.Quad(name="qf", ds=0.21, k=8.0, nslice=ns)
quad_d = elements.Quad(name="qd", ds=0.21, k=-8.0, nslice=ns)

monitor = elements.BeamMonitor("monitor", backend="h5")
sim.lattice.extend([
    monitor, drift, quad_f, drift, edge, sbend, edge, drift,
    quad_d, drift, edge, sbend, edge, drift, quad_f, drift, monitor,
])

# Number of turns to simulate
sim.periods = 100

# User-defined hook for turn-by-turn analysis
def before_period(sim):
    turn = sim.tracking_period
    beam = sim.beam
    print(f"Turn {turn}: {beam.ref.s:.3f} m")

sim.hook["before_period"] = before_period

sim.track_particles()
sim.finalize()
```

## Input File Format

Configure simulations using text-based input files with element definitions and parameters.

```
###############################################################################
# Particle Beam
###############################################################################
beam.npart = 10000
beam.units = static
beam.kin_energy = 2.0e3
beam.charge = 1.0e-9
beam.particle = electron
beam.distribution = waterbag
beam.lambdaX = 3.9984884770e-5
beam.lambdaY = beam.lambdaX
beam.lambdaT = 1.0e-3
beam.lambdaPx = 2.6623538760e-5
beam.lambdaPy = beam.lambdaPx
beam.lambdaPt = 2.0e-3
beam.muxpx = -0.846574929020762
beam.muypy = -beam.muxpx
beam.mutpt = 0.0

###############################################################################
# Beamline
###############################################################################
lattice.elements = monitor drift1 monitor quad1 monitor drift2 monitor quad2 monitor drift3 monitor
lattice.nslice = 25
lattice.periods = 1

monitor.type = beam_monitor
monitor.backend = h5

drift1.type = drift
drift1.ds = 0.25

quad1.type = quad
quad1.ds = 1.0
quad1.k = 1.0

drift2.type = drift
drift2.ds = 0.5

quad2.type = quad
quad2.ds = 1.0
quad2.k = -1.0

drift3.type = drift
drift3.ds = 0.25

###############################################################################
# Algorithms
###############################################################################
algo.space_charge = false

###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = true
```

Run with: `impactx input_fodo.in` or override parameters: `impactx input_fodo.in quad1.k=1.5`

## Apertures and Collimators

Define particle apertures with various shapes for beam scraping and loss modeling.

```python
from impactx import elements

# Rectangular aperture
rect_aperture = elements.Aperture(
    aperture_x=0.02,  # Horizontal half-aperture (m)
    aperture_y=0.01,  # Vertical half-aperture (m)
    repeat_x=0.0,     # Horizontal period for repeated masking (0=off)
    repeat_y=0.0,     # Vertical period for repeated masking (0=off)
    shift_odd_x=False,# Hexagonal pattern shift
    shape="rectangular",  # or "elliptical"
)

# Elliptical aperture
ellip_aperture = elements.Aperture(
    aperture_x=0.02,
    aperture_y=0.01,
    repeat_x=0.0,
    repeat_y=0.0,
    shift_odd_x=False,
    shape="elliptical",
)

# Polygon aperture (counter-clockwise, closed curve)
import numpy as np
vertices_x = [0.01, -0.01, -0.01, 0.01, 0.01]  # First = last
vertices_y = [0.02, 0.02, -0.02, -0.02, 0.02]
poly_aperture = elements.PolygonAperture(
    vertices_x=vertices_x,
    vertices_y=vertices_y,
    min_radius2=1e-4,  # Skip polygon calc for r^2 < this (m^2)
    repeat_x=0.0,
    repeat_y=0.0,
    shift_odd_x=False,
    action="transmit",  # "transmit" or "absorb"
)
```

## Programmable Element

Create custom beam optics elements with user-defined push functions for advanced physics models.

```python
from impactx import ImpactX, distribution, elements
import numpy as np

sim = ImpactX()
sim.space_charge = False
sim.init_grids()

ref = sim.beam.ref
ref.set_species("electron").set_kin_energy_MeV(100.0)

distr = distribution.Gaussian(
    lambdaX=1e-4, lambdaY=1e-4, lambdaT=1e-3,
    lambdaPx=1e-5, lambdaPy=1e-5, lambdaPt=1e-3,
)
sim.add_particles(1e-9, distr, 1000)

# Custom thin lens using programmable element
def custom_push(pc, step):
    """Apply thin quadrupole kick to all particles."""
    ref = pc.ref
    k = 1.0  # Quadrupole strength

    for pti in pc:
        soa = pti.soa()
        x = soa.real_data("x")
        y = soa.real_data("y")
        px = soa.real_data("px")
        py = soa.real_data("py")

        # Apply thin quad kick
        for i in range(len(x)):
            px[i] -= k * x[i]
            py[i] += k * y[i]

prog = elements.Programmable(name="custom_quad", ds=0.0, nslice=1)
prog.push = custom_push

sim.lattice.extend([
    elements.Drift(ds=1.0),
    prog,
    elements.Drift(ds=1.0),
])

sim.track_particles()
sim.finalize()
```

## Beam Moments and Analysis

Access beam statistical moments during and after simulation for analysis.

```python
from impactx import ImpactX, distribution, elements

sim = ImpactX()
sim.space_charge = False
sim.init_grids()

ref = sim.beam.ref
ref.set_species("electron").set_kin_energy_MeV(100.0)

distr = distribution.Gaussian(
    lambdaX=1e-4, lambdaY=1e-4, lambdaT=1e-3,
    lambdaPx=1e-5, lambdaPy=1e-5, lambdaPt=1e-3,
)
sim.add_particles(1e-9, distr, 10000)

# Enable beam moments history storage
sim.beam.store_beam_moments = True

sim.lattice.extend([
    elements.Drift(ds=1.0, nslice=10),
    elements.Quad(ds=0.5, k=2.0, nslice=10),
    elements.Drift(ds=1.0, nslice=10),
])

sim.track_particles()

# Get current beam moments
moments = sim.beam.beam_moments()
print(f"sig_x: {moments['sig_x']:.6e} m")
print(f"sig_y: {moments['sig_y']:.6e} m")
print(f"emittance_x: {moments['emittance_x']:.6e} m-rad")
print(f"emittance_y: {moments['emittance_y']:.6e} m-rad")

# Get full history as pandas DataFrame
history = sim.beam.beam_moments_history()
print(history[['s', 'sig_x', 'sig_y', 'emittance_x', 'emittance_y']])

# Position statistics
x_min, y_min, z_min, x_max, y_max, z_max = sim.beam.min_and_max_positions()
x_mean, x_std, y_mean, y_std, z_mean, z_std = sim.beam.mean_and_std_positions()

sim.finalize()
```

ImpactX is designed for high-performance simulations of particle accelerators, from single-pass linacs to multi-turn storage rings. The code supports both CPU and GPU execution, making it suitable for everything from rapid design studies on workstations to large-scale production runs on supercomputers. The Python interface enables seamless integration with optimization frameworks, machine learning workflows, and data analysis pipelines.

Common use cases include FODO channel matching, chicane and bunch compressor design with CSR effects, linac acceleration with space charge, storage ring dynamics with multiple turns, and envelope matching for injector design. The code's modular architecture allows users to combine standard beamline elements with custom programmable elements for specialized physics models, while the openPMD output format ensures compatibility with standard post-processing tools.