diff --git a/__EXAMPLES/main_files/EX_23_single_particle_with_oneturnmatrix_lhc.py b/__EXAMPLES/main_files/EX_23_single_particle_with_oneturnmatrix_lhc.py
new file mode 100644
index 0000000000000000000000000000000000000000..6edecffd1809620fdd5be4e3b9dddb62a806a441
--- /dev/null
+++ b/__EXAMPLES/main_files/EX_23_single_particle_with_oneturnmatrix_lhc.py
@@ -0,0 +1,225 @@
+'''
+LHC simulation of a single particle.
+
+:Author: **Thom Arnoldus van Rijswijk**
+'''
+
+# Imports
+import os
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+
+# Xsuite imports
+import xpart as xp
+import xtrack as xt
+import xplt
+
+# BLonD objects
+from blond.beam.beam import Beam, Proton
+from blond.input_parameters.ring import Ring
+from blond.input_parameters.rf_parameters import RFStation
+from blond.trackers.tracker import RingAndRFTracker
+
+# Interface
+from blond.interfaces.xsuite import (BlondElement, BlondObserver,
+ blond_to_xsuite_transform,
+ xsuite_to_blond_transform)
+
+# Monitor objects
+from blond.monitors.monitors import BunchMonitor
+from blond.plots.plot import Plot
+
+# Directory to save files
+this_directory = os.path.dirname(os.path.realpath(__file__)) + '/'
+results_directory = this_directory + '../output_files/EX_23'
+os.makedirs(results_directory, exist_ok=True)
+
+# Parameters ----------------------------------------------------------------------------------------------------------
+# Accelerator parameters
+C = 26658.8832 # Machine circumference [m]
+p_s = 450e9 # Synchronous momentum [eV/c]
+p_f = 450.0e9 # Synchronous momentum, final
+h = 35640 # Harmonic number [-]
+alpha = 0.00034849575112251314 # First order mom. comp. factor [-]
+V = 5e6 # RF voltage [V]
+dphi = 0 # Phase modulation/offset [rad]
+
+# Bunch parameters
+N_bunches = 1 # Number of bunches [-]
+N_m = 1 # Number of macroparticles [-]
+N_p = 1.15e11 # Intensity
+blen = 1.25e-9 # Bunch length [s]
+energy_err = 100e6 # Beamenergy error [eV/c]
+
+# Simulation parameters
+N_t = 481 # Number of (tracked) turns [-]
+N_buckets = 2 # Number of buckets [-]
+dt_plt = 50 # Timestep between plots [-]
+input_dt = 2 * blen - 0.1e-9 # Input particles dt [s]
+input_dE = 0.0 # Input particles dE [eV]
+
+# Preparing BLonD items -----------------------------------------------------------------------------------------------
+print(f"\nPreparing BLonD items...")
+# --- BLonD objects ---
+ring = Ring(C, alpha, np.linspace(p_s, p_f, N_t + 1), Proton(), n_turns=N_t)
+rfstation = RFStation(ring, [h], [V], [dphi])
+beam = Beam(ring, N_m, N_p)
+
+# --- Insert particle ---
+beam.dt = np.array([input_dt])
+beam.dE = np.array([input_dE])
+
+# --- RF tracker object ---
+blond_tracker = RingAndRFTracker(rfstation, beam, Profile=None, interpolation=False,
+ CavityFeedback=None, BeamFeedback=None,
+ TotalInducedVoltage=None,
+ with_xsuite=True)
+
+# --- Bunch monitor ---
+# Print parameters for plotting
+print(f"\n Calculating plotting parameters...")
+
+t_u = np.pi / rfstation.omega_rf[0, 0]
+Y_phi_s = np.sqrt(abs(-np.cos(rfstation.phi_s[0]) + (np.pi - 2 * rfstation.phi_s[0]) / 2 * np.sin(rfstation.phi_s[0])))
+dE_sep = np.sqrt(2 * abs(rfstation.charge) * rfstation.voltage[0, 0] * rfstation.beta[0] * rfstation.beta[0]
+ * rfstation.energy[0] / (np.pi * rfstation.harmonic[0, 0] * abs(ring.eta_0[0, 0]))) * Y_phi_s
+
+print('t_u: ' + str(2 * t_u))
+print('dE_sep,m: ' + str(dE_sep))
+
+# Make bunchmonitor
+bunchmonitor = BunchMonitor(ring, rfstation, beam,
+ results_directory + '/EX_23_output_data', Profile=None)
+
+format_options = {'dirname': results_directory + '/EX_23_fig'}
+plots = Plot(ring, rfstation, beam, dt_plt, N_t, 0, 2 * t_u,
+ -1.05 * dE_sep, 1.05 * dE_sep, xunit='s', separatrix_plot=True,
+ Profile=None, h5file=results_directory + '/EX_23_output_data',
+ format_options=format_options)
+
+# --- Creating interface elements ---
+# Blond tracker
+cavity = BlondElement(blond_tracker, beam)
+# BLonD monitor
+beam_monitor = BlondElement(bunchmonitor, beam)
+beam_plots = BlondElement(plots, beam)
+
+# Preparing Xsuite items -----------------------------------------------------------------------------------------------
+print(f"\nPreparing Xsuite items...")
+# --- Setup matrix ---
+# Make First order matrix map (takes care of drift in Xsuite)
+matrix = xt.LineSegmentMap(
+ longitudinal_mode='nonlinear',
+ qx=1.1, qy=1.2,
+ betx=1.,
+ bety=1.,
+ voltage_rf=0,
+ frequency_rf=0,
+ lag_rf=0,
+ momentum_compaction_factor=alpha,
+ length=C
+)
+
+# Create line
+line = xt.Line(elements=[matrix], element_names={'matrix'} )
+line['matrix'].length = C
+
+# Insert the BLonD elements
+line.insert_element(index=0, element=cavity, name='blond_cavity')
+line.append_element(element=beam_monitor, name='blond_monitor')
+line.append_element(element=beam_plots, name='blond_plots')
+
+# Add particles to line and build tracker
+line.particle_ref = xp.Particles(p0c=p_s, mass0=xp.PROTON_MASS_EV, q0=1.)
+line.build_tracker()
+
+# Show table
+line.get_table().show()
+# prints:
+
+# name s element_type isthick isreplica parent_name iscollective
+# blond_cavity 0 BlondElement False False None True
+# matrix 0 LineSegmentMap True False None False
+# blond_monitor 26658.9 BlondObserver False False None True
+# blond_plots 26658.9 BlondElement False False None True
+# _end_point 26658.9 False False None False
+
+# Simulating ----------------------------------------------------------------------------------------------------------
+print(f"\nSetting up simulation...")
+
+# --- Convert the initial BLonD distribution to xsuite coordinates ---
+zeta, ptau = blond_to_xsuite_transform(
+ beam.dt, beam.dE, line.particle_ref.beta0[0], line.particle_ref.energy0[0],
+ phi_s=rfstation.phi_s[0], omega_rf=rfstation.omega_rf[0, 0]
+)
+
+# --- Track matrix ---
+particles = line.build_particles(x=0, y=0, px=0, py=0, zeta=np.copy(zeta), ptau=np.copy(ptau))
+line.track(particles, num_turns=N_t, turn_by_turn_monitor=True, with_progress=True)
+mon = line.record_last_track
+
+# Saving turn-by-turn particle coordinates
+test_string = ''
+test_string += '{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\n'.format(
+ 'dE', 'dt', 'momentum', 'gamma', 'beta', 'energy', 'x', 'px', 'y', 'py')
+test_string += ('{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}'
+ '\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\n').format(
+ beam.dE[0], beam.dt[0], rfstation.momentum[0], rfstation.gamma[0], rfstation.beta[0],
+ rfstation.energy[0], mon.x[0, 0].T, mon.px[0, 0].T, mon.y[0, 0].T, mon.py[0, 0].T)
+
+# Saving turn-by-turn Xsuite particle coordinates
+xsuite_string = ''
+xsuite_string += '{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\n'.format(
+ 'ptau', 'zeta', 'momentum', 'gamma', 'beta')
+xsuite_string += '{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\n'.format(
+ mon.ptau[0, 0].T, mon.zeta[0, 0].T, np.float64(mon.p0c[0, 0]), np.float64(mon.gamma0[0, 0]),
+ np.float64(mon.beta0[0, 0]))
+
+# Convert the xsuite particle coordinates back to BLonD
+for i in range(N_t):
+ dt, dE = xsuite_to_blond_transform(
+ mon.zeta[:, i].T, mon.ptau[:, i].T,
+ rfstation.beta[i], rfstation.energy[i],
+ phi_s=rfstation.phi_s[i] - rfstation.phi_rf[0, 0], # Correct for RF phase
+ omega_rf=rfstation.omega_rf[0, i]
+ )
+
+ # Statistics
+ test_string += ('{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}'
+ '\t{:+10.10e}\t{:+10.10e}\n').format(
+ dE[0], dt[0], rfstation.momentum[i], rfstation.gamma[i], rfstation.beta[i],
+ rfstation.energy[i], mon.x[0, i].T, mon.px[0, i].T, mon.y[0, i].T, mon.py[0, i].T)
+ xsuite_string += '{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\n'.format(
+ mon.ptau[0, i].T, mon.zeta[0, i].T, np.float64(mon.p0c[0, i]), np.float64(mon.gamma0[0, i]),
+ np.float64(mon.beta0[0, i]))
+
+with open(results_directory + '/EX_23_output_data.txt', 'w') as f:
+ f.write(test_string)
+with open(results_directory + '/EX_23_output_data_xsuite.txt', 'w') as f:
+ f.write(xsuite_string)
+
+# Results --------------------------------------------------------------------------------------------------------------
+print(f"\nPlotting result...")
+
+# Plot Phasespace
+
+plt.figure()
+plt.plot(mon.zeta.T, mon.ptau.T)
+plt.scatter(np.mean(mon.zeta.T), np.mean(mon.ptau.T), label='mon')
+plt.grid()
+plt.title('Phasespace')
+plt.xlabel(r'$\zeta$')
+plt.ylabel(r'$p_{\tau}$')
+plt.savefig(results_directory + '/figure_Phasespace.png')
+
+# Use Xplt to plot Phasespace
+
+plot = xplt.PhaseSpacePlot(
+ mon
+)
+plot.fig.suptitle("Particle distribution for a single turn")
+plot.fig.savefig(results_directory + '/figure_ps_mon.png')
+
+# Show plots
+plt.show()
diff --git a/__EXAMPLES/main_files/EX_24_single_particle_ramp_lhc.py b/__EXAMPLES/main_files/EX_24_single_particle_ramp_lhc.py
new file mode 100644
index 0000000000000000000000000000000000000000..cfaa4644e1ee8134658a55112fc0b85dae37277c
--- /dev/null
+++ b/__EXAMPLES/main_files/EX_24_single_particle_ramp_lhc.py
@@ -0,0 +1,237 @@
+'''
+LHC simulation of a single particle with a ramp.
+
+: Author: **Birk Emil Karlsen-Baeck**, **Thom Arnoldus van Rijswijk**
+'''
+
+# Imports
+import os
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+#mpl.use('Agg')
+
+# Xsuite imports
+import xpart as xp
+import xtrack as xt
+import xplt
+
+# BLonD objects
+from blond.beam.beam import Beam, Proton
+from blond.input_parameters.ring import Ring
+from blond.input_parameters.rf_parameters import RFStation
+from blond.trackers.tracker import RingAndRFTracker
+
+# Interface
+from blond.interfaces.xsuite import (BlondElement, BlondObserver, EnergyUpdate,
+ blond_to_xsuite_transform,
+ xsuite_to_blond_transform)
+
+# Monitor objects
+from blond.monitors.monitors import BunchMonitor
+from blond.plots.plot import Plot
+
+# Directory to save files
+this_directory = os.path.dirname(os.path.realpath(__file__)) + '/'
+results_directory = this_directory + '../output_files/EX_24'
+os.makedirs(results_directory, exist_ok=True)
+
+# Parameters ----------------------------------------------------------------------------------------------------------
+# Accelerator parameters
+C = 26658.8832 # Machine circumference [m]
+p_s = 450e9 # Synchronous momentum [eV/c]
+p_f = 450.1e9 # Synchronous momentum, final
+h = 35640 # Harmonic number [-]
+alpha = 0.00034849575112251314 # First order mom. comp. factor [-]
+V = 5e6 # RF voltage [V]
+dphi = 0 # Phase modulation/offset [rad]
+
+# Bunch parameters
+N_bunches = 1 # Number of bunches [-]
+N_m = 1 # Number of macroparticles [-]
+N_p = 1.15e11 # Intensity
+blen = 1.25e-9 # Bunch length [s]
+energy_err = 100e6 # Beamenergy error [eV/c]
+
+# Simulation parameters
+N_t = 330 # Number of (tracked) turns [-]
+N_buckets = 2 # Number of buckets [-]
+dt_plt = 30 # Timestep between plots [-]
+input_dt = 2 * blen - 0.4e-9 # Input particles dt [s]
+input_dE = 0.0 # Input particles dE [eV]
+
+# Xsuite parameters
+cavity_name = 'acsca.d5l4.b1'
+
+# Preparing BLonD items -----------------------------------------------------------------------------------------------
+print(f"\nPreparing BLonD items...")
+# --- BLonD objects ---
+ring = Ring(C, alpha, np.linspace(p_s, p_f, N_t + 1), Proton(), n_turns=N_t)
+rfstation = RFStation(ring, [h], [V], [dphi])
+beam = Beam(ring, N_m, N_p)
+
+# --- Insert particle ---
+beam.dt = np.array([input_dt])
+beam.dE = np.array([input_dE])
+
+# --- RF tracker object ---
+blond_tracker = RingAndRFTracker(rfstation, beam, Profile=None, interpolation=False,
+ CavityFeedback=None, BeamFeedback=None,
+ TotalInducedVoltage=None,
+ with_xsuite=True)
+
+# --- Bunch monitor ---
+# Print parameters for plotting
+print(f"\n Calculating plotting parameters...")
+
+t_u = np.pi / rfstation.omega_rf[0, 0]
+Y_phi_s = np.sqrt(abs(-np.cos(rfstation.phi_s[0]) + (np.pi - 2 * rfstation.phi_s[0]) / 2 * np.sin(rfstation.phi_s[0])))
+dE_sep = np.sqrt(2 * abs(rfstation.charge) * rfstation.voltage[0, 0] * rfstation.beta[0] * rfstation.beta[0]
+ * rfstation.energy[0] / (np.pi * rfstation.harmonic[0, 0] * abs(ring.eta_0[0, 0]))) * Y_phi_s
+
+print('t_u: ' + str(2 * t_u))
+print('dE_sep,m: ' + str(dE_sep))
+
+# Make bunchmonitor
+bunchmonitor = BunchMonitor(ring, rfstation, beam,
+ results_directory + '/EX_24_output_data', Profile=None)
+
+format_options = {'dirname': results_directory + '/EX_24_fig'}
+plots = Plot(ring, rfstation, beam, dt_plt, N_t, 0, 2*t_u,
+ -1.05 * dE_sep, 1.05 * dE_sep, xunit='s', separatrix_plot=True,
+ Profile=None, h5file=results_directory + '/EX_24_output_data',
+ format_options=format_options)
+
+# --- Creating interface elements ---
+# Blond tracker
+cavity = BlondElement(blond_tracker, beam)
+# BLonD monitor
+beam_monitor = BlondElement(bunchmonitor, beam)
+beam_plots = BlondElement(plots, beam)
+
+# Preparing Xsuite items -----------------------------------------------------------------------------------------------
+print(f"\nPreparing Xsuite items...")
+# --- Setup matrix ---
+# Make First order matrix map (takes care of drift in Xsuite)
+matrix = xt.LineSegmentMap(
+ longitudinal_mode='nonlinear',
+ qx=1.1, qy=1.2,
+ betx=1.,
+ bety=1.,
+ voltage_rf=0,
+ frequency_rf=0,
+ lag_rf=0,
+ momentum_compaction_factor=alpha,
+ length=C)
+
+# Create line
+line = xt.Line(elements=[matrix], element_names={'matrix'})
+line['matrix'].length = C
+
+# Insert the BLonD elements
+line.insert_element(index=0, element=cavity, name='blond_cavity')
+line.append_element(element=beam_monitor, name='blond_monitor')
+line.append_element(element=beam_plots, name='blond_plots')
+
+# Insert energy ramp
+energy_update = EnergyUpdate(ring.momentum[0, :])
+line.insert_element(index='matrix', element=energy_update, name='energy_update')
+
+# Add particles to line and build tracker
+line.particle_ref = xp.Particles(p0c=p_s, mass0=xp.PROTON_MASS_EV, q0=1.)
+line.build_tracker()
+
+# Show table
+line.get_table().show()
+# prints:
+
+# name s element_type isthick isreplica parent_name iscollective
+# energy_update 0 EnergyUpdate False False None True
+# blond_cavity 0 BlondElement False False None True
+# matrix 0 LineSegmentMap True False None False
+# blond_monitor 26658.9 BlondObserver False False None True
+# blond_plots 26658.9 BlondElement False False None True
+# _end_point 26658.9 False False None False
+
+# Simulating ----------------------------------------------------------------------------------------------------------
+print(f"\nSetting up simulation...")
+
+# --- Convert the initial BLonD distribution to xsuite coordinates ---
+zeta, ptau = blond_to_xsuite_transform(
+ beam.dt, beam.dE, line.particle_ref.beta0[0],
+ line.particle_ref.energy0[0], phi_s=rfstation.phi_s[0],
+ omega_rf=rfstation.omega_rf[0, 0]
+)
+
+# --- Track matrix ---
+particles = line.build_particles(x=0, y=0, px=0, py=0, zeta=np.copy(zeta), ptau=np.copy(ptau))
+line.track(particles, num_turns=N_t, turn_by_turn_monitor=True, with_progress=True)
+mon = line.record_last_track
+
+# Saving turn-by-turn particle coordinates
+test_string = ''
+test_string += '{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\n'.format(
+ 'dE', 'dt', 'momentum', 'gamma', 'beta', 'energy', 'x', 'px', 'y', 'py')
+test_string += ('{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}'
+ '\t{:+10.10e}\t{:+10.10e}\n').format(beam.dE[0], beam.dt[0], rfstation.momentum[0],
+ rfstation.gamma[0], rfstation.beta[0], rfstation.energy[0],
+ mon.x[0, 0].T, mon.px[0, 0].T, mon.y[0, 0].T, mon.py[0, 0].T)
+
+# Saving turn-by-turn Xsuite particle coordinates
+xsuite_string = ''
+xsuite_string += '{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\n'.format(
+ 'ptau', 'zeta', 'momentum', 'gamma', 'beta')
+xsuite_string += '{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\n'.format(
+ mon.ptau[0, 0].T, mon.zeta[0, 0].T, np.float64(mon.p0c[0, 0]), np.float64(mon.gamma0[0, 0]),
+ np.float64(mon.beta0[0, 0])
+)
+
+# Convert the xsuite particle coordinates back to BLonD
+for i in range(N_t):
+ dt, dE = xsuite_to_blond_transform(
+ mon.zeta[:, i].T, mon.ptau[:, i].T,
+ rfstation.beta[i], rfstation.energy[i],
+ phi_s=rfstation.phi_s[i] - rfstation.phi_rf[0, 0], # Correct for RF phase
+ omega_rf=rfstation.omega_rf[0, i]
+ )
+
+ # Statistics
+ test_string += ('{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}'
+ '\t{:+10.10e}\t{:+10.10e}\n').format(dE[0], dt[0], rfstation.momentum[i], rfstation.gamma[i],
+ rfstation.beta[i], rfstation.energy[i], mon.x[0, i].T,
+ mon.px[0, i].T, mon.y[0, i].T, mon.py[0, i].T)
+
+ xsuite_string += '{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\n'.format(
+ mon.ptau[0, i].T, mon.zeta[0, i].T, np.float64(mon.p0c[0, i]), np.float64(mon.gamma0[0, i]),
+ np.float64(mon.beta0[0, i])
+ )
+
+with open(results_directory + '/EX_24_output_data.txt', 'w') as f:
+ f.write(test_string)
+with open(results_directory + '/EX_24_output_data_xsuite.txt', 'w') as f:
+ f.write(xsuite_string)
+
+# Results --------------------------------------------------------------------------------------------------------------
+print(f"\nPlotting result...")
+
+# Plot Phasespace
+
+plt.figure()
+plt.plot(mon.zeta.T,mon.ptau.T)
+plt.scatter(np.mean(mon.zeta.T), np.mean(mon.ptau.T), label='mon')
+plt.grid()
+plt.title('Phasespace')
+plt.xlabel(r'$\zeta$')
+plt.ylabel(r'$p_{\tau}$')
+plt.savefig(results_directory + '/figure_Phasespace.png')
+
+# Use Xplt to plot Phasespace
+
+plot = xplt.PhaseSpacePlot(
+ mon
+)
+plot.fig.suptitle("Particle distribution for a single turn")
+plot.fig.savefig(results_directory + '/figure_ps_mon.png')
+
+# Show plots
+plt.show()
diff --git a/__EXAMPLES/main_files/EX_25_single_particle_two_rfstations_lhc.py b/__EXAMPLES/main_files/EX_25_single_particle_two_rfstations_lhc.py
new file mode 100644
index 0000000000000000000000000000000000000000..9103a32f0d2c00928389cd0445ae94ddde72eb24
--- /dev/null
+++ b/__EXAMPLES/main_files/EX_25_single_particle_two_rfstations_lhc.py
@@ -0,0 +1,256 @@
+'''
+LHC simulation of a single particle with two RF stations.
+
+: Author: **Thom Arnoldus van Rijswijk**
+'''
+
+# Imports
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Xsuite imports
+import xpart as xp
+import xtrack as xt
+import xplt
+
+# BLonD objects
+from blond.beam.beam import Beam, Proton
+from blond.input_parameters.ring import Ring
+from blond.input_parameters.rf_parameters import RFStation
+from blond.trackers.tracker import RingAndRFTracker
+
+# Interface
+from blond.interfaces.xsuite import (BlondElement,
+ blond_to_xsuite_transform,
+ xsuite_to_blond_transform)
+
+# Monitor objects
+from blond.monitors.monitors import BunchMonitor
+from blond.plots.plot import Plot
+
+# Directory to save files
+this_directory = os.path.dirname(os.path.realpath(__file__)) + '/'
+results_directory = this_directory + '../output_files/EX_25'
+os.makedirs(results_directory, exist_ok=True)
+
+# Parameters ----------------------------------------------------------------------------------------------------------
+# Accelerator parameters
+C = 26658.8832 # Machine circumference [m]
+p_s = 450e9 # Synchronous momentum [eV/c]
+p_f = 450.0e9 # Synchronous momentum, final
+h = 35640 # Harmonic number [-]
+alpha = 0.00034849575112251314 # First order mom. comp. factor [-]
+V = 5e6 # RF voltage [V]
+dphi = 0 # Phase modulation/offset [rad]
+
+# Bunch parameters
+N_bunches = 1 # Number of bunches [-]
+N_m = 1 # Number of macroparticles [-]
+N_p = 1.15e11 # Intensity
+blen = 1.25e-9 # Bunch length [s]
+energy_err = 100e6 # Beamenergy error [eV/c]
+
+# Simulation parameters
+N_t = 481 # Number of (tracked) turns [-]
+N_buckets = 2 # Number of buckets [-]
+dt_plt = 50 # Timestep between plots [-]
+input_dt = 2 * blen - 0.1e-9 # Input particles dt [s]
+input_dE = 0.0 # Input particles dE [eV]
+
+# Additional input
+C_part = 0.3
+V_part = 0.5
+
+# Preparing BLonD items -----------------------------------------------------------------------------------------------
+print(f"\nPreparing BLonD items...")
+# --- BLonD objects ---
+ring = Ring([C_part * C, (1 - C_part) * C],
+ [[alpha], [alpha]],
+ [p_s, p_s],
+ Proton(), n_turns=N_t, n_sections=2)
+rfstation_1 = RFStation(ring, [h], [V * V_part], [dphi], section_index=1)
+rfstation_2 = RFStation(ring, [h], [V * (1 - V_part)], [dphi], section_index=2)
+beam = Beam(ring, N_m, N_p)
+
+# --- Insert particle ---
+beam.dt = np.array([input_dt])
+beam.dE = np.array([input_dE])
+
+# --- RF tracker object ---
+rftracker_1 = RingAndRFTracker(rfstation_1, beam, Profile=None, interpolation=False,
+ CavityFeedback=None, BeamFeedback=None,
+ TotalInducedVoltage=None,
+ with_xsuite=True)
+rftracker_2 = RingAndRFTracker(rfstation_2, beam, Profile=None, interpolation=False,
+ CavityFeedback=None, BeamFeedback=None,
+ TotalInducedVoltage=None,
+ with_xsuite=True)
+
+# --- Bunch monitor ---
+# Print parameters for plotting
+print(f"\n Calculating plotting parameters...")
+
+t_u = np.pi / rfstation_1.omega_rf[0,0]
+Y_phi_s = np.sqrt(abs(-np.cos(rfstation_1.phi_s[0]) + (np.pi - 2*rfstation_1.phi_s[0])
+ / 2 * np.sin(rfstation_1.phi_s[0])))
+dE_sep = np.sqrt(2 * abs(rfstation_1.charge) * V * rfstation_1.beta[0] * rfstation_1.beta[0] * rfstation_1.energy[0]
+ / (np.pi * rfstation_1.harmonic[0, 0] * abs(ring.eta_0[0, 0]))) * Y_phi_s
+
+print('t_u: ' + str(2 * t_u))
+print('dE_sep,m: ' + str(dE_sep))
+
+# Make bunchmonitor
+bunchmonitor = BunchMonitor(ring, rfstation_1, beam,
+ results_directory + '/EX_25_output_data', Profile=None)
+
+format_options = {'dirname': results_directory + '/EX_25_fig'}
+plots = Plot(ring, rfstation_1, beam, dt_plt, N_t, 0, 2 * t_u,
+ -1.05 * dE_sep, 1.05 * dE_sep, xunit='s', separatrix_plot=True,
+ Profile=None, h5file=results_directory + '/EX_25_output_data',
+ format_options=format_options)
+
+# --- Creating interface elements ---
+# Blond tracker
+cavity_1 = BlondElement(rftracker_1, beam)
+cavity_2 = BlondElement(rftracker_2, beam)
+# BLonD monitor
+beam_monitor = BlondElement(bunchmonitor, beam)
+beam_plots = BlondElement(plots, beam)
+
+# Xsuite initialization -----------------------------------------------------------------------------------------------
+print(f"\nXsuite initialization...")
+# --- Setup matrices ---
+# Make First order matrix maps (takes care of drift in Xsuite)
+matrix_1 = xt.LineSegmentMap(
+ longitudinal_mode='nonlinear',
+ qx=1.1, qy=1.2,
+ betx=1.,
+ bety=1.,
+ voltage_rf=0,
+ frequency_rf=0,
+ lag_rf=0,
+ momentum_compaction_factor=alpha,
+ length=C * C_part)
+matrix_2 = xt.LineSegmentMap(
+ longitudinal_mode='nonlinear',
+ qx=1.1, qy=1.2,
+ betx=1.,
+ bety=1.,
+ voltage_rf=0,
+ frequency_rf=0,
+ lag_rf=0,
+ momentum_compaction_factor=alpha,
+ length=C * (1 - C_part))
+
+# Create new line
+line = xt.Line(elements=[matrix_1], element_names={'matrix_1'} )
+line['matrix_1'].length = C_part * C
+line.append_element(element=matrix_2, name='matrix_2')
+line['matrix_2'].length = (1 - C_part) * C
+
+# Insert the BLonD elements
+line.insert_element(index='matrix_1', element=cavity_1, name='blond_cavity_1')
+line.insert_element(index='matrix_2', element=cavity_2, name='blond_cavity_2')
+line.append_element(element=beam_monitor, name='blond_monitor')
+line.append_element(element=beam_plots, name='blond_plots')
+
+# Add particles to line and build tracker
+line.particle_ref = xp.Particles(p0c=p_s, mass0=xp.PROTON_MASS_EV, q0=1.)
+line.build_tracker()
+
+# Show table
+line.get_table().show()
+# prints:
+
+# name s element_type isthick isreplica parent_name iscollective
+# blond_cavity_1 0 BlondElement False False None True
+# matrix_1 0 LineSegmentMap True False None False
+# blond_cavity_2 7997.66 BlondElement False False None True
+# matrix_2 7997.66 LineSegmentMap True False None False
+# blond_monitor 26658.9 BlondObserver False False None True
+# blond_plots 26658.9 BlondElement False False None True
+# _end_point 26658.9 False False None False
+
+# Simulating ----------------------------------------------------------------------------------------------------------
+print(f"\nSetting up simulation...")
+
+# --- Convert the initial BLonD distribution to xsuite coordinates ---
+zeta, ptau = blond_to_xsuite_transform(
+ beam.dt, beam.dE, line.particle_ref.beta0[0], line.particle_ref.energy0[0],
+ phi_s=rfstation_1.phi_s[0], omega_rf=rfstation_1.omega_rf[0, 0]
+)
+
+# --- Track matrix ---
+particles = line.build_particles(x=0, y=0, px=0, py=0, zeta=np.copy(zeta), ptau=np.copy(ptau))
+line.track(particles, num_turns=N_t, turn_by_turn_monitor=True, with_progress=True)
+mon = line.record_last_track
+
+# Saving turn-by-turn particle coordinates
+test_string = ''
+test_string += '{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\n'.format(
+ 'dE', 'dt', 'momentum', 'gamma', 'beta', 'energy', 'x', 'px', 'y', 'py')
+test_string += ('{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}'
+ '\t{:+10.10e}\t{:+10.10e}\n').format(beam.dE[0], beam.dt[0], rfstation_1.momentum[0],
+ rfstation_1.gamma[0], rfstation_1.beta[0], rfstation_1.energy[0],
+ mon.x[0, 0].T, mon.px[0, 0].T, mon.y[0, 0].T, mon.py[0, 0].T)
+
+# Saving turn-by-turn Xsuite particle coordinates
+xsuite_string = ''
+xsuite_string += '{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\n'.format(
+ 'ptau', 'zeta', 'momentum', 'gamma', 'beta')
+xsuite_string += '{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\n'.format(
+ mon.ptau[0, 0].T, mon.zeta[0, 0].T, np.float64(mon.p0c[0, 0]), np.float64(mon.gamma0[0, 0]),
+ np.float64(mon.beta0[0, 0])
+)
+
+# Convert the xsuite particle coordinates back to BLonD
+for i in range(N_t):
+ dt, dE = xsuite_to_blond_transform(
+ mon.zeta[:,i].T, mon.ptau[:,i].T,
+ rfstation_1.beta[i], rfstation_1.energy[i],
+ phi_s=rfstation_1.phi_s[i] - rfstation_1.phi_rf[0, 0],
+ omega_rf=rfstation_1.omega_rf[0, i]
+ )
+
+ # Statistics
+ test_string += ('{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}'
+ '\t{:+10.10e}\t{:+10.10e}\n').format(dE[0], dt[0], rfstation_1.momentum[i],
+ rfstation_1.gamma[i], rfstation_1.beta[i],
+ rfstation_1.energy[i], mon.x[0, i].T, mon.px[0, i].T,
+ mon.y[0, i].T, mon.py[0, i].T)
+
+ xsuite_string += '{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\n'.format(
+ mon.ptau[0, i].T, mon.zeta[0, i].T, np.float64(mon.p0c[0, i]), np.float64(mon.gamma0[0, i]),
+ np.float64(mon.beta0[0, i])
+ )
+
+with open(results_directory + '/EX_25_output_data.txt', 'w') as f:
+ f.write(test_string)
+with open(results_directory + '/EX_25_output_data_xsuite.txt', 'w') as f:
+ f.write(xsuite_string)
+
+# Results --------------------------------------------------------------------------------------------------------------
+print(f"\nPlotting result...")
+
+# Plot Phasespace
+
+plt.figure()
+plt.plot(mon.zeta.T,mon.ptau.T)
+plt.scatter(np.mean(mon.zeta.T), np.mean(mon.ptau.T), label='mon')
+plt.grid()
+plt.title('Phasespace')
+plt.xlabel(r'$\zeta$')
+plt.ylabel(r'$p_{\tau}$')
+plt.savefig(results_directory + '/figure_Phasespace.png')
+
+# Use Xplt to plot Phasespace
+
+plot = xplt.PhaseSpacePlot(
+ mon
+)
+plot.fig.suptitle("Particle distribution for a single turn")
+plot.fig.savefig(results_directory + '/figure_ps_mon.png')
+
+# Show plots
+plt.show()
diff --git a/__EXAMPLES/main_files/EX_26_single_particle_oneturnmatrix_psb.py b/__EXAMPLES/main_files/EX_26_single_particle_oneturnmatrix_psb.py
new file mode 100644
index 0000000000000000000000000000000000000000..845ffa3c5fe5830780c9e7ba90d14a202d255ab1
--- /dev/null
+++ b/__EXAMPLES/main_files/EX_26_single_particle_oneturnmatrix_psb.py
@@ -0,0 +1,233 @@
+'''
+PSB simulation of a single particle.
+
+:Author: **Thom Arnoldus van Rijswijk**
+'''
+
+# Imports
+import os
+import numpy as np
+from scipy.constants import c, e, m_p
+import matplotlib.pyplot as plt
+
+# Xsuite imports
+import xpart as xp
+import xtrack as xt
+import xplt
+
+# BLonD objects
+from blond.beam.beam import Beam, Proton
+from blond.input_parameters.ring import Ring
+from blond.input_parameters.rf_parameters import RFStation
+from blond.trackers.tracker import RingAndRFTracker
+
+# Interface
+from blond.interfaces.xsuite import (BlondElement,
+ blond_to_xsuite_transform,
+ xsuite_to_blond_transform)
+
+# Monitor objects
+from blond.monitors.monitors import BunchMonitor
+from blond.plots.plot import Plot
+
+# Directory to save files
+this_directory = os.path.dirname(os.path.realpath(__file__)) + '/'
+results_directory = this_directory + '../output_files/EX_26'
+os.makedirs(results_directory, exist_ok=True)
+
+# Parameters ----------------------------------------------------------------------------------------------------------
+# Accelerator parameters
+C = 2 * np.pi * 25.000057945260654 # Machine circumference, radius 25 [m]
+gamma_t = 4.11635447373496 # Transition gamma [-]
+alpha = 1./gamma_t/gamma_t # First order mom. comp. factor [-]
+h = 1 # Harmonic number [-]
+V = 8e3 # RF voltage [V]
+dphi = np.pi # Phase modulation/offset [rad]
+
+# Beam parameters
+N_m = 1 # Number of macroparticles [-]
+N_p = 1e11 # Intensity / number of particles
+sigma_dt = 180e-9 / 4 # [s]
+blen = sigma_dt * 4 # Initial bunch length, 4 sigma [s]
+kin_beam_energy = 1.4e9 # [eV]
+
+
+# Simulation parameters
+N_t = 6800 # Number of (tracked) turns [-]
+N_buckets = 1 # Number of buckets [-]
+dt_plt = 500 # Timestep between plots [-]
+input_dt = 5.720357923415153e-07 - sigma_dt # Input particles dt [s]
+input_dE = 0.0 # Input particles dE [eV]
+
+# Derived parameters
+E_0 = m_p * c**2 / e # [eV]
+tot_beam_energy = E_0 + kin_beam_energy # [eV]
+sync_momentum = np.sqrt(tot_beam_energy**2 - E_0**2) # [eV / c]
+
+# Preparing BLonD items -----------------------------------------------------------------------------------------------
+print(f"\nPreparing BLonD items...")
+# --- BLonD objects ---
+ring = Ring(C, alpha, np.linspace(sync_momentum, sync_momentum, N_t + 1), Proton(), n_turns=N_t)
+rfstation = RFStation(ring, [h], [V], [dphi])
+beam = Beam(ring, N_m, N_p)
+
+# --- Insert particle ---
+beam.dt = np.array([input_dt])
+beam.dE = np.array([input_dE])
+
+# --- RF tracker object ---
+blond_tracker = RingAndRFTracker(rfstation, beam, Profile=None, interpolation=False,
+ CavityFeedback=None, BeamFeedback=None,
+ TotalInducedVoltage=None,
+ with_xsuite=True)
+
+# --- Bunch monitor ---
+# Print parameters for plotting
+print(f"\n Calculating plotting parameters...")
+
+t_u = np.pi / rfstation.omega_rf[0, 0]
+Y_phi_s = np.sqrt(abs(-np.cos(rfstation.phi_s[0]) + (np.pi - 2 * rfstation.phi_s[0]) / 2 * np.sin(rfstation.phi_s[0])))
+dE_sep = np.sqrt(2 * abs(rfstation.charge) * rfstation.voltage[0, 0] * rfstation.beta[0] * rfstation.beta[0]
+ * rfstation.energy[0] / (np.pi * rfstation.harmonic[0, 0] * abs(ring.eta_0[0, 0]))) * Y_phi_s
+
+print('t_u: ' + str(2 * t_u))
+print('dE_sep,m: ' + str(dE_sep))
+
+# Make bunchmonitor
+bunchmonitor = BunchMonitor(ring, rfstation, beam,
+ results_directory + '/EX_26_output_data', Profile=None)
+
+format_options = {'dirname': results_directory + '/EX_26_fig'}
+plots = Plot(ring, rfstation, beam, dt_plt, N_t, 0, 2 * t_u,
+ -1.05 * dE_sep, 1.05 * dE_sep, xunit='s', separatrix_plot=True,
+ Profile=None, h5file=results_directory + '/EX_26_output_data',
+ format_options=format_options)
+
+# --- Creating interface elements ---
+# Blond tracker
+cavity = BlondElement(blond_tracker, beam)
+# BLonD monitor
+beam_monitor = BlondElement(bunchmonitor, beam)
+beam_plots = BlondElement(plots, beam)
+
+# Preparing Xsuite items -----------------------------------------------------------------------------------------------
+print(f"\nPreparing Xsuite items...")
+# --- Setup matrix ---
+# Make First order matrix map (takes care of drift in Xsuite)
+matrix = xt.LineSegmentMap(
+ longitudinal_mode='nonlinear',
+ qx=1.1, qy=1.2,
+ betx=1.,
+ bety=1.,
+ voltage_rf=0,
+ frequency_rf=0,
+ lag_rf=0,
+ momentum_compaction_factor=alpha,
+ length=C)
+
+# Create line
+line = xt.Line(elements=[matrix], element_names={'matrix'})
+line['matrix'].length = C
+
+# Insert the BLonD elements
+line.insert_element(index=0, element=cavity, name='blond_cavity')
+line.append_element(element=beam_monitor, name='blond_monitor')
+line.append_element(element=beam_plots, name='blond_plots')
+
+# Add particles to line and build tracker
+line.particle_ref = xp.Particles(p0c=sync_momentum, mass0=xp.PROTON_MASS_EV, q0=1.)
+line.build_tracker()
+
+# Show table
+line.get_table().show()
+# prints:
+
+# name s element_type isthick isreplica parent_name iscollective
+# blond_cavity 0 BlondElement False False None True
+# matrix 0 LineSegmentMap True False None False
+# blond_monitor 157.08 BlondObserver False False None True
+# blond_plots 157.08 BlondElement False False None True
+# _end_point 157.08 False False None False
+
+# Simulating ----------------------------------------------------------------------------------------------------------
+print(f"\nSetting up simulation...")
+
+# --- Convert the initial BLonD distribution to xsuite coordinates ---
+zeta, ptau = blond_to_xsuite_transform(
+ beam.dt, beam.dE, line.particle_ref.beta0[0],
+ line.particle_ref.energy0[0], phi_s=rfstation.phi_s[0] - rfstation.phi_rf[0, 0], # Correct for RF phase
+ omega_rf=rfstation.omega_rf[0, 0]
+)
+
+# --- Track matrix ---
+particles = line.build_particles(x=0, y=0, px=0, py=0, zeta=np.copy(zeta), ptau=np.copy(ptau))
+line.track(particles, num_turns=N_t, turn_by_turn_monitor=True, with_progress=True)
+mon = line.record_last_track
+
+# Saving turn-by-turn particle coordinates
+test_string = ''
+test_string += '{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\n'.format(
+ 'dE', 'dt', 'momentum', 'gamma', 'beta', 'energy', 'x', 'px', 'y', 'py')
+test_string += ('{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}'
+ '\t{:+10.10e}\t{:+10.10e}\n').format(beam.dE[0], beam.dt[0], rfstation.momentum[0],
+ rfstation.gamma[0], rfstation.beta[0], rfstation.energy[0],
+ mon.x[0, 0].T, mon.px[0, 0].T, mon.y[0, 0].T, mon.py[0, 0].T)
+
+# Saving turn-by-turn Xsuite particle coordinates
+xsuite_string = ''
+xsuite_string += '{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\n'.format(
+ 'ptau', 'zeta', 'momentum', 'gamma', 'beta')
+xsuite_string += '{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\n'.format(
+ mon.ptau[0, 0].T, mon.zeta[0, 0].T, np.float64(mon.p0c[0, 0]), np.float64(mon.gamma0[0, 0]),
+ np.float64(mon.beta0[0, 0])
+)
+
+# Convert the xsuite particle coordinates back to BLonD
+for i in range(N_t):
+ dt, dE = xsuite_to_blond_transform(
+ mon.zeta[:, i].T, mon.ptau[:, i].T,
+ rfstation.beta[i], rfstation.energy[i],
+ phi_s=rfstation.phi_s[i] - rfstation.phi_rf[0, 0], # Correct for RF phase
+ omega_rf=rfstation.omega_rf[0, i]
+ )
+
+ # Statistics
+ test_string += ('{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}'
+ '\t{:+10.10e}\t{:+10.10e}\n').format(dE[0], dt[0], rfstation.momentum[i], rfstation.gamma[i],
+ rfstation.beta[i], rfstation.energy[i], mon.x[0, i].T,
+ mon.px[0, i].T, mon.y[0, i].T, mon.py[0, i].T)
+
+ xsuite_string += '{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\n'.format(
+ mon.ptau[0, i].T, mon.zeta[0, i].T, np.float64(mon.p0c[0, i]), np.float64(mon.gamma0[0, i]),
+ np.float64(mon.beta0[0, i])
+ )
+
+with open(results_directory + '/EX_26_output_data.txt', 'w') as f:
+ f.write(test_string)
+with open(results_directory + '/EX_26_output_data_xsuite.txt', 'w') as f:
+ f.write(xsuite_string)
+
+# Results --------------------------------------------------------------------------------------------------------------
+print(f"\nPlotting result...")
+
+# Plot Phasespace
+
+plt.figure()
+plt.plot(mon.zeta.T,mon.ptau.T)
+plt.scatter(np.mean(mon.zeta.T), np.mean(mon.ptau.T), label='mon')
+plt.grid()
+plt.title('Phasespace')
+plt.xlabel(r'$\zeta$')
+plt.ylabel(r'$p_{\tau}$')
+plt.savefig(results_directory + '/figure_Phasespace.png')
+
+# Use Xplt to plot Phasespace
+
+plot = xplt.PhaseSpacePlot(
+ mon
+)
+plot.fig.suptitle("Particle distribution for a single turn")
+plot.fig.savefig(results_directory + '/figure_ps_mon.png')
+
+# Show plots
+plt.show()
diff --git a/__EXAMPLES/main_files/EX_27_single_particle_ramp_psb.py b/__EXAMPLES/main_files/EX_27_single_particle_ramp_psb.py
new file mode 100644
index 0000000000000000000000000000000000000000..dd82280ebed2cf979bbd7dda6af67df4d44312f8
--- /dev/null
+++ b/__EXAMPLES/main_files/EX_27_single_particle_ramp_psb.py
@@ -0,0 +1,239 @@
+'''
+PSB simulation of a single particle with a ramp.
+
+:Author: **Thom Arnoldus van Rijswijk**
+'''
+
+# Imports
+import os
+import numpy as np
+from scipy.constants import c, e, m_p
+import matplotlib.pyplot as plt
+
+# Xsuite imports
+import xpart as xp
+import xtrack as xt
+import xplt
+
+# BLonD objects
+from blond.beam.beam import Beam, Proton
+from blond.input_parameters.ring import Ring
+from blond.input_parameters.rf_parameters import RFStation
+from blond.trackers.tracker import RingAndRFTracker
+
+# Interface
+from blond.interfaces.xsuite import (BlondElement, EnergyUpdate,
+ blond_to_xsuite_transform,
+ xsuite_to_blond_transform)
+
+# Monitor objects
+from blond.monitors.monitors import BunchMonitor
+from blond.plots.plot import Plot
+
+# Directory to save files
+this_directory = os.path.dirname(os.path.realpath(__file__)) + '/'
+results_directory = this_directory + '../output_files/EX_27'
+os.makedirs(results_directory, exist_ok=True)
+
+# Parameters ----------------------------------------------------------------------------------------------------------
+# Accelerator parameters
+C = 2 * np.pi * 25.000057945260654 # Machine circumference, radius 25 [m]
+gamma_t = 4.11635447373496 # Transition gamma [-]
+alpha = 1./gamma_t/gamma_t # First order mom. comp. factor [-]
+h = 1 # Harmonic number [-]
+V = 8e3 # RF voltage [V]
+dphi = np.pi # Phase modulation/offset [rad]
+
+# Beam parameters
+N_m = 1 # Number of macroparticles [-]
+N_p = 1e11 # Intensity / number of particles
+sigma_dt = 180e-9 / 4 # [s]
+blen = sigma_dt * 4 # Initial bunch length, 4 sigma [s]
+kin_beam_energy = 1.4e9 # Kinetic energy [eV]
+
+
+# Simulation parameters
+N_t = 6000 # Number of (tracked) turns [-]
+N_buckets = 1 # Number of buckets [-]
+dt_plt = 500 # Timestep between plots [-]
+input_dt = 5.720357923415153e-07 - sigma_dt # Input particles dt [s]
+input_dE = 0.0 # Input particles dE [eV]
+
+# Derived parameters
+E_0 = m_p * c**2 / e # [eV]
+tot_beam_energy = E_0 + kin_beam_energy # [eV]
+sync_momentum = np.sqrt(tot_beam_energy**2 - E_0**2) # [eV / c]
+momentum_increase = 0.001e9 # [eV / c]
+
+# Preparing BLonD items -----------------------------------------------------------------------------------------------
+print(f"\nPreparing BLonD items...")
+# --- BLonD objects ---
+ring = Ring(C, alpha, np.linspace(sync_momentum, sync_momentum + momentum_increase, N_t + 1), Proton(), n_turns=N_t)
+rfstation = RFStation(ring, [h], [V], [dphi])
+beam = Beam(ring, N_m, N_p)
+
+# --- Insert particle ---
+beam.dt = np.array([input_dt])
+beam.dE = np.array([input_dE])
+
+# --- RF tracker object ---
+blond_tracker = RingAndRFTracker(rfstation, beam, Profile=None, interpolation=False,
+ CavityFeedback=None, BeamFeedback=None,
+ TotalInducedVoltage=None,
+ with_xsuite=True)
+
+# --- Bunch monitor ---
+# Print parameters for plotting
+print(f"\n Calculating plotting parameters...")
+
+t_u = np.pi / rfstation.omega_rf[0, 0]
+Y_phi_s = np.sqrt(abs(-np.cos(rfstation.phi_s[0]) + (np.pi - 2 * rfstation.phi_s[0]) / 2 * np.sin(rfstation.phi_s[0])))
+dE_sep = np.sqrt(2 * abs(rfstation.charge) * rfstation.voltage[0, 0] * rfstation.beta[0] * rfstation.beta[0]
+ * rfstation.energy[0] / (np.pi * rfstation.harmonic[0, 0] * abs(ring.eta_0[0, 0]))) * Y_phi_s
+
+print('t_u: ' + str(2 * t_u))
+print('dE_sep,m: ' + str(dE_sep))
+
+# Make bunchmonitor
+bunchmonitor = BunchMonitor(ring, rfstation, beam,
+ results_directory + '/EX_27_output_data', Profile=None)
+
+format_options = {'dirname': results_directory + '/EX_27_fig'}
+plots = Plot(ring, rfstation, beam, dt_plt, N_t, 0, 2 * t_u,
+ -1.05 * dE_sep, 1.05 * dE_sep, xunit='s', separatrix_plot=True,
+ Profile=None, h5file=results_directory + '/EX_27_output_data',
+ format_options=format_options)
+
+# --- Creating interface elements ---
+# Blond tracker
+cavity = BlondElement(blond_tracker, beam)
+# BLonD monitor
+beam_monitor = BlondElement(bunchmonitor, beam)
+beam_plots = BlondElement(plots, beam)
+
+# Preparing Xsuite items -----------------------------------------------------------------------------------------------
+print(f"\nPreparing Xsuite items...")
+# --- Setup matrix ---
+# Make First order matrix map (takes care of drift in Xsuite)
+matrix = xt.LineSegmentMap(
+ longitudinal_mode='nonlinear',
+ qx=1.1, qy=1.2,
+ betx=1.,
+ bety=1.,
+ voltage_rf=0,
+ frequency_rf=0,
+ lag_rf=0,
+ momentum_compaction_factor=alpha,
+ length=C)
+
+# Create line
+line = xt.Line(elements=[matrix], element_names={'matrix'})
+line['matrix'].length = C
+
+# Insert the BLonD elements
+line.insert_element(index=0, element=cavity, name='blond_cavity')
+line.append_element(element=beam_monitor, name='blond_monitor')
+line.append_element(element=beam_plots, name='blond_plots')
+
+# Insert energy ramp
+energy_update = EnergyUpdate(ring.momentum[0, :])
+line.insert_element(index=0, element=energy_update, name='energy_update')
+
+# Add particles to line and build tracker
+line.particle_ref = xp.Particles(p0c=sync_momentum, mass0=xp.PROTON_MASS_EV, q0=1.)
+line.build_tracker()
+
+# Show table
+line.get_table().show()
+# prints:
+
+# name s element_type isthick isreplica parent_name iscollective
+# energy_update 0 EnergyUpdate False False None True
+# blond_cavity 0 BlondElement False False None True
+# matrix 0 LineSegmentMap True False None False
+# blond_monitor 157.08 BlondObserver False False None True
+# blond_plots 157.08 BlondElement False False None True
+# _end_point 157.08 False False None False
+
+# Simulating ----------------------------------------------------------------------------------------------------------
+print(f"\nSetting up simulation...")
+
+# --- Convert the initial BLonD distribution to xsuite coordinates ---
+zeta, ptau = blond_to_xsuite_transform(
+ beam.dt, beam.dE, line.particle_ref.beta0[0],
+ line.particle_ref.energy0[0], phi_s=rfstation.phi_s[0] - rfstation.phi_rf[0, 0], # Correct for RF phase
+ omega_rf=rfstation.omega_rf[0, 0]
+)
+
+# --- Track matrix ---
+particles = line.build_particles(x=0, y=0, px=0, py=0, zeta=np.copy(zeta), ptau=np.copy(ptau))
+line.track(particles, num_turns=N_t, turn_by_turn_monitor=True, with_progress=True)
+mon = line.record_last_track
+
+# Saving turn-by-turn particle coordinates
+test_string = ''
+test_string += '{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\n'.format(
+ 'dE', 'dt', 'momentum', 'gamma', 'beta', 'energy', 'x', 'px', 'y', 'py')
+test_string += ('{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}'
+ '\t{:+10.10e}\t{:+10.10e}\n').format(beam.dE[0], beam.dt[0], rfstation.momentum[0],
+ rfstation.gamma[0], rfstation.beta[0], rfstation.energy[0],
+ mon.x[0, 0].T, mon.px[0, 0].T, mon.y[0, 0].T, mon.py[0, 0].T)
+
+# Saving turn-by-turn Xsuite particle coordinates
+xsuite_string = ''
+xsuite_string += '{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\n'.format(
+ 'ptau', 'zeta', 'momentum', 'gamma', 'beta')
+xsuite_string += '{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\n'.format(
+ mon.ptau[0, 0].T, mon.zeta[0, 0].T, np.float64(mon.p0c[0, 0]), np.float64(mon.gamma0[0, 0]),
+ np.float64(mon.beta0[0, 0])
+)
+
+# Convert the xsuite particle coordinates back to BLonD
+for i in range(N_t):
+ dt, dE = xsuite_to_blond_transform(
+ mon.zeta[:, i].T, mon.ptau[:, i].T,
+ rfstation.beta[i], rfstation.energy[i],
+ phi_s=rfstation.phi_s[i] - rfstation.phi_rf[0, 0], # Correct for RF phase
+ omega_rf=rfstation.omega_rf[0, i]
+ )
+
+ # Statistics
+ test_string += ('{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}'
+ '\t{:+10.10e}\t{:+10.10e}\n').format(dE[0], dt[0], rfstation.momentum[i], rfstation.gamma[i],
+ rfstation.beta[i], rfstation.energy[i], mon.x[0, i].T,
+ mon.px[0, i].T, mon.y[0, i].T, mon.py[0, i].T)
+
+ xsuite_string += '{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\n'.format(
+ mon.ptau[0, i].T, mon.zeta[0, i].T, np.float64(mon.p0c[0, i]), np.float64(mon.gamma0[0, i]),
+ np.float64(mon.beta0[0, i])
+ )
+
+with open(results_directory + '/EX_27_output_data.txt', 'w') as f:
+ f.write(test_string)
+with open(results_directory + '/EX_27_output_data_xsuite.txt', 'w') as f:
+ f.write(xsuite_string)
+
+# Results --------------------------------------------------------------------------------------------------------------
+print(f"\nPlotting result...")
+
+# Plot Phasespace
+
+plt.figure()
+plt.plot(mon.zeta.T,mon.ptau.T)
+plt.scatter(np.mean(mon.zeta.T), np.mean(mon.ptau.T), label='mon')
+plt.grid()
+plt.title('Phasespace')
+plt.xlabel(r'$\zeta$')
+plt.ylabel(r'$p_{\tau}$')
+plt.savefig(results_directory + '/figure_Phasespace.png')
+
+# Use Xplt to plot Phasespace
+
+plot = xplt.PhaseSpacePlot(
+ mon
+)
+plot.fig.suptitle("Particle distribution for a single turn")
+plot.fig.savefig(results_directory + '/figure_ps_mon.png')
+
+# Show plots
+plt.show()
diff --git a/__EXAMPLES/main_files/EX_28_phase_loop_lhc.py b/__EXAMPLES/main_files/EX_28_phase_loop_lhc.py
new file mode 100644
index 0000000000000000000000000000000000000000..bcb349ff863f0551d6386602af3dc35c3e1043b3
--- /dev/null
+++ b/__EXAMPLES/main_files/EX_28_phase_loop_lhc.py
@@ -0,0 +1,247 @@
+'''
+LHC simulation using the LHC phase loop.
+
+:Author: **Birk Emil Karlsen-Baeck**, **Thom Arnoldus van Rijswijk**
+'''
+
+# Imports
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Xsuite imports
+import xpart as xp
+import xtrack as xt
+import xplt
+
+# BLonD objects
+from blond.beam.beam import Beam, Proton
+from blond.beam.profile import Profile, CutOptions, FitOptions
+from blond.beam.distributions import parabolic
+from blond.input_parameters.ring import Ring
+from blond.input_parameters.rf_parameters import RFStation
+from blond.trackers.tracker import RingAndRFTracker
+from blond.llrf.beam_feedback import BeamFeedback
+
+# Interface
+from blond.interfaces.xsuite import (BlondElement, BlondObserver,
+ blond_to_xsuite_transform,
+ xsuite_to_blond_transform)
+
+# Monitor objects
+from blond.monitors.monitors import BunchMonitor
+from blond.plots.plot import Plot
+
+# Directory to save files
+this_directory = os.path.dirname(os.path.realpath(__file__)) + '/'
+results_directory = this_directory + '../output_files/EX_28'
+os.makedirs(results_directory, exist_ok=True)
+
+# Parameters ----------------------------------------------------------------------------------------------------------
+# Accelerator parameters
+C = 26658.8832 # Machine circumference [m]
+p_s = 450e9 # Synchronous momentum [eV/c]
+p_f = 450.0e9 # Synchronous momentum, final
+h = 35640 # Harmonic number [-]
+alpha = 0.00034849575112251314 # First order mom. comp. factor [-]
+V = 5e6 # RF voltage [V]
+dphi = 0 # Phase modulation/offset [rad]
+
+# Bunch parameters
+N_bunches = 1 # Number of bunches [-]
+N_m = int(1e3) # Number of macroparticles [-]
+N_p = 1.15e11 # Intensity
+blen = 1.25e-9 # Bunch length [s]
+energy_err = 100e6 # Beamenergy error [eV/c]
+
+# Simulation parameters
+N_t = 2 * 481 # Number of (tracked) turns [-]
+N_buckets = 2 # Number of buckets [-]
+dt_plt = 50 # Timestep between plots [-]
+
+# Preparing BLonD items -----------------------------------------------------------------------------------------------
+print(f"\nPreparing BLonD items...")
+# --- BLonD objects ---
+ring = Ring(C, alpha, np.linspace(p_s, p_f, N_t + 1), Proton(), n_turns=N_t)
+rfstation = RFStation(ring, [h], [V], [dphi])
+beam = Beam(ring, N_m, N_p)
+
+# --- Insert particles ---
+parabolic(ring, rfstation, beam, bunch_length=blen)
+beam.dE += energy_err
+
+# Create Profile
+profile = Profile(beam,
+ CutOptions(cut_left=-1.5 * rfstation.t_rf[0, 0],
+ cut_right=(N_buckets + 1.5) * rfstation.t_rf[0, 0],
+ n_slices=(N_buckets + 3) * 2**5),
+ FitOptions(fit_option='gaussian'))
+profile.track()
+
+# --- LHC beam-phase loop and synchronization loop ---
+PL_gain = 1 / (5 * ring.t_rev[0])
+SL_gain = PL_gain / 10
+
+bl_config = {'machine': 'LHC',
+ 'PL_gain': PL_gain,
+ 'SL_gain': SL_gain}
+BL = BeamFeedback(ring, rfstation, profile, bl_config)
+
+# --- RF tracker object ---
+blond_tracker = RingAndRFTracker(rfstation, beam, Profile=profile, interpolation=False,
+ CavityFeedback=None, BeamFeedback=BL,
+ TotalInducedVoltage=None,
+ with_xsuite=True)
+
+# --- Bunch monitor ---
+# Print parameters for plotting
+print(f"\n Calculating plotting parameters...")
+
+t_u = np.pi / rfstation.omega_rf[0, 0]
+Y_phi_s = np.sqrt(abs(-np.cos(rfstation.phi_s[0]) + (np.pi - 2 * rfstation.phi_s[0]) / 2 * np.sin(rfstation.phi_s[0])))
+dE_sep = np.sqrt(2 * abs(rfstation.charge) * rfstation.voltage[0, 0] * rfstation.beta[0] * rfstation.beta[0]
+ * rfstation.energy[0] / (np.pi * rfstation.harmonic[0, 0] * abs(ring.eta_0[0, 0]))) * Y_phi_s
+
+print('t_u: ' + str(2 * t_u))
+print('dE_sep,m: ' + str(dE_sep))
+
+# Make bunchmonitor
+bunchmonitor = BunchMonitor(ring, rfstation, beam,
+ results_directory + '/EX_28_output_data', Profile=profile, PhaseLoop=BL)
+
+format_options = {'dirname': results_directory + '/EX_28_fig'}
+plots = Plot(ring, rfstation, beam, dt_plt, N_t, 0, 2*t_u,
+ -1.05 * dE_sep, 1.05 * dE_sep, xunit='s', separatrix_plot=True,
+ Profile=profile, PhaseLoop=BL, h5file=results_directory + '/EX_28_output_data',
+ format_options=format_options)
+
+# --- Creating interface elements ---
+# Blond tracker
+cavity = BlondElement(blond_tracker, beam)
+# BLonD monitor
+beam_monitor = BlondElement(bunchmonitor, beam)
+beam_profile = BlondElement(profile, beam)
+beam_plots = BlondElement(plots, beam)
+
+# Preparing Xsuite items -----------------------------------------------------------------------------------------------
+print(f"\nPreparing Xsuite items...")
+# --- Setup matrix ---
+# Make First order matrix map (takes care of drift in Xsuite)
+matrix = xt.LineSegmentMap(
+ longitudinal_mode='nonlinear',
+ qx=1.1, qy=1.2,
+ betx=1.,
+ bety=1.,
+ voltage_rf=0,
+ frequency_rf=0,
+ lag_rf=0,
+ momentum_compaction_factor=alpha,
+ length=C)
+
+# Create line
+line = xt.Line(elements=[matrix], element_names={'matrix'})
+line['matrix'].length = C
+
+# Insert the BLonD elements
+line.insert_element(index=0, element=cavity, name='blond_cavity')
+line.append_element(element=beam_monitor, name='blond_monitor')
+line.append_element(element=beam_profile, name='blond_profile')
+line.append_element(element=beam_plots, name='blond_plots')
+
+# Add particles to line and build tracker
+line.particle_ref = xp.Particles(p0c=p_s, mass0=xp.PROTON_MASS_EV, q0=1.)
+line.build_tracker()
+
+# Show table
+line.get_table().show()
+# prints:
+
+# name s element_type isthick isreplica parent_name iscollective
+# blond_cavity 0 BlondElement False False None True
+# matrix 0 LineSegmentMap True False None False
+# blond_monitor 26658.9 BlondObserver False False None True
+# blond_plots 26658.9 BlondElement False False None True
+# _end_point 26658.9 False False None False
+
+# Simulating ----------------------------------------------------------------------------------------------------------
+print(f"\nSetting up simulation...")
+
+# --- Convert the initial BLonD distribution to xsuite coordinates ---
+zeta, ptau = blond_to_xsuite_transform(
+ beam.dt, beam.dE, line.particle_ref.beta0[0],
+ line.particle_ref.energy0[0], phi_s=rfstation.phi_s[0],
+ omega_rf=rfstation.omega_rf[0, 0]
+)
+
+# --- Track matrix ---
+particles = line.build_particles(x=0, y=0, px=0, py=0, zeta=np.copy(zeta), ptau=np.copy(ptau))
+line.track(particles, num_turns=N_t, turn_by_turn_monitor=True, with_progress=True)
+mon = line.record_last_track
+
+# Saving turn-by-turn particle coordinates
+test_string = ''
+test_string += '{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\n'.format(
+ 'dE', 'dt', 'momentum', 'gamma', 'beta', 'energy', 'x', 'px', 'y', 'py')
+test_string += ('{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}'
+ '\t{:+10.10e}\t{:+10.10e}\n').format(beam.dE[0], beam.dt[0], rfstation.momentum[0],
+ rfstation.gamma[0], rfstation.beta[0], rfstation.energy[0],
+ mon.x[0, 0].T, mon.px[0, 0].T, mon.y[0, 0].T, mon.py[0, 0].T)
+
+# Saving turn-by-turn Xsuite particle coordinates
+xsuite_string = ''
+xsuite_string += '{:<17}\t{:<17}\t{:<17}\t{:<17}\t{:<17}\n'.format(
+ 'ptau', 'zeta', 'momentum', 'gamma', 'beta')
+xsuite_string += '{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\n'.format(
+ mon.ptau[0, 0].T, mon.zeta[0, 0].T, np.float64(mon.p0c[0, 0]), np.float64(mon.gamma0[0, 0]),
+ np.float64(mon.beta0[0, 0])
+)
+
+# Convert the xsuite particle coordinates back to BLonD
+for i in range(N_t):
+ dt, dE = xsuite_to_blond_transform(
+ mon.zeta[:, i].T, mon.ptau[:, i].T,
+ rfstation.beta[i], rfstation.energy[i],
+ phi_s=rfstation.phi_s[i] - rfstation.phi_rf[0, 0], # Correct for RF phase
+ omega_rf=rfstation.omega_rf[0, i]
+ )
+
+ # Statistics
+ test_string += ('{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}'
+ '\t{:+10.10e}\t{:+10.10e}\n').format(dE[0], dt[0], rfstation.momentum[i], rfstation.gamma[i],
+ rfstation.beta[i], rfstation.energy[i], mon.x[0, i].T,
+ mon.px[0, i].T, mon.y[0, i].T, mon.py[0, i].T)
+
+ xsuite_string += '{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\t{:+10.10e}\n'.format(
+ mon.ptau[0, i].T, mon.zeta[0, i].T, np.float64(mon.p0c[0, i]), np.float64(mon.gamma0[0, i]),
+ np.float64(mon.beta0[0, i])
+ )
+
+with open(results_directory + '/EX_28_output_data.txt', 'w') as f:
+ f.write(test_string)
+with open(results_directory + '/EX_28_output_data_xsuite.txt', 'w') as f:
+ f.write(xsuite_string)
+
+# Results --------------------------------------------------------------------------------------------------------------
+print(f"\nPlotting result...")
+
+# Plot Phasespace
+
+plt.figure()
+plt.plot(mon.zeta.T,mon.ptau.T)
+plt.scatter(np.mean(mon.zeta.T), np.mean(mon.ptau.T), label='mon')
+plt.grid()
+plt.title('Phasespace')
+plt.xlabel(r'$\zeta$')
+plt.ylabel(r'$p_{\tau}$')
+plt.savefig(results_directory + '/figure_Phasespace.png')
+
+# Use Xplt to plot Phasespace
+
+plot = xplt.PhaseSpacePlot(
+ mon
+)
+plot.fig.suptitle("Particle distribution for a single turn")
+plot.fig.savefig(results_directory + '/figure_ps_mon.png')
+
+# Show plots
+plt.show()
diff --git a/blond/interfaces/__init__.py b/blond/interfaces/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/blond/interfaces/xsuite.py b/blond/interfaces/xsuite.py
new file mode 100644
index 0000000000000000000000000000000000000000..bd30be38886831b8d2597239353660936382ec73
--- /dev/null
+++ b/blond/interfaces/xsuite.py
@@ -0,0 +1,404 @@
+'''
+Functions and classes to interface BLonD with xsuite.
+
+:Authors: **Birk Emil Karlsen-Baeck**, **Thom Arnoldus van Rijswijk**, **Helga Timko**
+'''
+
+from __future__ import annotations
+from typing import TYPE_CHECKING
+import numpy as np
+from scipy.constants import c as clight
+
+from xtrack import ReferenceEnergyIncrease, ZetaShift
+from blond.trackers.tracker import RingAndRFTracker
+from blond.impedances.impedance import TotalInducedVoltage
+
+if TYPE_CHECKING:
+ from xtrack import Particles
+ from xtrack import Line
+ from blond.beam.beam import Beam
+ from blond.beam.profile import Profile
+ from typing import Any, Union
+
+
+def blond_to_xsuite_transform(dt: Union[float, np.ndarray], de: Union[float, np.ndarray], beta0: float,
+ energy0: float, omega_rf: float, phi_s: float = 0):
+ r'''
+ Coordinate transformation from Xsuite to BLonD at a given turn or multiple turns if numpy arrays ar given.
+
+ :param dt: float or NDArray.
+ The deviation in time from the reference clock in BLonD.
+ :param de: float or NDArray.
+ The deviation in energy from the synchronous particle.
+ :param beta0: float.
+ Synchronous beta.
+ :param energy0: float.
+ Synchronous energy.
+ :param omega_rf: float.
+ The rf angular frequency.
+ :param phi_s: float.
+ Synchronous phase in radians equivalent to Xsuite's phi_rf (below transition energy input should be phi_s-phi_rf).
+ The default value is 0.
+
+ :return zeta, ptau: as numpy-arrays (or single variable)
+ '''
+
+ ptau = de / (beta0 * energy0)
+ zeta = -(dt - phi_s / omega_rf) * beta0 * clight
+ return zeta, ptau
+
+
+def xsuite_to_blond_transform(zeta: Union[float, np.ndarray], ptau: Union[float, np.ndarray],
+ beta0: float, energy0: float, omega_rf: float, phi_s: float = 0):
+ r'''
+ Coordinate transformation from Xsuite to BLonD.
+
+ :param zeta: float or numpy-array.
+ The zeta coordinate as defined in Xsuite.
+ :param ptau: float or numpy-array.
+ The ptau coordinate as defined in Xsuite.
+ :param beta0: float.
+ The synchronous beta.
+ :param energy0: float.
+ The synchronous energy in eV.
+ :param omega_rf: float.
+ The rf angular frequency.
+ :param phi_s: float.
+ The synchronous phase in radians equivalent to Xsuite's phi_rf
+ (below transition energy input should be phi_s-phi_rf)
+
+ :return dt, dE: as numpy-arrays (or single variable)
+ '''
+
+ dE = ptau * beta0 * energy0
+ dt = -zeta / (beta0 * clight) + phi_s / omega_rf
+ return dt, dE
+
+
+class BlondElement:
+
+ def __init__(self, trackable: Any, beam: Beam, update_zeta: bool = False) -> None:
+ r"""
+ The BlondElement class contains a trackable object from the BLonD simulation suite and the Beam object from
+ BLonD. The class is intended to be an element to be added to the Line in XTrack instead of the default
+ RF cavities in XTrack.
+
+ The behavior of the tracking of the class depends on what trackable object is passed.
+ If the RingAndRFTracker class is passed then a coordinate transformation will be performed and the energy
+ deviation of the particles will be updated. If the TotalInducedVoltage class is passed then the
+ TotalInducedVoltage.induced_voltage_sum() method will be called. Lastly, if any other class is passed then
+ the tracking of BlondElement will track this class.
+
+ trackable : BLonD trackable object
+ BLonD object to be tracked, e.g. the RingAndRFTracker.
+ beam : blond.beam.beam.Beam class
+ BLonD Beam class used for tracking.
+ """
+
+ self.trackable = trackable
+ self.beam = beam
+ self.update_zeta = update_zeta
+
+ # Initialize time- and orbit-shift to BLonD coordinates
+ self._dt_shift = None
+ self.orbit_shift = ZetaShift(dzeta=0)
+
+ # Check what BLonD trackable has been passed and track the object accordingly
+ if isinstance(self.trackable, RingAndRFTracker):
+ self.track = getattr(self, "rf_track")
+ elif isinstance(self.trackable, TotalInducedVoltage):
+ self.track = getattr(self, "ind_track")
+ else:
+ self.track = getattr(self, "obs_track")
+
+ def rf_track(self, particles: Particles) -> None:
+ r"""
+ Tracking method which is called if the trackable BLonD class is the RingAndRFTracker.
+
+ particles : xtrack.Particles
+ Particles class from xtrack
+ """
+
+ # Compute the shift to BLonD coordinates
+ self._get_time_shift()
+
+ # Convert the Xsuite coordinates to BLonD coordinates
+ self.xsuite_part_to_blond_beam(particles)
+
+ # Track with only the energy kick
+ self.trackable.track()
+
+ # Convert the BLonD energy coordinate to the equivalent Xsuite coordinate
+ self.blond_beam_to_xsuite_part(particles, self.update_zeta)
+
+ # Update the zeta shift due to potential frequency shifts in the RF
+ self._orbit_shift(particles)
+
+ def ind_track(self, particles: Particles) -> None:
+ r"""
+ Tracking method which is called if the trackable BLonD class is the TotalInducedVoltage.
+
+ particles : xtrack.Particles
+ Particles class from xtrack
+ """
+ self.trackable.induced_voltage_sum()
+
+ def obs_track(self, particles: Particles) -> None:
+ r"""
+ Tracking method which is called if the trackable BLonD class is not the RingAndRFTracker or TotalInducedVoltage.
+
+ particles : xtrack.Particles
+ Particles class from xtrack
+ """
+ self.trackable.track()
+
+ def xsuite_part_to_blond_beam(self, particles: Particles) -> None:
+ r"""
+ Coordinate transformation from Xsuite to BLonD.
+ It uses the initial particle coordinates and beam properties in stored in the Particles class from xtrack.
+
+ particles : xtrack.Particles
+ Particles class from xtrack
+ """
+ # Convert Xsuite momentum to BLonD energy deviation
+ self.beam.dE[:] = particles.beta0 * particles.energy0 * particles.ptau
+
+ # Convert Xsuite zeta coordinate to BLonD time deviation
+ self.beam.dt[:] = -particles.zeta / particles.beta0 / clight + self._dt_shift
+
+ # Check what particles are still alive
+ self.beam.id[:] = np.int_(particles.state > 0)
+
+ def blond_beam_to_xsuite_part(self, particles: Particles, update_zeta: bool = False) -> None:
+ r"""
+ Coordinate transformation from BLonD to Xsuite.
+ It uses the particle coordinates stored in Beam class of BLonD
+ It uses the beam properties in stored in the Particles class from xtrack
+
+ particles : xtrack.Particles
+ Particles class from xtrack
+ """
+ # Subtract the given acceleration kick in BLonD, in Xsuite this is dealt with differently
+ if isinstance(self.trackable, RingAndRFTracker):
+ self.beam.dE = self.beam.dE - self.trackable.acceleration_kick[self.trackable.counter[0] - 1]
+
+ # Convert BLonD energy deviation to Xsuite momentum
+ particles.ptau = self.beam.dE / (particles.beta0 * particles.energy0)
+
+ # Convert BLonD time deviation to Xsuite zeta.
+ # This step is not needed usually because the BLonD simulation only does the kick, so dt is not changed.
+ if update_zeta:
+ particles.zeta = - (self.beam.dt - self._dt_shift) * particles.beta0 * clight
+
+ # Check what particles are still alive after the BLonD track
+ mask_lost = (self.beam.id <= 0) & particles.state > 0
+
+ # If the particle is lost its state is set to -500 by convention
+ particles.state[mask_lost] = -500
+
+ def _get_time_shift(self) -> None:
+ r'''
+ Computes the time-shift between the Xsuite and BLonD coordinate systems.
+ '''
+ # Get turn counter from the RingAndRFTracker
+ counter = self.trackable.rf_params.counter[0]
+
+ # Compute the time-shift based on the synchronous phase
+ self._dt_shift = ((self.trackable.rf_params.phi_s[counter] - self.trackable.rf_params.phi_rf[0, counter])
+ / self.trackable.rf_params.omega_rf[0, counter])
+
+ def _orbit_shift(self, particles) -> None:
+ r'''
+ Computes the radial steering due to rf periods which are not an integer multiple of the revolution period.
+ This is for example needed when tracking with global LLRF feedback loops.
+
+ particles : xtrack.Particles
+ Particles class from xtrack
+ '''
+ # Get turn counter from the RingAndRFTracker
+ counter = self.trackable.counter[0]
+
+ # Compute the orbit shift due to the difference in rf frequency
+ dzeta = self.trackable.rf_params.ring_circumference
+ omega_rf = self.trackable.rf_params.omega_rf[:, counter]
+ omega_rf_design = (2 * np.pi * self.trackable.rf_params.harmonic[:, counter]
+ / self.trackable.rf_params.t_rev[counter])
+ domega = omega_rf - omega_rf_design
+
+ dzeta *= domega / omega_rf_design
+
+ # Apply the shift
+ self.orbit_shift.dzeta = dzeta
+
+ # Track the shift
+ self.orbit_shift.track(particles)
+
+
+class EnergyUpdate:
+
+ def __init__(self, momentum: Union[list, np.ndarray]) -> None:
+ r"""
+ Class to update the synchronous energy from the momentum program in BLonD.
+
+ :param momentum: list or numpy-array.
+ Momentum program from BLonD in units of eV
+ """
+
+ # Load momentum program
+ self.momentum = momentum
+
+ # Find initial momentum update
+ init_p0c = self.momentum[1] - self.momentum[0]
+
+ # Enter the initial momentum update in the ReferenceEnergyIncrease class in xsuite
+ self.xsuite_energy_update = ReferenceEnergyIncrease(Delta_p0c=init_p0c)
+
+ def track(self, particles: Particles) -> None:
+ r'''
+ Track method for the class to update the synchronous energy.
+
+ particles : xtrack.Particles
+ Particles class from xtrack
+ '''
+ # Check for particles which are still alive
+ mask_alive = particles.state > 0
+
+ # Use the still alive particles to find the current turn momentum
+ p0c_before = particles.p0c[mask_alive]
+
+ # Find the momentum for the next turn
+ p0c_after = self.momentum[particles.at_turn[mask_alive][0]]
+
+ # Update the energy increment
+ self.xsuite_energy_update.Delta_p0c = p0c_after - p0c_before[0]
+
+ # Apply the energy increment to the particles
+ self.xsuite_energy_update.track(particles)
+
+
+class EnergyFrequencyUpdate:
+
+ def __init__(self, momentum: np.ndarray, f_rf: np.ndarray, line: Line, cavity_name: str) -> None:
+ r"""
+ Class to update energy of Particles class turn-by-turn with the ReferenceEnergyIncrease function
+ from xtrack. Additionally it updates the frequency of the xtrack cavity in the line.
+ Intended to be used without BLonD-Xsuite interface.
+
+ :param momentum: numpy-array
+ The momentum program from BLonD in eV
+ :param f_rf: numpy-array
+ The frequency program from BLonD in Hz.
+ particles : xtrack.Line
+ Line class from xtrack
+ :param cavity_name: string
+ Name of cavity to update frequency.
+ """
+
+ # Load the parameters
+ self.momentum = momentum
+ self.f_rf = f_rf
+ self.line = line
+ self.cavity_name = cavity_name
+
+ # Find initial momentum update
+ init_p0c = self.momentum[1] - self.momentum[0]
+
+ # Enter the initial momentum update in the ReferenceEnergyIncrease class in xsuite
+ self.xsuite_energy_update = ReferenceEnergyIncrease(Delta_p0c=init_p0c)
+
+ def track(self, particles: Particles) -> None:
+ r'''
+ Track-method from for the class. This method updates the synchronous momentum and the rf frequency.
+
+ particles : xtrack.Particles
+ Particles class from xtrack
+ '''
+ # Check for particles which are still alive
+ mask_alive = particles.state > 0
+
+ # Use the still alive particles to find the current turn momentum
+ p0c_before = particles.p0c[mask_alive]
+
+ # Find the momentum for the next turn
+ p0c_after = self.momentum[particles.at_turn[mask_alive][0]]
+
+ # Update the energy increment
+ self.xsuite_energy_update.Delta_p0c = p0c_after - p0c_before[0]
+
+ # Apply the energy increment to the particles
+ self.xsuite_energy_update.track(particles)
+
+ # Update the rf frequency
+ self.line[self.cavity_name].frequency = self.f_rf[particles.at_turn[mask_alive][0]]
+
+
+class BlondObserver(BlondElement):
+
+ def __init__(self, trackable: Any, beam: Beam, blond_cavity: bool, update_zeta: bool = False,
+ profile: Profile = None) -> None:
+ r'''
+ Child-class of the BlondElement, except that it updates the coordinates
+ in BLonD when an observing element is used such as BunchMonitor.
+
+ :param trackable: BLonD trackable object
+ BLonD object to be tracked, e.g. the RingAndRFTracker.
+ :param beam: blond.beam.beam.Beam class
+ BLonD Beam class used for tracking.
+ :param blond_cavity: bool.
+ If there is no BlondCavity (bool = False), it updates its own turn-counter.
+ :param update_zeta: bool.
+ Boolean that decides whether zeta is converter back to dt after tracking object or not.
+ Usually not necessary so default is False.
+ :param profile: blond.beam.profile.Profile class
+ BLonD Profile class used for tracking.
+ '''
+ # Initialize the parent class
+ super().__init__(trackable, beam, update_zeta)
+
+ # Load the parameters
+ self.blond_cavity = blond_cavity
+ # For bunch monitoring we need the profile separately added
+ self.profile = profile
+
+ # Initializing arrays for storing some turn-by-turn properties of xtrack.Particles class
+ self.xsuite_ref_energy = np.zeros(self.trackable.rf_params.n_turns + 1)
+ self.xsuite_trev = np.zeros(self.trackable.rf_params.n_turns + 1)
+
+ def obs_track(self, particles: Particles) -> None:
+ r'''
+ observation tracker which performs the coordinate transformations.
+
+ particles : xtrack.Particles
+ Particles class from xtrack
+ '''
+ # Compute the shift to BLonD coordinates
+ self._get_time_shift()
+
+ # Convert the Xsuite coordinates to BLonD coordinates
+ self.xsuite_part_to_blond_beam(particles)
+
+ # Track profile if given
+ if self.profile is not None:
+ self.profile.track()
+
+ # Track
+ self.trackable.track()
+
+ if not self.blond_cavity:
+ # Updating the beam synchronous momentum etc.
+ turn = self.trackable.rf_params.counter[0]
+
+ self.beam.beta = self.trackable.rf_params.beta[turn + 1]
+ self.beam.gamma = self.trackable.rf_params.gamma[turn + 1]
+ self.beam.energy = self.trackable.rf_params.energy[turn + 1]
+ self.beam.momentum = self.trackable.rf_params.momentum[turn + 1]
+
+ # Update custom counter
+ self.trackable.rf_params.counter[0] += 1
+
+ # Convert the BLonD energy coordinate to the equivalent Xsuite coordinate
+ self.blond_beam_to_xsuite_part(particles, self.update_zeta)
+
+ # Track properties of xtrack.Particles
+ self.xsuite_ref_energy[self.trackable.rf_params.counter[0] - 1] = particles.energy0[0]
+ self.xsuite_trev[self.trackable.rf_params.counter[0] - 1] = particles.t_sim # Does not update per turn!
diff --git a/unittests/interfaces/__init__.py b/unittests/interfaces/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/unittests/interfaces/test_xsuite.py b/unittests/interfaces/test_xsuite.py
new file mode 100644
index 0000000000000000000000000000000000000000..c5a0e57d4e8ed6721ee5a53220c14107a23c51d7
--- /dev/null
+++ b/unittests/interfaces/test_xsuite.py
@@ -0,0 +1,629 @@
+# coding: utf8
+# Copyright 2014-2017 CERN. This software is distributed under the
+# terms of the GNU General Public Licence version 3 (GPL Version 3),
+# copied verbatim in the file LICENCE.md.
+# In applying this licence, CERN does not waive the privileges and immunities
+# granted to it by virtue of its status as an Intergovernmental Organization or
+# submit itself to any jurisdiction.
+# Project website: http://blond.web.cern.ch/
+
+"""
+Unittest for interfaces.xsuite
+
+:Authors: **Birk Emil Karlsen-Bæck**, **Thom Arnoldus van Rijswijk**
+"""
+
+import os
+import unittest
+
+import numpy as np
+from scipy.constants import c, e, m_p
+
+import xpart as xp
+import xtrack as xt
+
+from blond.beam.beam import Beam, Proton
+from blond.beam.profile import CutOptions, Profile
+from blond.input_parameters.rf_parameters import RFStation
+from blond.input_parameters.ring import Ring
+from blond.trackers.tracker import RingAndRFTracker, FullRingAndRF
+from blond.monitors.monitors import BunchMonitor
+
+from blond.interfaces.xsuite import (BlondElement, BlondObserver,
+ EnergyUpdate,
+ blond_to_xsuite_transform,
+ xsuite_to_blond_transform)
+
+# TODO: Make a test of the EnergyFrequencyUpdate and BlondObserver classes
+
+
+class TestXsuiteBLonDTransforms(unittest.TestCase):
+
+ def setUp(self):
+ # Accelerator parameters
+ C = 26658.8832
+ p_s = 450e9
+ h = 35640
+ alpha = 0.00034849575112251314
+ V = 5e6
+ dphi = 0
+
+ # BLonD ring and RF station
+ self.ring = Ring(C, alpha, p_s, Proton(), n_turns=1)
+ self.rfstation = RFStation(self.ring, [h], [V], [dphi])
+
+ def testBlondToXsuiteTransform(self):
+ # Arrays of time and energy offsets
+ dts = np.linspace(
+ self.rfstation.t_rf[0, 0] * 1e-3,
+ self.rfstation.t_rf[0, 0], 10
+ )
+ dEs = np.linspace(-200e6, 200e6, 10)
+
+ # Transform BLonD coordinates to Xsuite
+ zeta, ptau = blond_to_xsuite_transform(
+ dts, dEs, self.ring.beta[0, 0], self.ring.energy[0, 0],
+ self.rfstation.omega_rf[0, 0], self.rfstation.phi_s[0]
+ )
+
+ # Correct values
+ act_zeta = np.array(
+ [0.37325428, 0.29022578, 0.20719727, 0.12416876, 0.04114025,
+ -0.04188826, -0.12491676, -0.20794527, -0.29097378, -0.37400229]
+ )
+ act_ptau = np.array(
+ [-4.44444444e-04, -3.45679012e-04, -2.46913580e-04, -1.48148148e-04, -4.93827160e-05,
+ 4.93827160e-05, 1.48148148e-04, 2.46913580e-04, 3.45679012e-04, 4.44444444e-04]
+ )
+
+ # Check if the coordinates are the same
+ np.testing.assert_allclose(
+ act_zeta, zeta,
+ err_msg='In testBlondToXsuiteTransform the zeta-coordinate differs'
+ )
+ np.testing.assert_allclose(
+ act_ptau, ptau,
+ err_msg='In testBlondToXsuiteTransform the ptau-coordinate differs'
+ )
+
+ def testXsuiteToBlondTransform(self):
+ # Arrays of z and ptau from xsuite
+ zeta = np.linspace(
+ -0.37, 0.37, 10
+ )
+ ptau = np.linspace(
+ -4e-4, 4e-4, 10
+ )
+
+ # Transform BLonD coordinates to Xsuite
+ dt, de = xsuite_to_blond_transform(
+ zeta, ptau, self.ring.beta[0, 0], self.ring.energy[0, 0],
+ self.rfstation.omega_rf[0, 0], self.rfstation.phi_s[0]
+ )
+
+ # Correct values
+ act_dt = np.array(
+ [2.48172990e-09, 2.20746549e-09, 1.93320108e-09, 1.65893668e-09, 1.38467227e-09,
+ 1.11040786e-09, 8.36143453e-10, 5.61879046e-10, 2.87614638e-10, 1.33502300e-11]
+ )
+ act_de = np.array(
+ [-1.8e+08, -1.4e+08, -1.0e+08, -6.0e+07, -2.0e+07,
+ 2.0e+07, 6.0e+07, 1.0e+08, 1.4e+08, 1.8e+08]
+ )
+
+ # Check if the coordinates are the same
+ np.testing.assert_allclose(
+ act_dt, dt,
+ err_msg='In testXsuiteToBlondTransform the time-coordinate differs'
+ )
+ np.testing.assert_allclose(
+ act_de, de,
+ err_msg='In testXsuiteToBlondTransform the energy-coordinate differs'
+ )
+
+ def testTransformInverse(self):
+ # Arrays of time and energy offsets
+ dts = np.linspace(
+ self.rfstation.t_rf[0, 0] * 1e-3,
+ self.rfstation.t_rf[0, 0], 20
+ )
+ dEs = np.linspace(-200e6, 200e6, 20)
+
+ # Perform transforms to xsuite coordinates and back
+ zeta, ptau = blond_to_xsuite_transform(
+ dts, dEs, self.ring.beta[0, 0], self.ring.energy[0, 0],
+ self.rfstation.omega_rf[0, 0], self.rfstation.phi_s[0]
+ )
+
+ _dts, _dEs = xsuite_to_blond_transform(
+ zeta, ptau, self.ring.beta[0, 0], self.ring.energy[0, 0],
+ self.rfstation.omega_rf[0, 0], self.rfstation.phi_s[0]
+ )
+
+ # Check if the coordinates are the same
+ np.testing.assert_allclose(
+ dts, _dts,
+ err_msg='In testTransformInverse the time-coordinate differs'
+ )
+ np.testing.assert_allclose(
+ dEs, _dEs,
+ err_msg='In testTransformInverse the energy-coordinate differs'
+ )
+
+class TestXsuiteLHC(unittest.TestCase):
+
+ def setUp(self):
+ # Accelerator parameters
+ self.C = 26658.8832 # Machine circumference [m]
+ self.p_s = 450e9 # Synchronous momentum [eV/c]
+ self.p_f = 450.1e9 # Synchronous momentum, final
+ self.h = 35640 # Harmonic number [-]
+ self.alpha = 0.00034849575112251314 # First order mom. comp. factor [-]
+ self.V = 5e6 # RF voltage [V]
+ self.dphi = 0 # Phase modulation/offset [rad]
+
+ def testSingleParticle(self):
+ r'''Test of a single particle simulation in the LHC at constant energy'''
+ rtol = 4e-9 # relative tolerance
+ atol = 5e-9 # absolute tolerance
+
+ # ----- Interface Simulation -----
+ # Initialize the BLonD cavity
+ blond_tracker, beam = self.singleParticleBLonDSimulation()
+ cavity = BlondElement(blond_tracker, beam)
+
+ line = self.initializeXsuiteLine()
+
+ # Insert the BLonD elements
+ line.insert_element(index=0, element=cavity, name='blond_cavity')
+
+ dt_inter, dE_inter = self.performSingleParticleInterfaceSimulation(line, blond_tracker, beam)
+ dt_inter = dt_inter / np.max(dt_inter)
+ dE_inter = dE_inter / np.max(dE_inter)
+
+ # ----- Pure BLonD Simulation -----
+ blond_tracker, beam = self.singleParticleBLonDSimulation()
+ dt_blond, dE_blond = self.performSingleParticleBLonDSimulation(blond_tracker, beam)
+ dt_blond = dt_blond / np.max(dt_blond)
+ dE_blond = dE_blond / np.max(dE_blond)
+
+ # ----- Perform test -----
+ np.testing.assert_allclose(
+ dt_inter, dt_blond, rtol=rtol, atol=atol,
+ err_msg='In testSingleParticle the time-coordinate differs'
+ )
+ np.testing.assert_allclose(
+ dE_inter, dE_blond, rtol=rtol, atol=atol,
+ err_msg='In testSingleParticle the energy-coordinate differs'
+ )
+
+ def testSingleParticleRamp(self):
+ r'''Test of acceleration of a single particle simulation in the LHC'''
+ rtol = 4e-9 # relative tolerance
+ atol = 5e-9 # absolute tolerance
+
+ # Initialize the BLonD cavity
+ blond_tracker, beam = self.singleParticleBLonDSimulation(ramp=True)
+ cavity = BlondElement(blond_tracker, beam)
+
+ line = self.initializeXsuiteLine()
+
+ # Insert the BLonD elements
+ line.insert_element(index=0, element=cavity, name='blond_cavity')
+
+ # Insert energy ramp
+ energy_update = EnergyUpdate(blond_tracker.rf_params.momentum)
+ line.insert_element(index='blond_cavity', element=energy_update, name='energy_update')
+
+ dt_inter, dE_inter = self.performSingleParticleInterfaceSimulation(line, blond_tracker, beam)
+ dt_inter = dt_inter / np.max(dt_inter)
+ dE_inter = dE_inter / np.max(dE_inter)
+ line.get_table().show()
+
+ # ----- Pure BLonD Simulation -----
+ blond_tracker, beam = self.singleParticleBLonDSimulation(ramp=True, pure_blond=True)
+ dt_blond, dE_blond = self.performSingleParticleBLonDSimulation(blond_tracker, beam)
+ dt_blond = dt_blond / np.max(dt_blond)
+ dE_blond = dE_blond / np.max(dE_blond)
+
+ # ----- Perform test -----
+ np.testing.assert_allclose(
+ dt_inter, dt_blond, rtol=rtol, atol=atol,
+ err_msg='In testSingleParticleRamp the time-coordinate differs'
+ )
+ np.testing.assert_allclose(
+ dE_inter, dE_blond, rtol=rtol, atol=atol,
+ err_msg='In testSingleParticleRamp the energy-coordinate differs'
+ )
+
+ def testSingleParticleTwoRFStations(self):
+ r'''Test of a single particle simulation in the LHC at constant energy with two RF stations'''
+ rtol = 4e-9 # relative tolerance
+ atol = 5e-9 # absolute tolerance
+
+ # Initialize the BLonD cavities
+ blond_tracker_1, blond_tracker_2, beam = self.singleParticleBLonDSimulation(two_rfstations=True)
+ cavity_1 = BlondElement(blond_tracker_1, beam)
+ cavity_2 = BlondElement(blond_tracker_2, beam)
+
+ line = self.initializeXsuiteLine(two_rfstations=True)
+
+ # Insert the BLonD elements
+ line.insert_element(index='matrix', element=cavity_1, name='blond_cavity_1')
+ line.insert_element(index='matrix_2', element=cavity_2, name='blond_cavity_2')
+
+ dt_inter, dE_inter = self.performSingleParticleInterfaceSimulation(line, blond_tracker_1, beam)
+ dt_inter = dt_inter / np.max(dt_inter)
+ dE_inter = dE_inter / np.max(dE_inter)
+
+ # ----- Pure BLonD Simulation -----
+ blond_tracker_1, blond_tracker_2, beam = self.singleParticleBLonDSimulation(two_rfstations=True)
+ dt_blond, dE_blond = self.performSingleParticleBLonDSimulation(blond_tracker_1, beam, blond_tracker_2)
+ dt_blond = dt_blond / np.max(dt_blond)
+ dE_blond = dE_blond / np.max(dE_blond)
+
+ # ----- Perform test -----
+ np.testing.assert_allclose(
+ dt_inter, dt_blond, rtol=rtol, atol=atol,
+ err_msg='In testSingleParticleTwoRFStations the time-coordinate differs'
+ )
+ np.testing.assert_allclose(
+ dE_inter, dE_blond, rtol=rtol, atol=atol,
+ err_msg='In testSingleParticleTwoRFStations the energy-coordinate differs'
+ )
+
+ def singleParticleBLonDSimulation(self, ramp: bool = False, two_rfstations: bool = False,
+ pure_blond: bool = False):
+ r'''
+ Method to generate a BLonD simulation in the LHC with a single particle with a time offset with respect to the
+ RF bucket.
+ '''
+ # Bunch parameters
+ N_m = 1 # Number of macroparticles [-]
+ N_p = 1.15e11 # Intensity
+
+ # Simulation parameters
+ N_t = 481 # Number of (tracked) turns [-]
+ dt_offset = 0.15 # Input particles dt [s]
+ v_part = 1
+
+ # --- BLonD objects ---
+ if ramp:
+ cycle = np.linspace(self.p_s, self.p_f, N_t + 1)
+ if pure_blond:
+ cycle = np.concatenate(
+ (np.array([self.p_s]), np.linspace(self.p_s, self.p_f, N_t + 1)[:-1])
+ )
+ ring = Ring(self.C, self.alpha, cycle, Proton(), n_turns=N_t)
+ else:
+ if two_rfstations:
+ # Additional input
+ c_part = 0.3
+ v_part = 0.5
+
+ ring = Ring(
+ [self.C * c_part, self.C * (1 - c_part)],
+ [[self.alpha], [self.alpha]],
+ [self.p_s, self.p_s],
+ Proton(), n_turns=N_t, n_sections=2
+ )
+ else:
+ ring = Ring(self.C, self.alpha, self.p_s, Proton(), n_turns=N_t)
+
+ beam = Beam(ring, N_m, N_p)
+ rfstation = RFStation(ring, [self.h], [self.V * v_part], [self.dphi])
+
+ # --- Insert particle ---
+ beam.dt = np.array([rfstation.t_rf[0, 0] * (1 - dt_offset)])
+ beam.dE = np.array([0.0])
+
+ if two_rfstations:
+ second_rfstation = RFStation(ring, [self.h], [self.V * (1 - v_part)],
+ [self.dphi], section_index=2)
+
+ return (RingAndRFTracker(rfstation, beam),
+ RingAndRFTracker(second_rfstation, beam),
+ beam)
+
+ return RingAndRFTracker(rfstation, beam), beam
+
+ def initializeXsuiteLine(self, two_rfstations: bool = False):
+ r'''
+ Method to generate an Xsuite Line using the one-turn matrix.
+ '''
+ c_part = 1
+
+ if two_rfstations:
+ c_part = 0.3
+ matrix_2 = xt.LineSegmentMap(
+ longitudinal_mode='nonlinear',
+ qx=1.1, qy=1.2,
+ betx=1.,
+ bety=1.,
+ voltage_rf=0,
+ frequency_rf=0,
+ lag_rf=0,
+ momentum_compaction_factor=self.alpha,
+ length=self.C * (1 - c_part)
+ )
+
+ matrix = xt.LineSegmentMap(
+ longitudinal_mode='nonlinear',
+ qx=1.1, qy=1.2,
+ betx=1.,
+ bety=1.,
+ voltage_rf=0,
+ frequency_rf=0,
+ lag_rf=0,
+ momentum_compaction_factor=self.alpha,
+ length=self.C * c_part
+ )
+
+ # Create line
+ line = xt.Line(elements=[matrix], element_names={'matrix'})
+ line['matrix'].length = self.C * c_part
+
+ if two_rfstations:
+ line.append_element(element=matrix_2, name='matrix_2')
+ line['matrix_2'].length = (1 - c_part) * self.C
+
+ return line
+
+ def performSingleParticleInterfaceSimulation(self, line, blond_track, beam):
+ r'''
+ Perform simulation with the interface.
+ '''
+ # Add particles to line and build tracker
+ line.particle_ref = xp.Particles(p0c=self.p_s, mass0=xp.PROTON_MASS_EV, q0=1.)
+ line.build_tracker()
+
+ # --- Convert the initial BLonD distribution to xsuite coordinates ---
+ zeta, ptau = blond_to_xsuite_transform(
+ beam.dt, beam.dE, line.particle_ref.beta0[0],
+ line.particle_ref.energy0[0], phi_s=blond_track.rf_params.phi_s[0],
+ omega_rf=blond_track.rf_params.omega_rf[0, 0]
+ )
+
+ # --- Track matrix ---
+ N_t = len(blond_track.rf_params.phi_s) - 1
+ particles = line.build_particles(x=0, y=0, px=0, py=0, zeta=np.copy(zeta), ptau=np.copy(ptau))
+ line.track(particles, num_turns=N_t, turn_by_turn_monitor=True, with_progress=False)
+ mon = line.record_last_track
+
+ dt_array = np.zeros(N_t)
+ dE_array = np.zeros(N_t)
+ # Convert the xsuite particle coordinates back to BLonD
+ for i in range(N_t):
+ dt_array[i], dE_array[i] = xsuite_to_blond_transform(
+ mon.zeta[:, i].T, mon.ptau[:, i].T, blond_track.rf_params.beta[i], blond_track.rf_params.energy[i],
+ phi_s=blond_track.rf_params.phi_s[i] - blond_track.rf_params.phi_rf[0, 0],
+ omega_rf=blond_track.rf_params.omega_rf[0, i]
+ )
+
+ return dt_array, dE_array
+
+ @staticmethod
+ def performSingleParticleBLonDSimulation(blond_track, beam, second_blond_tracker: RingAndRFTracker = None):
+ r'''
+ Perform simulation in pure BLonD.
+ '''
+ # The number of turns to track
+ N_t = len(blond_track.rf_params.phi_s) - 1
+
+ # Save particle coordinates
+ dt_array = np.zeros(N_t)
+ dE_array = np.zeros(N_t)
+
+ if second_blond_tracker is not None:
+ full_tracker = FullRingAndRF([blond_track, second_blond_tracker])
+ else:
+ full_tracker = blond_track
+
+ for i in range(N_t):
+ dt_array[i] = beam.dt[0]
+ dE_array[i] = beam.dE[0]
+ full_tracker.track()
+
+ return dt_array, dE_array
+
+
+class TestXsuitePSB(unittest.TestCase):
+
+ def setUp(self):
+ # Accelerator parameters
+ self.C = 2 * np.pi * 25.00005794526065 # Machine circumference, radius 25 [m]
+ gamma_t = 4.11635447373496 # Transition gamma [-]
+ self.alpha = 1. / gamma_t / gamma_t # First order mom. comp. factor [-]
+ self.h = 1 # Harmonic number [-]
+ self.V = 8e3 # RF voltage [V]
+ self.dphi = np.pi # Phase modulation/offset [rad]
+
+ # Derived parameters
+ kin_energy = 1.4e9 # Kinetic energy [eV]
+ E_0 = m_p * c ** 2 / e # [eV]
+ tot_energy = E_0 + kin_energy # [eV]
+
+ self.p_s = np.sqrt(tot_energy ** 2 - E_0 ** 2)
+ self.p_f = self.p_s + 0.001e9
+
+ def testSingleParticle(self):
+ r'''Test of a single particle simulation in the PSB at constant energy'''
+ rtol = 4e-4 # relative tolerance
+ atol = 4.1e-4 # absolute tolerance
+
+ # ----- Interface Simulation -----
+ # Initialize the BLonD cavity
+ blond_tracker, beam = self.singleParticleBLonDSimulation()
+ cavity = BlondElement(blond_tracker, beam)
+
+ line = self.initializeXsuiteLine()
+
+ # Insert the BLonD elements
+ line.insert_element(index=0, element=cavity, name='blond_cavity')
+
+ dt_inter, dE_inter = self.performSingleParticleInterfaceSimulation(line, blond_tracker, beam)
+ dt_inter = dt_inter / np.max(dt_inter)
+ dE_inter = dE_inter / np.max(dE_inter)
+
+ # ----- Pure BLonD Simulation -----
+ blond_tracker, beam = self.singleParticleBLonDSimulation()
+ dt_blond, dE_blond = self.performSingleParticleBLonDSimulation(blond_tracker, beam)
+ dt_blond = dt_blond / np.max(dt_blond)
+ dE_blond = dE_blond / np.max(dE_blond)
+
+ # ----- Perform test -----
+ np.testing.assert_allclose(
+ dt_inter, dt_blond, rtol=rtol, atol=atol,
+ err_msg='In testSingleParticle the time-coordinate differs'
+ )
+ np.testing.assert_allclose(
+ dE_inter, dE_blond, rtol=rtol, atol=atol,
+ err_msg='In testSingleParticle the energy-coordinate differs'
+ )
+
+ def testSingleParticleRamp(self):
+ r'''Test of acceleration of a single particle simulation in the PSB'''
+ rtol = 4e-4 # relative tolerance
+ atol = 4.1e-4 # absolute tolerance
+
+ # Initialize the BLonD cavity
+ blond_tracker, beam = self.singleParticleBLonDSimulation(ramp=True)
+ cavity = BlondElement(blond_tracker, beam)
+
+ line = self.initializeXsuiteLine()
+
+ # Insert the BLonD elements
+ line.insert_element(index=0, element=cavity, name='blond_cavity')
+
+ # Insert energy ramp
+ energy_update = EnergyUpdate(blond_tracker.rf_params.momentum)
+ line.insert_element(index='matrix', element=energy_update, name='energy_update')
+
+ dt_inter, dE_inter = self.performSingleParticleInterfaceSimulation(line, blond_tracker, beam)
+ dt_inter = dt_inter / np.max(dt_inter)
+ dE_inter = dE_inter / np.max(dE_inter)
+
+ # ----- Pure BLonD Simulation -----
+ blond_tracker, beam = self.singleParticleBLonDSimulation(ramp=True, pure_blond=True)
+ dt_blond, dE_blond = self.performSingleParticleBLonDSimulation(blond_tracker, beam)
+ dt_blond = dt_blond / np.max(dt_blond)
+ dE_blond = dE_blond / np.max(dE_blond)
+
+ # ----- Perform test -----
+ np.testing.assert_allclose(
+ dt_inter, dt_blond, rtol=rtol, atol=atol,
+ err_msg='In testSingleParticleRamp the time-coordinate differs'
+ )
+ np.testing.assert_allclose(
+ dE_inter, dE_blond, rtol=rtol, atol=atol,
+ err_msg='In testSingleParticleRamp the energy-coordinate differs'
+ )
+
+ def singleParticleBLonDSimulation(self, ramp: bool = False, pure_blond: bool = False):
+ r'''
+ Method to generate a BLonD simulation in the LHC with a single particle with a time offset with respect to the
+ RF bucket.
+ '''
+ # Bunch parameters
+ N_m = 1 # Number of macroparticles [-]
+ N_p = 1e11 # Intensity
+
+ # Simulation parameters
+ N_t = 6000 # Number of (tracked) turns [-]
+ dt_offset = 0.15 # Input particles dt [s]
+
+ # --- BLonD objects ---
+ if ramp:
+ cycle = np.linspace(self.p_s, self.p_f, N_t + 1)
+ if pure_blond:
+ cycle = np.concatenate(
+ (np.array([self.p_s]), np.linspace(self.p_s, self.p_f, N_t + 1)[:-1])
+ )
+ ring = Ring(self.C, self.alpha, cycle, Proton(), n_turns=N_t)
+ else:
+ ring = Ring(self.C, self.alpha, self.p_s, Proton(), n_turns=N_t)
+
+ beam = Beam(ring, N_m, N_p)
+ rfstation = RFStation(ring, [self.h], [self.V], [self.dphi])
+
+ # --- Insert particle ---
+ beam.dt = np.array([rfstation.t_rf[0, 0] * (1 - dt_offset)])
+ beam.dE = np.array([0.0])
+
+ return RingAndRFTracker(rfstation, beam), beam
+
+ def initializeXsuiteLine(self):
+ r'''
+ Method to generate an Xsuite Line using the one-turn matrix.
+ '''
+ matrix = xt.LineSegmentMap(
+ longitudinal_mode='nonlinear',
+ qx=1.1, qy=1.2,
+ betx=1.,
+ bety=1.,
+ voltage_rf=0,
+ frequency_rf=0,
+ lag_rf=0,
+ momentum_compaction_factor=self.alpha,
+ length=self.C
+ )
+
+ # Create line
+ line = xt.Line(elements=[matrix], element_names={'matrix'})
+ line['matrix'].length = self.C
+
+ return line
+
+ def performSingleParticleInterfaceSimulation(self, line, blond_track, beam):
+ r'''
+ Perform simulation with the interface.
+ '''
+ # Add particles to line and build tracker
+ line.particle_ref = xp.Particles(p0c=self.p_s, mass0=xp.PROTON_MASS_EV, q0=1.)
+ line.build_tracker()
+
+ # --- Convert the initial BLonD distribution to xsuite coordinates ---
+ zeta, ptau = blond_to_xsuite_transform(
+ beam.dt, beam.dE, line.particle_ref.beta0[0], line.particle_ref.energy0[0],
+ phi_s=blond_track.rf_params.phi_s[0] - blond_track.rf_params.phi_rf[0, 0],
+ omega_rf=blond_track.rf_params.omega_rf[0, 0]
+ )
+
+ # --- Track matrix ---
+ N_t = len(blond_track.rf_params.phi_s) - 1
+ particles = line.build_particles(x=0, y=0, px=0, py=0, zeta=np.copy(zeta), ptau=np.copy(ptau))
+ line.track(particles, num_turns=N_t, turn_by_turn_monitor=True, with_progress=False)
+ mon = line.record_last_track
+
+ dt_array = np.zeros(N_t)
+ dE_array = np.zeros(N_t)
+ # Convert the xsuite particle coordinates back to BLonD
+ for i in range(N_t):
+ dt_array[i], dE_array[i] = xsuite_to_blond_transform(
+ mon.zeta[:, i].T, mon.ptau[:, i].T, blond_track.rf_params.beta[i], blond_track.rf_params.energy[i],
+ phi_s=blond_track.rf_params.phi_s[i] - blond_track.rf_params.phi_rf[0, 0],
+ omega_rf=blond_track.rf_params.omega_rf[0, i]
+ )
+
+ return dt_array, dE_array
+
+ @staticmethod
+ def performSingleParticleBLonDSimulation(blond_track, beam):
+ r'''
+ Perform simulation in pure BLonD.
+ '''
+ # The number of turns to track
+ N_t = len(blond_track.rf_params.phi_s) - 1
+
+ # Save particle coordinates
+ dt_array = np.zeros(N_t)
+ dE_array = np.zeros(N_t)
+
+ full_tracker = blond_track
+
+ for i in range(N_t):
+ dt_array[i] = beam.dt[0]
+ dE_array[i] = beam.dE[0]
+ full_tracker.track()
+
+ return dt_array, dE_array
\ No newline at end of file