Automatic Kick Coefficients Re-Computation

This example shows how to use the the auto recomputing functionality for kick coefficients in the kick classes. To follow this example please first have a look at the Kinetic Kicks example as well as the relevant FAQ section.

This script will showcase the functionality by following an identical flow to the Kinetic kicks example linked above and expanding on the relevant parts. Demonstration will be done using the CLIC Damping Ring line.

import logging
from dataclasses import dataclass
from typing import Self

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xobjects as xo
import xpart as xp
import xtrack as xt

from xibs.formulary import _bunch_length, _geom_epsx, _geom_epsy, _percent_change, _sigma_delta
from xibs.inputs import BeamParameters, OpticsParameters
from xibs.kicks import KineticKickIBS

logging.basicConfig(
    level=logging.WARNING,
    format="[%(asctime)s] [%(levelname)s] - %(module)s.%(funcName)s:%(lineno)d - %(message)s",
    datefmt="%H:%M:%S",
)
plt.rcParams.update(
    {
        "font.family": "serif",
        "font.size": 20,
        "axes.titlesize": 20,
        "axes.labelsize": 20,
        "xtick.labelsize": 20,
        "ytick.labelsize": 20,
        "legend.fontsize": 15,
        "figure.titlesize": 20,
    }
)

Let’s start by defining the line and particle information, as well as some parameters for later use. Note that bunch_intensity is the actual number of particles in the bunch, and its value influences the IBS kick coefficients calculation while n_part is the number of generated particles for tracking, which is much lower.

line_file = "lines/chrom-corr_DR.newlattice_2GHz.json"
bunch_intensity = int(4.5e9)
n_part = int(1.5e3)
sigma_z = 1.5e-3
nemitt_x = 5.5e-7
nemitt_y = 3.5e-9

Setting up line and particles

Let’s now load the line from file, activate the accelerating cavities and create a context for multithreading with OpenMP, since tracking particles is going to take some time:

line = xt.Line.from_json(line_file)
context = xo.ContextCpu(omp_num_threads="auto")
line.particle_ref = xp.Particles(mass0=xp.ELECTRON_MASS_EV, q0=1, p0c=2.86e9)

# ----- Power accelerating cavities ----- #
for cavity in [element for element in line.elements if isinstance(element, xt.Cavity)]:
    cavity.lag = 180  # we are above transition

line.build_tracker(context, extra_headers=["#define XTRACK_MULTIPOLE_NO_SYNRAD"])
line.optimize_for_tracking()
twiss = line.twiss()

particles = xp.generate_matched_gaussian_bunch(
    num_particles=n_part,
    total_intensity_particles=bunch_intensity,
    nemitt_x=nemitt_x,
    nemitt_y=nemitt_y,
    sigma_z=sigma_z,
    line=line,
)
particles2 = particles.copy()

Loading line from dict:   0%|          | 0/71796 [00:00<?, ?it/s]
Loading line from dict:   4%|▍         | 3024/71796 [00:00<00:02, 30217.05it/s]
Loading line from dict:   9%|▊         | 6110/71796 [00:00<00:02, 30585.03it/s]
Loading line from dict:  13%|█▎        | 9169/71796 [00:00<00:02, 30307.12it/s]
Loading line from dict:  17%|█▋        | 12201/71796 [00:00<00:02, 23869.11it/s]
Loading line from dict:  21%|██▏       | 15287/71796 [00:00<00:02, 26009.02it/s]
Loading line from dict:  26%|██▌       | 18454/71796 [00:00<00:01, 27728.64it/s]
Loading line from dict:  30%|███       | 21698/71796 [00:00<00:01, 29154.49it/s]
Loading line from dict:  35%|███▍      | 24912/71796 [00:00<00:01, 30055.63it/s]
Loading line from dict:  39%|███▉      | 27982/71796 [00:00<00:01, 29399.06it/s]
Loading line from dict:  43%|████▎     | 30968/71796 [00:01<00:01, 28372.53it/s]
Loading line from dict:  47%|████▋     | 33842/71796 [00:01<00:01, 27735.92it/s]
Loading line from dict:  51%|█████     | 36641/71796 [00:01<00:01, 27548.12it/s]
Loading line from dict:  56%|█████▌    | 39848/71796 [00:01<00:01, 28847.85it/s]
Loading line from dict:  60%|██████    | 43100/71796 [00:01<00:00, 29917.56it/s]
Loading line from dict:  64%|██████▍   | 46264/71796 [00:01<00:00, 30423.48it/s]
Loading line from dict:  69%|██████▉   | 49520/71796 [00:01<00:00, 31044.28it/s]
Loading line from dict:  73%|███████▎  | 52635/71796 [00:01<00:00, 24938.70it/s]
Loading line from dict:  78%|███████▊  | 55720/71796 [00:01<00:00, 26436.25it/s]
Loading line from dict:  82%|████████▏ | 58803/71796 [00:02<00:00, 27605.57it/s]
Loading line from dict:  86%|████████▌ | 61890/71796 [00:02<00:00, 28504.11it/s]
Loading line from dict:  90%|█████████ | 64840/71796 [00:02<00:00, 27964.10it/s]
Loading line from dict:  94%|█████████▍| 67707/71796 [00:02<00:00, 27619.45it/s]
Loading line from dict:  98%|█████████▊| 70518/71796 [00:02<00:00, 27410.97it/s]
Loading line from dict: 100%|██████████| 71796/71796 [00:02<00:00, 28052.18it/s]
Done loading line from dict.
Compiling ContextCpu kernels...
Done compiling ContextCpu kernels.
Disable xdeps expressions
Remove markers
Remove inactive multipoles
Merge consecutive multipoles
Remove redundant apertures
Remove zero length drifts
Merge consecutive drifts
Use simple bends
Use simple quadrupoles
Rebuild tracker data
Compiling ContextCpu kernels...
Done compiling ContextCpu kernels.

We can compute initial (geometrical) emittances as well as the bunch length from the xtrack.Particles object:

geom_epsx = _geom_epsx(particles, twiss.betx[0], twiss.dx[0])
geom_epsy = _geom_epsy(particles, twiss.bety[0], twiss.dy[0])
bunch_l = _bunch_length(particles)
sig_delta = _sigma_delta(particles)

Generating the IBS Kick Class

Just as all user-facing classes in xibs, the KineticKickIBS [POA06] class is instantiated by providing a BeamParameters and an OpticsParameters objects:

beamparams = BeamParameters.from_line(line, n_part=bunch_intensity)
opticsparams = OpticsParameters.from_line(line)

# This is the threshold that would trigger the auto-recomping of the growth rates.
AUTO_PERCENT = 8  # threshold: 8% change from kick application
IBS = KineticKickIBS(beamparams, opticsparams)
AUTO_IBS = KineticKickIBS(beamparams, opticsparams, auto_recompute_coefficients_percent=AUTO_PERCENT)

Compiling ContextCpu kernels...
Done compiling ContextCpu kernels.

Computing and Applying IBS Kicks

The so-called “kick coefficients” are computed directly from the particle distribution by calling the dedicated method. They correspond to the diffusion minus friction terms.

IBS.compute_kick_coefficients(particles2)
AUTO_IBS.compute_kick_coefficients(particles2)

Let’s see the evolution of the transverse emittances, \(\sigma_{\delta}\) and bunch length from a kick. We can have a look at the relative change of the particles’s distribution properties from before to after.

IBS.apply_ibs_kick(particles2)
new_geom_epsx = _geom_epsx(particles2, twiss.betx[0], twiss.dx[0])
new_geom_epsy = _geom_epsy(particles2, twiss.bety[0], twiss.dy[0])
new_sig_delta = _sigma_delta(particles2)
new_bunch_length = _bunch_length(particles2)

# Let's see the relative change
print(f"Geom. epsx: {geom_epsx:.2e} -> {new_geom_epsx:.2e} | ({100 * _percent_change(geom_epsx, new_geom_epsx):.2e}% change)")
print(f"Geom. epsy: {geom_epsy:.2e} -> {new_geom_epsy:.2e} | ({100 * _percent_change(geom_epsy, new_geom_epsy):.2e}% change)")
print(f"Sigma delta: {sig_delta:.2e} -> {new_sig_delta:.2e} | ({100 * _percent_change(sig_delta, new_sig_delta):.2e}% change)")
print(f"Bunch length: {bunch_l:.4e} -> {new_bunch_length:.4e} | ({100 * _percent_change(bunch_l, new_bunch_length):.2e}% change)")

# Let's reset these
particles2 = particles.copy()

Geom. epsx: 1.02e-10 -> 1.02e-10 | (-1.26e-06% change)
Geom. epsy: 5.88e-13 -> 5.88e-13 | (0.00e+00% change)
Sigma delta: 1.63e-03 -> 1.63e-03 | (6.77e-01% change)
Bunch length: 1.5163e-03 -> 1.5163e-03 | (0.00e+00% change)

Preparing for the Tracking Simulation

We will loop over turns for 1000 turns, and compare results from the regular and auto-recomputing scenarios. Let’s set up the utilities needed for this:

nturns = 1000  # number of turns to loop for
ibs_step = 50  # frequency at which to re-compute coefficients in [turns]
turns = np.linspace(0, nturns, nturns, dtype=int)  # array of tracked turns


# Set up a dataclass to store the results
@dataclass
class Records:
    epsilon_x: np.ndarray
    epsilon_y: np.ndarray
    sigma_delta: np.ndarray
    bunch_length: np.ndarray

    def update_at_turn(self, turn: int, parts: xp.Particles, twiss: xt.TwissTable):
        self.epsilon_x[turn] = _geom_epsx(parts, twiss.betx[0], twiss.dx[0])
        self.epsilon_y[turn] = _geom_epsy(parts, twiss.bety[0], twiss.dy[0])
        self.sigma_delta[turn] = _sigma_delta(parts)
        self.bunch_length[turn] = _bunch_length(parts)

    @classmethod
    def init_zeroes(cls, n_turns: int) -> Self:  # noqa: F821
        return cls(
            epsilon_x=np.zeros(n_turns, dtype=float),
            epsilon_y=np.zeros(n_turns, dtype=float),
            sigma_delta=np.zeros(n_turns, dtype=float),
            bunch_length=np.zeros(n_turns, dtype=float),
        )


# Initialize the dataclass & store the initial values
regular = Records.init_zeroes(nturns)
regular.update_at_turn(0, particles, twiss)

auto = Records.init_zeroes(nturns)
auto.update_at_turn(0, particles, twiss)

# These arrays we will use to see when the auto-recomputing kicks in
auto_recomputes = np.zeros(nturns, dtype=float)
fixed_recomputes = np.zeros(nturns, dtype=float)

Tracking Evolution Over Turns

Let’s now loop and let the auto-recomputing do its job.

for turn in range(1, nturns):
    # This is not necessary, just for showcasing in this tutorial
    n_recomputes_for_auto = AUTO_IBS._number_of_coefficients_computations

    # ----- Potentially re-compute the IBS kick coefficients ----- #
    if (turn % ibs_step == 0) or (turn == 1):
        print(f"Turn {turn:d}: Fixed interval re-computing of coefficients")
        # Below always re-computes the coefficients every 'ibs_step' turns
        IBS.compute_kick_coefficients(particles)

    # ----- Manually Apply IBS Kick and Track Turn ----- #
    IBS.apply_ibs_kick(particles)
    AUTO_IBS.apply_ibs_kick(particles2)  # auto-recomputes if necessary
    line.track(particles, num_turns=1)
    line.track(particles2, num_turns=1)

    # ----- Update records for tracked particles ----- #
    regular.update_at_turn(turn, particles, twiss)
    auto.update_at_turn(turn, particles2, twiss)

    # ----- Check if the rates were auto-recomputed ----- #
    # This is also not necessary, just for showcasing in this tutorial
    if AUTO_IBS._number_of_coefficients_computations > n_recomputes_for_auto:
        print(f"At turn {turn} - Auto re-computed coefficients")
        auto_recomputes[turn - 1] = 1


# Here we also aggregate the fixed recomputes
for turn in turns:
    if (turn % ibs_step == 0) or (turn == 1):
        fixed_recomputes[turn - 1] = 1

# And these are simply 1D arrays of the turns at which re-computing happened
where_auto_recomputes = np.flatnonzero(auto_recomputes)
where_fixed_recomputes = np.flatnonzero(fixed_recomputes)

Turn 1: Fixed interval re-computing of coefficients
At turn 9 - Auto re-computed coefficients
At turn 10 - Auto re-computed coefficients
At turn 12 - Auto re-computed coefficients
At turn 14 - Auto re-computed coefficients
At turn 18 - Auto re-computed coefficients
At turn 35 - Auto re-computed coefficients
At turn 39 - Auto re-computed coefficients
At turn 43 - Auto re-computed coefficients
At turn 44 - Auto re-computed coefficients
At turn 45 - Auto re-computed coefficients
At turn 47 - Auto re-computed coefficients
Turn 50: Fixed interval re-computing of coefficients
At turn 53 - Auto re-computed coefficients
At turn 54 - Auto re-computed coefficients
At turn 57 - Auto re-computed coefficients
At turn 58 - Auto re-computed coefficients
At turn 66 - Auto re-computed coefficients
At turn 75 - Auto re-computed coefficients
At turn 83 - Auto re-computed coefficients
At turn 91 - Auto re-computed coefficients
At turn 99 - Auto re-computed coefficients
Turn 100: Fixed interval re-computing of coefficients
At turn 115 - Auto re-computed coefficients
At turn 118 - Auto re-computed coefficients
At turn 119 - Auto re-computed coefficients
At turn 122 - Auto re-computed coefficients
At turn 125 - Auto re-computed coefficients
At turn 132 - Auto re-computed coefficients
At turn 135 - Auto re-computed coefficients
At turn 141 - Auto re-computed coefficients
At turn 146 - Auto re-computed coefficients
Turn 150: Fixed interval re-computing of coefficients
At turn 152 - Auto re-computed coefficients
At turn 158 - Auto re-computed coefficients
At turn 159 - Auto re-computed coefficients
At turn 165 - Auto re-computed coefficients
At turn 166 - Auto re-computed coefficients
At turn 167 - Auto re-computed coefficients
At turn 170 - Auto re-computed coefficients
At turn 173 - Auto re-computed coefficients
At turn 177 - Auto re-computed coefficients
At turn 199 - Auto re-computed coefficients
Turn 200: Fixed interval re-computing of coefficients
At turn 204 - Auto re-computed coefficients
At turn 217 - Auto re-computed coefficients
At turn 227 - Auto re-computed coefficients
At turn 238 - Auto re-computed coefficients
At turn 243 - Auto re-computed coefficients
At turn 244 - Auto re-computed coefficients
At turn 245 - Auto re-computed coefficients
Turn 250: Fixed interval re-computing of coefficients
At turn 264 - Auto re-computed coefficients
At turn 265 - Auto re-computed coefficients
At turn 297 - Auto re-computed coefficients
Turn 300: Fixed interval re-computing of coefficients
At turn 318 - Auto re-computed coefficients
At turn 331 - Auto re-computed coefficients
At turn 333 - Auto re-computed coefficients
At turn 338 - Auto re-computed coefficients
At turn 340 - Auto re-computed coefficients
At turn 343 - Auto re-computed coefficients
At turn 344 - Auto re-computed coefficients
Turn 350: Fixed interval re-computing of coefficients
At turn 356 - Auto re-computed coefficients
At turn 362 - Auto re-computed coefficients
At turn 374 - Auto re-computed coefficients
At turn 377 - Auto re-computed coefficients
At turn 389 - Auto re-computed coefficients
At turn 391 - Auto re-computed coefficients
At turn 392 - Auto re-computed coefficients
Turn 400: Fixed interval re-computing of coefficients
At turn 401 - Auto re-computed coefficients
At turn 405 - Auto re-computed coefficients
At turn 415 - Auto re-computed coefficients
At turn 422 - Auto re-computed coefficients
At turn 424 - Auto re-computed coefficients
At turn 430 - Auto re-computed coefficients
At turn 432 - Auto re-computed coefficients
Turn 450: Fixed interval re-computing of coefficients
At turn 453 - Auto re-computed coefficients
At turn 454 - Auto re-computed coefficients
At turn 460 - Auto re-computed coefficients
At turn 476 - Auto re-computed coefficients
At turn 478 - Auto re-computed coefficients
At turn 481 - Auto re-computed coefficients
At turn 497 - Auto re-computed coefficients
Turn 500: Fixed interval re-computing of coefficients
At turn 504 - Auto re-computed coefficients
At turn 505 - Auto re-computed coefficients
At turn 507 - Auto re-computed coefficients
At turn 511 - Auto re-computed coefficients
At turn 514 - Auto re-computed coefficients
At turn 517 - Auto re-computed coefficients
At turn 530 - Auto re-computed coefficients
At turn 534 - Auto re-computed coefficients
At turn 543 - Auto re-computed coefficients
Turn 550: Fixed interval re-computing of coefficients
At turn 552 - Auto re-computed coefficients
At turn 558 - Auto re-computed coefficients
At turn 561 - Auto re-computed coefficients
At turn 574 - Auto re-computed coefficients
At turn 577 - Auto re-computed coefficients
At turn 581 - Auto re-computed coefficients
At turn 584 - Auto re-computed coefficients
At turn 593 - Auto re-computed coefficients
Turn 600: Fixed interval re-computing of coefficients
At turn 602 - Auto re-computed coefficients
At turn 606 - Auto re-computed coefficients
At turn 627 - Auto re-computed coefficients
At turn 629 - Auto re-computed coefficients
At turn 638 - Auto re-computed coefficients
At turn 643 - Auto re-computed coefficients
At turn 644 - Auto re-computed coefficients
Turn 650: Fixed interval re-computing of coefficients
At turn 650 - Auto re-computed coefficients
At turn 651 - Auto re-computed coefficients
At turn 662 - Auto re-computed coefficients
At turn 675 - Auto re-computed coefficients
At turn 681 - Auto re-computed coefficients
At turn 684 - Auto re-computed coefficients
At turn 687 - Auto re-computed coefficients
Turn 700: Fixed interval re-computing of coefficients
At turn 700 - Auto re-computed coefficients
At turn 703 - Auto re-computed coefficients
At turn 729 - Auto re-computed coefficients
Turn 750: Fixed interval re-computing of coefficients
At turn 757 - Auto re-computed coefficients
At turn 761 - Auto re-computed coefficients
At turn 790 - Auto re-computed coefficients
At turn 793 - Auto re-computed coefficients
At turn 795 - Auto re-computed coefficients
Turn 800: Fixed interval re-computing of coefficients
At turn 818 - Auto re-computed coefficients
At turn 819 - Auto re-computed coefficients
Turn 850: Fixed interval re-computing of coefficients
At turn 853 - Auto re-computed coefficients
At turn 857 - Auto re-computed coefficients
At turn 860 - Auto re-computed coefficients
At turn 864 - Auto re-computed coefficients
At turn 872 - Auto re-computed coefficients
At turn 873 - Auto re-computed coefficients
At turn 892 - Auto re-computed coefficients
Turn 900: Fixed interval re-computing of coefficients
At turn 913 - Auto re-computed coefficients
At turn 920 - Auto re-computed coefficients
At turn 937 - Auto re-computed coefficients
At turn 947 - Auto re-computed coefficients
Turn 950: Fixed interval re-computing of coefficients
At turn 955 - Auto re-computed coefficients
At turn 957 - Auto re-computed coefficients
At turn 968 - Auto re-computed coefficients
At turn 978 - Auto re-computed coefficients
At turn 993 - Auto re-computed coefficients
At turn 997 - Auto re-computed coefficients

Let’s see how many recomputes happened in total for each scenario:

print(f"Fixed re-computes: {IBS._number_of_coefficients_computations}")
print(f"Auto re-computes: {AUTO_IBS._number_of_coefficients_computations}")

Fixed re-computes: 21
Auto re-computes: 136

That’s almost four times as many updates of the coeffifients that were deemed necessary by the auto-recomputing mechanism. Let’s see the effect it has had on our results by having a look at the evolution of emittances over time:

AVG_TURNS = 3  # for rolling average plotting
fig, axs = plt.subplot_mosaic([["epsx", "epsy"], ["sigd", "bl"]], sharex=True, figsize=(13, 7))

# We will add vertical lines at the times where recomputing of the kick
# coefficients happened (do this first so they show up in the background)
for axis in axs.values():
    for turn in where_fixed_recomputes:
        axis.axvline(turn, color="gray", linestyle="--", lw=1, alpha=0.7)
    for turn in where_auto_recomputes:
        axis.axvline(turn, color="C1", linestyle="-", alpha=0.2)

axs["epsx"].plot(turns, 1e10 * pd.Series(regular.epsilon_x).rolling(AVG_TURNS, closed="both", min_periods=1).mean(), lw=2, label=f"Fixed ({ibs_step}) turns")
axs["epsy"].plot(turns, 1e13 * pd.Series(regular.epsilon_y).rolling(AVG_TURNS, closed="both", min_periods=1).mean(), lw=2, label=f"Fixed ({ibs_step}) turns")
axs["sigd"].plot(turns, 1e3 * pd.Series(regular.sigma_delta).rolling(AVG_TURNS, closed="both", min_periods=1).mean(), lw=2, label=f"Fixed ({ibs_step}) turns")
axs["bl"].plot(turns, 1e3 * pd.Series(regular.bunch_length).rolling(AVG_TURNS, closed="both", min_periods=1).mean(), lw=2, label=f"Fixed ({ibs_step}) turns")

axs["epsx"].plot(turns, 1e10 * pd.Series(auto.epsilon_x).rolling(AVG_TURNS, closed="both", min_periods=1).mean(), lw=1.9, label=f"Auto ({AUTO_PERCENT:.0f}% change)")
axs["epsy"].plot(turns, 1e13 * pd.Series(auto.epsilon_y).rolling(AVG_TURNS, closed="both", min_periods=1).mean(), lw=1.9, label=f"Auto ({AUTO_PERCENT:.0f}% change)")
axs["sigd"].plot(turns, 1e3 * pd.Series(auto.sigma_delta).rolling(AVG_TURNS, closed="both", min_periods=1).mean(), lw=1.9, label=f"Auto ({AUTO_PERCENT:.0f}% change)")
axs["bl"].plot(turns, 1e3 * pd.Series(auto.bunch_length).rolling(AVG_TURNS, closed="both", min_periods=1).mean(), lw=1.9, label=f"Auto ({AUTO_PERCENT:.0f}% change)")

# Axes parameters
axs["epsx"].set_ylabel(r"$\varepsilon_x$ [$10^{-10}$m]")
axs["epsy"].set_ylabel(r"$\varepsilon_y$ [$10^{-13}$m]")
axs["sigd"].set_ylabel(r"$\sigma_{\delta}$ [$10^{-3}$]")
axs["bl"].set_ylabel(r"Bunch length [mm]")

for axis in (axs["epsy"], axs["bl"]):
    axis.yaxis.set_label_position("right")
    axis.yaxis.tick_right()

for axis in (axs["sigd"], axs["bl"]):
    axis.set_xlabel("Turn Number")

for axis in axs.values():
    axis.yaxis.set_major_locator(plt.MaxNLocator(3))

fig.align_ylabels((axs["epsx"], axs["sigd"]))
fig.align_ylabels((axs["epsy"], axs["bl"]))
fig.suptitle(f"Evolution of Emittances in CLIC DR (Kinetic Kicks)\nAverage of Last {AVG_TURNS} Turns")

plt.legend(title="Recompute Coefficients")
plt.tight_layout()
plt.show()

We can have a look at the relative change from turn to turn:

fig, axs = plt.subplot_mosaic([["epsx", "epsy"], ["sigd", "bl"]], sharex=True, figsize=(13, 7))

# We will add vertical lines at the times where recomputing of the kick
# coefficients happened (do this first so they show up in the background)
for axis in axs.values():
    for turn in where_fixed_recomputes:
        axis.axvline(turn, color="gray", linestyle="--", lw=1, alpha=0.7)
    for turn in where_auto_recomputes:
        axis.axvline(turn, color="C1", linestyle="-", alpha=0.2)

axs["epsx"].plot(turns, 1e2 * pd.Series(regular.epsilon_x).pct_change(), lw=1, label=f"Fixed ({ibs_step}) turns")
axs["epsy"].plot(turns, 1e2 * pd.Series(regular.epsilon_y).pct_change(), lw=1, label=f"Fixed ({ibs_step}) turns")
axs["sigd"].plot(turns, 1e2 * pd.Series(regular.sigma_delta).pct_change(), lw=1, label=f"Fixed ({ibs_step}) turns")
axs["bl"].plot(turns, 1e2 * pd.Series(regular.bunch_length).pct_change(), lw=1, label=f"Fixed ({ibs_step}) turns")

axs["epsx"].plot(turns, 1e2 * pd.Series(auto.epsilon_x).pct_change(), lw=1, label=f"Auto ({AUTO_PERCENT:.0f}% change)")
axs["epsy"].plot(turns, 1e2 * pd.Series(auto.epsilon_y).pct_change(), lw=1, label=f"Auto ({AUTO_PERCENT:.0f}% change)")
axs["sigd"].plot(turns, 1e2 * pd.Series(auto.sigma_delta).pct_change(), lw=1, label=f"Auto ({AUTO_PERCENT:.0f}% change)")
axs["bl"].plot(turns, 1e2 * pd.Series(auto.bunch_length).pct_change(), lw=1, label=f"Auto ({AUTO_PERCENT:.0f}% change)")

# Axes parameters
axs["epsx"].set_ylabel(r"$\varepsilon_x$")
axs["epsy"].set_ylabel(r"$\varepsilon_y$")
axs["sigd"].set_ylabel(r"$\sigma_{\delta}$")
axs["bl"].set_ylabel(r"Bunch length")

for axis in (axs["epsy"], axs["bl"]):
    axis.yaxis.set_label_position("right")
    axis.yaxis.tick_right()

for axis in (axs["sigd"], axs["bl"]):
    axis.set_xlabel("Duration [h]")

for axis in axs.values():
    axis.axhline(AUTO_PERCENT, color="black", linestyle="--", alpha=0.5, label="Recompute Threshold")
    axis.yaxis.set_major_locator(plt.MaxNLocator(3))
    axis.set_yscale("log")

fig.align_ylabels((axs["epsx"], axs["sigd"]))
fig.align_ylabels((axs["epsy"], axs["bl"]))
fig.suptitle("Percent change from previous turn")

plt.legend(title="Recompute Coefficients")
plt.tight_layout()
plt.show()

We can see that the auto-recompute feature leads to correct results of the IBS effects in addition to making the tracking loop simpler to write. It also avoids for the user to have to determine a proper interval for re-computing the kicks - which could be correct for part of the simulation but not all of it - by adapting itself to the actual changes in the particle distribution.

References

The use of the following functions, methods, classes and modules is shown in this example:

• analytical: NagaitsevIBS, growth_rates, emittance_evolution

• inputs: BeamParameters, OpticsParameters

• kicks: DiffusionCoefficients, FrictionCoefficients, IBSKickCoefficients, KineticKickIBS, apply_ibs_kick, compute_kick_coefficients