Tune Calculation in a Periodic FODO Channel
This is identical to the FODO example, except that tracking for 100 periods is used to extract the horizontal tune.
Stable FODO cell with a zero-current phase advance of 67.8 degrees, corresponding to Qx = Qy = 0.1883.
The matched Twiss parameters at entry are:
\(\beta_\mathrm{x} = 2.82161941\) m
\(\alpha_\mathrm{x} = -1.59050035\)
\(\beta_\mathrm{y} = 2.82161941\) m
\(\alpha_\mathrm{y} = 1.59050035\)
We use a 2 GeV electron beam with initial unnormalized rms emittance of 2 nm.
The horizontal tune of a single particle is obtained from period-by-period tracking data using several algorithms.
In this test, the computed horizontal tune must agree with the nominal value to within acceptable tolerance.
This example requires installation of PyNAFF.
Run
This example can be run either as:
Python script:
python3 run_fodo_tune.pyorImpactX executable using an input file:
impactx input_fodo_tune.in
For MPI-parallel runs, prefix these lines with mpiexec -n 4 ... or srun -n 4 ..., depending on the system.
#!/usr/bin/env python3
#
# Copyright 2022-2023 ImpactX contributors
# Authors: Axel Huebl, Chad Mitchell
# License: BSD-3-Clause-LBNL
#
# -*- coding: utf-8 -*-
import amrex.space3d as amr
from impactx import ImpactX, distribution, elements
sim = ImpactX()
# set numerical parameters and IO control
sim.particle_shape = 2 # B-spline order
sim.space_charge = False
# sim.diagnostics = False # benchmarking
sim.slice_step_diagnostics = False
# domain decomposition & space charge mesh
sim.init_grids()
# load a 2 GeV electron beam with an initial
# unnormalized rms emittance of 2 nm
kin_energy_MeV = 2.0e3 # 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(0.510998950).set_kin_energy_MeV(kin_energy_MeV)
# particle bunch
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, distr, npart)
# add beam diagnostics
monitor = elements.BeamMonitor("monitor", backend="h5")
# design the accelerator lattice)
ns = 1 # number of slices per ds in the element
fodo = [
monitor,
elements.Drift(ds=0.25, nslice=ns),
elements.Quad(ds=1.0, k=1.0, nslice=ns),
elements.Drift(ds=0.5, nslice=ns),
elements.Quad(ds=1.0, k=-1.0, nslice=ns),
elements.Drift(ds=0.25, nslice=ns),
]
# assign a fodo segment
sim.lattice.extend(fodo)
# number of periods
sim.periods = 100
# run simulation
sim.evolve()
# clean shutdown
del sim
amr.finalize()
###############################################################################
# Particle Beam(s)
###############################################################################
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 = 3.9984884770e-5
beam.lambdaT = 1.0e-3
beam.lambdaPx = 2.6623538760e-5
beam.lambdaPy = 2.6623538760e-5
beam.lambdaPt = 2.0e-3
beam.muxpx = -0.846574929020762
beam.muypy = 0.846574929020762
beam.mutpt = 0.0
###############################################################################
# Beamline: lattice elements and segments
###############################################################################
lattice.periods = 100
lattice.elements = monitor drift1 quad1 drift2 quad2 drift3
lattice.nslice = 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.particle_shape = 2
algo.space_charge = false
###############################################################################
# Diagnostics
###############################################################################
diag.slice_step_diagnostics = false
Analyze
We run the following script to analyze correctness:
Script analysis_fodo_tune.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
import PyNAFF as pnf
# Collect beam data series
series = io.Series("diags/openPMD/monitor.h5", io.Access.read_only)
# Specify time series for particle j
j = 5
print(f"output for particle index = {j}")
# Create array of TBT data values
x = []
px = []
n = 0
for k_i, i in series.iterations.items():
beam = i.particles["beam"]
turn = beam.to_df()
x.append(turn["position_x"][j])
px.append(turn["momentum_x"][j])
n = n + 1
# Output number of periods in data series
nturns = len(x)
print(f"number of periods = {nturns}")
print()
# Approximate the tune and closed orbit using the 4-turn formula:
# from x data only
argument = (x[0] - x[1] + x[2] - x[3]) / (2.0 * (x[1] - x[2]))
tune1 = np.arccos(argument) / (2.0 * np.pi)
print(f"tune output from 4-turn formula, using x data = {tune1}")
# from px data only
argument = (px[0] - px[1] + px[2] - px[3]) / (2.0 * (px[1] - px[2]))
tune2 = np.arccos(argument) / (2.0 * np.pi)
print(f"tune output from 4-turn formula, using px data = {tune2}")
# orbit offset
xco = (x[1] ** 2 - x[2] ** 2 + x[1] * x[3] - x[2] * x[0]) / (
3.0 * (x[1] - x[2]) + x[3] - x[0]
)
pxco = (px[1] ** 2 - px[2] ** 2 + px[1] * px[3] - px[2] * px[0]) / (
3.0 * (px[1] - px[2]) + px[3] - px[0]
)
print(f"orbit offset from 4-turn formula: (x[m], px[rad]) = ({xco},{pxco})")
# matched Twiss functions from the 4-turn formula:
C1 = px[0] * (x[1] - x[2]) + px[1] * (x[2] - x[0]) + px[2] * (x[0] - x[1])
C2 = -px[0] * (x[1] - x[2]) + px[1] * (x[2] + x[0]) + px[2] * (-x[0] - x[1])
C3 = px[0] * (-x[1] - x[2]) - px[1] * (x[2] - x[0]) + px[2] * (x[0] + x[1])
C4 = px[0] * (x[1] + x[2]) + px[1] * (-x[2] - x[0]) - px[2] * (x[0] - x[1])
C = -C1 * C2 * C3 * C4
alphax = (
px[0] ** 2 * (x[1] ** 2 - x[2] ** 2)
+ px[1] ** 2 * (x[2] ** 2 - x[0] ** 2)
+ px[2] ** 2 * (x[0] ** 2 - x[1] ** 2)
) / np.sqrt(C)
betax = (
2.0
* (
px[0] * x[0] * (x[2] ** 2 - x[1] ** 2)
+ px[1] * x[1] * (x[0] ** 2 - x[2] ** 2)
+ px[2] * x[2] * (x[1] ** 2 - x[0] ** 2)
)
/ np.sqrt(C)
)
alphax = np.sign(betax) * alphax # ensure selection of correct sign
betax = np.sign(betax) * betax # ensure selection of correct sign
print(f"Twiss from 4-turn formula: alphax, betax [m] = {alphax}, {betax}")
# Normalize TBT data using Twiss functions:
z = []
for n in range(0, nturns):
xn = x[n] / np.sqrt(betax)
pxn = px[n] * np.sqrt(betax) + x[n] * alphax / np.sqrt(betax)
z.append(complex(xn, -pxn))
print()
# Approximate the tune using average phase advance (APA) - 1 period
tune3 = np.angle(z[1] / z[0]) / (2.0 * np.pi)
print(f"tune from phase advance = {tune3}")
print()
# Approximate the tune by using NAFF on the entire data series:
output = pnf.naff(
x, turns=nturns, nterms=4, skipTurns=0, getFullSpectrum=True, window=1
)
print("tune output from NAFF, using x data:")
print(output[0, 1], output[1, 1])
output = pnf.naff(
px, turns=nturns, nterms=4, skipTurns=0, getFullSpectrum=True, window=1
)
print("tune output from NAFF, using px data:")
print(output[0, 1], output[1, 1])
output = pnf.naff(
z, turns=nturns, nterms=4, skipTurns=0, getFullSpectrum=True, window=1
)
tune4 = output[0, 1]
print(f"tune output from NAFF, using normalized x-ipx data: {tune4}")
rtol = 1.0e-3
print(f" rtol={rtol}")
assert np.allclose(
[alphax, betax, tune1, tune2, tune3, tune4],
[
-1.59050035,
2.82161941,
0.1883,
0.1883,
0.1883,
0.1883,
],
rtol=rtol,
)