Soft-edge solenoid

Proton beam propagating through a 6 m region containing a soft-edge solenoid.

The solenoid model used is the default thin-shell model described in: P. Granum et al, “Efficient calculations of magnetic fields of solenoids for simulations,” NIMA 1034, 166706 (2022) DOI:10.1016/j.nima.2022.166706

The solenoid is a cylindrical current sheet with a length of 1 m and a radius of 0.1667 m, corresponding to an aspect ratio diameter/length = 1/3. The peak magnetic field on-axis is 3 T.

We use a 250 MeV proton beam with initial unnormalized rms emittance of 1 micron in the horizontal plane, and 2 micron in the vertical plane.

In this test, the initial and final values of \(\sigma_x\), \(\sigma_y\), \(\sigma_t\), \(\epsilon_x\), \(\epsilon_y\), and \(\epsilon_t\) must agree with nominal values.

Run

This example can be run either as:

  • Python script: python3 run_solenoid_softedge.py or

  • ImpactX executable using an input file: impactx input_solenoid_softedge.in

For MPI-parallel runs, prefix these lines with mpiexec -n 4 ... or srun -n 4 ..., depending on the system.

Listing 47 You can copy this file from examples/solenoid_softedge/run_solenoid_softedge.py.
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Chad Mitchell, Axel Huebl
# License: BSD-3-Clause-LBNL
#
# -*- coding: utf-8 -*-

from impactx import ImpactX, distribution, elements

sim = ImpactX()

# set numerical parameters and IO control
sim.space_charge = False
# sim.diagnostics = False  # benchmarking
sim.slice_step_diagnostics = False

# domain decomposition & space charge mesh
sim.init_grids()

# load a 250 MeV proton beam with an initial
# horizontal rms emittance of 1 um and an
# initial vertical rms emittance of 2 um
kin_energy_MeV = 250.0  # reference energy
bunch_charge_C = 1.0e-9  # used with space charge
npart = 10000  # number of macro particles

#   reference particle
ref = sim.particle_container().ref_particle()
ref.set_charge_qe(1.0).set_mass_MeV(938.27208816).set_kin_energy_MeV(kin_energy_MeV)

#   particle bunch
distr = distribution.Waterbag(
    lambdaX=1.559531175539e-3,
    lambdaY=2.205510139392e-3,
    lambdaT=1.0e-3,
    lambdaPx=6.41218345413e-4,
    lambdaPy=9.06819680526e-4,
    lambdaPt=1.0e-3,
)
sim.add_particles(bunch_charge_C, distr, npart)

# design the accelerator lattice
sol = elements.SoftSolenoid(
    name="sol1",
    ds=6.0,
    bscale=1.233482899483985,
    cos_coefficients=[
        0.350807812299706,
        0.323554693720069,
        0.260320578919415,
        0.182848575294969,
        0.106921016050403,
        4.409581845710694e-002,
        -9.416427163897508e-006,
        -2.459452716865687e-002,
        -3.272762575737291e-002,
        -2.936414401076162e-002,
        -1.995780078926890e-002,
        -9.102893342953847e-003,
        -2.456410658713271e-006,
        5.788233017324325e-003,
        8.040408292420691e-003,
        7.480064552867431e-003,
        5.230254569468851e-003,
        2.447614547094685e-003,
        -1.095525090532255e-006,
        -1.614586867387170e-003,
        -2.281365457438345e-003,
        -2.148709081338292e-003,
        -1.522541739363011e-003,
        -7.185505862719508e-004,
        -6.171194824600157e-007,
        4.842109305036943e-004,
        6.874508102002901e-004,
        6.535550288205728e-004,
        4.648795813759210e-004,
        2.216564722797528e-004,
        -4.100982995210341e-007,
        -1.499332112463395e-004,
        -2.151538438342482e-004,
        -2.044590946652016e-004,
        -1.468242784844341e-004,
    ],
    sin_coefficients=[
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
    ],
    mapsteps=200,
    nslice=4,
)

# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")

sim.lattice.extend(
    [
        monitor,
        sol,
        monitor,
    ]
)

# run simulation
sim.track_particles()

# clean shutdown
sim.finalize()
Listing 48 You can copy this file from examples/solenoid_softedge/input_solenoid_softedge.in.
###############################################################################
# Particle Beam(s)
###############################################################################
beam.npart = 10000
beam.units = static
beam.kin_energy = 250.0
beam.charge = 1.0e-9
beam.particle = proton
beam.distribution = waterbag
beam.lambdaX = 1.559531175539e-3
beam.lambdaY = 2.205510139392e-3
beam.lambdaT = 1.0e-3
beam.lambdaPx = 6.41218345413e-4
beam.lambdaPy = 9.06819680526e-4
beam.lambdaPt = 1.0e-3
beam.muxpx = 0.0
beam.muypy = 0.0
beam.mutpt = 0.0


###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.elements = monitor sol1 monitor
lattice.nslice = 1

monitor.type = beam_monitor
monitor.backend = h5

sol1.type = solenoid_softedge
sol1.ds = 6.0
sol1.bscale = 1.233482899483985
sol1.mapsteps = 800


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


###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = false

Analyze

We run the following script to analyze correctness:

Script analysis_solenoid_softedge.py
Listing 49 You can copy this file from examples/solenoid_softedge/analysis_solenoid_softedge.py.
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Chad Mitchell, Axel Huebl
# License: BSD-3-Clause-LBNL
#

import numpy as np
import openpmd_api as io
from scipy.stats import moment


def get_moments(beam):
    """Calculate standard deviations of beam position & momenta
    and emittance values

    Returns
    -------
    sigx, sigy, sigt, emittance_x, emittance_y, emittance_t
    """
    sigx = moment(beam["position_x"], moment=2) ** 0.5  # variance -> std dev.
    sigpx = moment(beam["momentum_x"], moment=2) ** 0.5
    sigy = moment(beam["position_y"], moment=2) ** 0.5
    sigpy = moment(beam["momentum_y"], moment=2) ** 0.5
    sigt = moment(beam["position_t"], moment=2) ** 0.5
    sigpt = moment(beam["momentum_t"], moment=2) ** 0.5

    epstrms = beam.cov(ddof=0)
    emittance_x = (sigx**2 * sigpx**2 - epstrms["position_x"]["momentum_x"] ** 2) ** 0.5
    emittance_y = (sigy**2 * sigpy**2 - epstrms["position_y"]["momentum_y"] ** 2) ** 0.5
    emittance_t = (sigt**2 * sigpt**2 - epstrms["position_t"]["momentum_t"] ** 2) ** 0.5

    return (sigx, sigy, sigt, emittance_x, emittance_y, emittance_t)


# initial/final beam
series = io.Series("diags/openPMD/monitor.h5", io.Access.read_only)
last_step = list(series.iterations)[-1]
initial = series.iterations[1].particles["beam"].to_df()
final = series.iterations[last_step].particles["beam"].to_df()

# compare number of particles
num_particles = 10000
assert num_particles == len(initial)
assert num_particles == len(final)

print("Initial Beam:")
sigx, sigy, sigt, emittance_x, emittance_y, emittance_t = get_moments(initial)
print(f"  sigx={sigx:e} sigy={sigy:e} sigt={sigt:e}")
print(
    f"  emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
)

atol = 0.0  # ignored
rtol = 2.0 * num_particles**-0.5  # from random sampling of a smooth distribution
print(f"  rtol={rtol} (ignored: atol~={atol})")

assert np.allclose(
    [sigx, sigy, sigt, emittance_x, emittance_y, emittance_t],
    [
        1.559531175539e-3,
        2.205510139392e-3,
        1.0e-3,
        1.0e-6,
        2.0e-6,
        1.0e-6,
    ],
    rtol=rtol,
    atol=atol,
)


print("")
print("Final Beam:")
sigx, sigy, sigt, emittance_x, emittance_y, emittance_t = get_moments(final)
print(f"  sigx={sigx:e} sigy={sigy:e} sigt={sigt:e}")
print(
    f"  emittance_x={emittance_x:e} emittance_y={emittance_y:e} emittance_t={emittance_t:e}"
)

atol = 0.0  # ignored
rtol = 2.0 * num_particles**-0.5  # from random sampling of a smooth distribution
print(f"  rtol={rtol} (ignored: atol~={atol})")

assert np.allclose(
    [sigx, sigy, sigt, emittance_x, emittance_y, emittance_t],
    [
        2.425578906459e-3,
        2.654015302646e-3,
        9.985897906860e-3,
        1.365357890e-6,
        1.634641555e-6,
        1.000000000e-6,
    ],
    rtol=rtol,
    atol=atol,
)

Exactly-solvable (non-uniform) soft-edge solenoid

This benchmark checks the calculation of the linear map for a soft-edge solenoid with a non-uniform longitudinal on-axis magnetic field profile, in the special case of a field profile for which the map is exactly-solvable.

The test involves 250 MeV protons propagating through a 2 m region with a solenoidal magnetic field.

The user can run this test for various values of the magnetic field strength by modifying the parameter “bscale” in both the input file and in the analysis script. (The values in the two files must be consistent.)

In this test, all 36 elements of the 6x6 transport matrix must agree with predicted values to within numerical tolerance (currently 1e-7).

Details

Introduction

Consider a Hamiltonian \(H\) quadratic in the phase space variables \(\zeta=(x,p_x,y,p_y,t,p_t)\). Then \(H\) can be written in the form:

\[H(\zeta,z)=\frac{1}{2}\zeta^TS(z)\zeta,\]

where \(S\) is a \(6\times 6\) symmetric matrix. The linear symplectic map \(M\) describing transport from \(z_0\mapsto z\) is the unique solution of:

\[\frac{dM(z)}{dz}=JS(z)M(z),\quad\quad M(z_0)=I,\]

where \(J\) is the \(6\times 6\) matrix of the symplectic form.

From R. Ryne’s USPAS notes, the Hamiltonian for linear transport in a solenoid is given by:

\[H=H_2^{foc}+H_2^{rot},\]

where

\[\begin{split}H_2^{foc}&=\frac{p_t^2}{2\beta_0^2\gamma_0^2}+\frac{1}{2}(p_x^2+p_y^2)+\frac{\alpha^2}{2}(x^2+y^2), \\ H_2^{rot}&=-\alpha(xp_y-yp_x),\end{split}\]

and \(\alpha\) is related to the on-axis (solenoidal) magnetic field \(B_z\) by:

\[\alpha(z)=\frac{1}{2}\frac{B_z(z)}{B\rho}.\]

We consider an on-axis solenoid field profile of the following special form:

\[B_z(z)=\frac{B_0}{1+(z/g)^2},\]

where \(B_0\) is just the peak magnetic field on-axis, and \(g\) is a length parameter describing the rate of decay of the field away from its peak at \(z=0\).

Since \(H_2^{foc}\) and \(H_2^{rot}\) Poisson-commute, the linear map \(M\) can be written in the factorized form:

\[M=M^0R,\]

where \(M^0\) is the linear map associated with \(H_2^{foc}\) , and \(R\) is the linear map associated with \(H_2^{rot}\) (see ref).

To express the linear map between points \(z_0\) and \(z\), we define the following dimensionless parameters:

\[\Delta\phi=\tan^{-1}\left(\frac{z}{g}\right)-\tan^{-1}\left(\frac{z_0}{g}\right),\quad\quad b=\frac{1}{2}\frac{B_0g}{B\rho},\quad\quad \beta=\sqrt{1+b^2}.\]

Then we obtain for the focusing matrix:

\[\begin{split}M^0_{11}&=M^0_{33}=\sqrt{\frac{g^2+z^2}{g^2+z_0^2}}\left\{\cos(\beta\Delta\phi)-\frac{z_0\sin(\beta\Delta\phi)}{g\beta}\right\}, \\ M^0_{12}&=M^0_{34}=\sqrt{(g^2+z^2)(g^2+z_0^2)}\frac{\sin(\beta\Delta\phi)}{g\beta}, \\ M^0_{21}&=M^0_{43}=\frac{\beta(z-z_0)\cos(\beta\Delta\phi)-g(\beta^2+zz_0/g^2)\sin(\beta\Delta\phi)}{\beta\sqrt{(g^2+z^2)(g^2+z_0^2)}}, \\ M^0_{22}&=M^0_{44}=\sqrt{\frac{g^2+z_0^2}{g^2+z^2}}\left\{\cos(\beta\Delta\phi)+\frac{z\sin(\beta\Delta\phi)}{g\beta}\right\}, \\ M^0_{56}&=\frac{z-z_0}{\beta_0^2\gamma_0^2},\quad\quad M^0_{55}=M^0_{66}=1.\end{split}\]

Likewise, for the rotation matrix we obtain:

\[\begin{split}R_{11}&=R_{22}=\cos(b\Delta\phi),\quad\quad R_{13}=R_{24}=\sin(b\Delta\phi), \\ R_{31}&=R_{42}=-\sin(b\Delta\phi),\quad\quad R_{33}=R_{44}=\cos(b\Delta\phi), \\ R_{55}&=R_{66}=1.\end{split}\]

The reader can verify that the map equations and the symplectic condition are satisfied:

\[\begin{split}&\frac{dM^0(z)}{dz}=JS^{foc}(z)M^0(z),\quad\quad M^0(z_0)=I, \quad\quad (M^0)^TJM^0=J\\ &\frac{dR(z)}{dz}=JS^{rot}(z)R(z),\quad\quad R(z_0)=I, \quad\quad R^TJR=J.\end{split}\]

Here the \(6\times 6\) symmetric matrices \(S^{foc}\) and \(S^{rot}\) are defined such that:

\[H^{foc}(\zeta,z)=\frac{1}{2}\zeta^TS^{foc}(z)\zeta,\quad\quad H^{rot}(\zeta,z)=\frac{1}{2}\zeta^TS^{rot}(z)\zeta.\]

End-to-end map

A case of special interest occurs when the map is taken over a region symmetric about the solenoid midpoint at \(z=0\). In this case, we define \(z_0=-\lambda g\) and \(z=\lambda g\). Then we obtain:

\[\Delta\phi=2\tan^{-1}\lambda,\quad\quad \lambda=\frac{L}{2g},\]

where \(L\) is the length of the field region over which the map is computed. The map takes the form:

\[\begin{split}M^0_{11}&=M^0_{33}=\cos(\beta\Delta\phi)+\frac{\lambda\sin(\beta\Delta\phi)}{\beta}, \\ M^0_{12}&=M^0_{34}=g(1+\lambda^2)\frac{\sin(\beta\Delta\phi)}{\beta}, \\ M^0_{21}&=M^0_{43}=\frac{2\beta\lambda\cos(\beta\Delta\phi)+(\lambda^2-\beta^2)\sin(\beta\Delta\phi)}{\beta g(1+\lambda^2)}, \\ M^0_{22}&=M^0_{44}=\cos(\beta\Delta\phi)+\frac{\lambda\sin(\beta\Delta\phi)}{\beta}, \\ M^0_{56}&=\frac{2\lambda g}{\beta_0^2\gamma_0^2},\quad\quad M^0_{55}=M^0_{66}=1.\end{split}\]

The expression for the rotation matrix \(R\) is unchanged. Note from the form of \(R\) that the Larmor rotation angle \(\Omega\) is given by:

\[\Omega=b\Delta\phi=2b\tan^{-1}\lambda.\]

This has a well defined limit as \(\lambda\rightarrow\infty\), which is:

\[\lim_{\lambda\rightarrow\infty}\Omega= b\pi =\frac{\pi}{2}\frac{B_0g}{B\rho}.\]

This is consistent with the Larmor angle we expect, namely:

\[\Omega_{\infty}=\int_{-\infty}^{\infty}\alpha(z)dz=\frac{\pi}{2}\frac{B_0g}{B\rho}.\]

References

  • R.D. Ryne, “Computational Methods in Accelerator Physics,” USPAS Lecture Notes (2009)

Run

This example can be run as:

  • Python script: python3 run_solenoid_softedge.py

For MPI-parallel runs, prefix these lines with mpiexec -n 4 ... or srun -n 4 ..., depending on the system.

Listing 50 You can copy this file from examples/solenoid_softedge/run_solenoid_softedge_solvable.py.
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Chad Mitchell, Axel Huebl
# License: BSD-3-Clause-LBNL
#
# -*- coding: utf-8 -*-

from impactx import ImpactX, distribution, elements

sim = ImpactX()

# set numerical parameters and IO control
sim.space_charge = False
# sim.diagnostics = False  # benchmarking
sim.slice_step_diagnostics = False

# domain decomposition & space charge mesh
sim.init_grids()

# load a 250 MeV proton beam with an initial
# horizontal rms emittance of 1 um and an
# initial vertical rms emittance of 2 um
kin_energy_MeV = 250.0  # reference energy
bunch_charge_C = 1.0e-9  # used with space charge
npart = 10000  # number of macro particles

#   reference particle
ref = sim.particle_container().ref_particle()
ref.set_charge_qe(1.0).set_mass_MeV(938.27208816).set_kin_energy_MeV(kin_energy_MeV)

#   particle bunch
distr = distribution.Waterbag(
    lambdaX=1.559531175539e-3,
    lambdaY=2.205510139392e-3,
    lambdaT=1.0e-3,
    lambdaPx=6.41218345413e-4,
    lambdaPy=9.06819680526e-4,
    lambdaPt=1.0e-3,
)
sim.add_particles(bunch_charge_C, distr, npart)

# design the accelerator lattice
sol = elements.SoftSolenoid(
    name="sol1",
    ds=6.0,
    bscale=1.233482899483985,
    cos_coefficients=[
        0.350807812299706,
        0.323554693720069,
        0.260320578919415,
        0.182848575294969,
        0.106921016050403,
        4.409581845710694e-002,
        -9.416427163897508e-006,
        -2.459452716865687e-002,
        -3.272762575737291e-002,
        -2.936414401076162e-002,
        -1.995780078926890e-002,
        -9.102893342953847e-003,
        -2.456410658713271e-006,
        5.788233017324325e-003,
        8.040408292420691e-003,
        7.480064552867431e-003,
        5.230254569468851e-003,
        2.447614547094685e-003,
        -1.095525090532255e-006,
        -1.614586867387170e-003,
        -2.281365457438345e-003,
        -2.148709081338292e-003,
        -1.522541739363011e-003,
        -7.185505862719508e-004,
        -6.171194824600157e-007,
        4.842109305036943e-004,
        6.874508102002901e-004,
        6.535550288205728e-004,
        4.648795813759210e-004,
        2.216564722797528e-004,
        -4.100982995210341e-007,
        -1.499332112463395e-004,
        -2.151538438342482e-004,
        -2.044590946652016e-004,
        -1.468242784844341e-004,
    ],
    sin_coefficients=[
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
        0,
    ],
    mapsteps=200,
    nslice=4,
)

# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")

sim.lattice.extend(
    [
        monitor,
        sol,
        monitor,
    ]
)

# run simulation
sim.track_particles()

# clean shutdown
sim.finalize()

Analyze

We run the following script to analyze correctness:

Script analysis_solenoid_softedge_solvable.py
Listing 51 You can copy this file from examples/solenoid_softedge/analysis_solenoid_softedge_solvable.py.
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#

import numpy as np
import openpmd_api as io

# initial/final beam
series = io.Series("diags/openPMD/monitor.h5", io.Access.read_only)
last_step = list(series.iterations)[-1]
initial = series.iterations[1].particles["beam"].to_df()
beam_final = series.iterations[last_step].particles["beam"]
final = beam_final.to_df()

# compare number of particles
num_particles = 6
assert num_particles == len(initial)
assert num_particles == len(final)

# initial data
xi = np.array([1, 0, 0, 0, 0, 0])
pxi = np.array([0, 1, 0, 0, 0, 0])
yi = np.array([0, 0, 1, 0, 0, 0])
pyi = np.array([0, 0, 0, 1, 0, 0])
ti = np.array([0, 0, 0, 0, 1, 0])
pti = np.array([0, 0, 0, 0, 0, 1])

# problem parameters
g = 0.1
L = 2.0
bscale = -1.0  # This coincides with the value set in the input file.
gamma_ref = beam_final.get_attribute("gamma_ref")

# derived parameters
lbda = L / (2.0 * g)
phi = 2.0 * np.arctan(lbda)
b = 0.5 * bscale * g
beta = np.sqrt(1.0 + b**2)
cos1 = np.cos(b * phi)
cos2 = np.cos(beta * phi)
sin1 = np.sin(b * phi)
sin2 = np.sin(beta * phi)
bg = np.sqrt(gamma_ref**2 - 1.0)

# exact transfer matrix
r11 = cos1 * (cos2 + lbda * sin2 / beta)
r12 = g * (1.0 + lbda**2) * cos1 * sin2 / beta
r13 = sin1 * (cos2 + lbda * sin2 / beta)
r14 = g * (1.0 + lbda**2) * sin1 * sin2 / beta
r21 = (
    cos1
    * (2.0 * beta * lbda * cos2 + (lbda**2 - beta**2) * sin2)
    / (g * beta * (1.0 + lbda**2))
)
r22 = cos1 * (cos2 + lbda * sin2 / beta)
r23 = (
    sin1
    * (2.0 * beta * lbda * cos2 + (lbda**2 - beta**2) * sin2)
    / (g * beta * (1.0 + lbda**2))
)
r24 = sin1 * (cos2 + lbda * sin2 / beta)
r31 = -r13
r32 = -r14
r33 = r22
r34 = r12
r41 = -r23
r42 = -r24
r43 = r21
r44 = r22
r55 = 1.0
r56 = 2.0 * lbda * g / bg**2
r66 = 1.0

# final data
xf = np.array([r11, r12, r13, r14, 0, 0])
pxf = np.array([r21, r22, r23, r24, 0, 0])
yf = np.array([r31, r32, r33, r34, 0, 0])
pyf = np.array([r41, r42, r43, r44, 0, 0])
tf = np.array([0, 0, 0, 0, 1.0, r56])
ptf = np.array([0, 0, 0, 0, 0, 1.0])

# compute differences
error_xi = np.max(np.abs(xi - initial["position_x"].to_numpy()))
error_pxi = np.max(np.abs(pxi - initial["momentum_x"].to_numpy()))
error_yi = np.max(np.abs(yi - initial["position_y"].to_numpy()))
error_pyi = np.max(np.abs(pyi - initial["momentum_y"].to_numpy()))
error_ti = np.max(np.abs(ti - initial["position_t"].to_numpy()))
error_pti = np.max(np.abs(pti - initial["momentum_t"].to_numpy()))

error_xf = np.max(np.abs(xf - final["position_x"].to_numpy()))
error_pxf = np.max(np.abs(pxf - final["momentum_x"].to_numpy()))
error_yf = np.max(np.abs(yf - final["position_y"].to_numpy()))
error_pyf = np.max(np.abs(pyf - final["momentum_y"].to_numpy()))
error_tf = np.max(np.abs(tf - final["position_t"].to_numpy()))
error_ptf = np.max(np.abs(ptf - final["momentum_t"].to_numpy()))

xf_max = np.max(np.abs(xf))
pxf_max = np.max(np.abs(pxf))
yf_max = np.max(np.abs(yf))
pyf_max = np.max(np.abs(pyf))
tf_max = np.max(np.abs(tf))
ptf_max = np.max(np.abs(ptf))

print("Initial maximum absolute difference in each column:")
print("Difference x, px, y, py, t, pt:")
print(error_xi)
print(error_pxi)
print(error_yi)
print(error_pyi)
print(error_ti)
print(error_pti)

atol = 1.0e-13
print(f"  atol={atol}")

assert np.allclose(
    [error_xi, error_pxi, error_yi, error_pyi, error_ti, error_pti],
    [
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
    ],
    atol=atol,
)

print("")
print("Final maximum relative difference in each column:")
print("Difference x, px, y, py, t, pt:")
print(error_xf / xf_max)
print(error_pxf / pxf_max)
print(error_yf / yf_max)
print(error_pyf / pyf_max)
print(error_tf / tf_max)
print(error_ptf / ptf_max)

atol = 1.0e-7
print(f"  tol={atol}")

assert np.allclose(
    [
        error_xf / xf_max,
        error_pxf / pxf_max,
        error_yf / yf_max,
        error_pyf / pyf_max,
        error_tf / tf_max,
        error_ptf / ptf_max,
    ],
    [
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
    ],
    atol=atol,
)