Bunch-by-bunch damper feedback

This notebook introduces the TransverseExponentialDamper and FIRDamper classes for mbtrack2:

• The TransverseExponentialDamper is a very simple bunch-by-bunch damper feedback class based on exponential damping.

• The FIRDamper is a more realistic bunch-by-bunch damper feedback class based on FIR filters.

Installation

! git clone -b develop https://gitlab.synchrotron-soleil.fr/PA/collective-effects/mbtrack2.git

Cloning into 'mbtrack2'...
remote: Enumerating objects: 1526, done.
remote: Counting objects: 100% (186/186), done.
remote: Compressing objects: 100% (186/186), done.
remote: Total 1526 (delta 118), reused 0 (delta 0), pack-reused 1340
Receiving objects: 100% (1526/1526), 709.26 KiB | 1024.00 KiB/s, done.
Resolving deltas: 100% (1047/1047), done.

%cd mbtrack2

/content/mbtrack2

Tracking set-up

import numpy as np
import matplotlib.pyplot as plt
import h5py as hp
from mbtrack2.tracking import Synchrotron, Electron
from mbtrack2.utilities import Optics
from mbtrack2.tracking import Bunch
from mbtrack2.tracking import LongitudinalMap, RFCavity, SynchrotronRadiation, TransverseMap

mbtrack2 version 0.8.0.37
(dirty git work tree)
--------------------------------------------------
If used in a publication, please cite mbtrack2 paper and the zenodo archive for the corresponding code version (and other papers for more specific features).
[1] A. Gamelin, W. Foosang, N. Yamamoto, V. Gubaidulin and R. Nagaoka, “mbtrack2”. Zenodo, Dec. 16, 2024. doi: 10.5281/zenodo.14418989.
[2] A. Gamelin, W. Foosang, and R. Nagaoka, “mbtrack2, a Collective Effect Library in Python”, presented at the 12th Int. Particle Accelerator Conf. (IPAC'21), Campinas, Brazil, May 2021, paper MOPAB070.

h = 20 # Harmonic number of the accelerator.
L = 100 # Ring circumference in [m].
E0 = 1.5e9 # Nominal (total) energy of the ring in [eV].
particle = Electron() # Particle considered.
ac = 1e-3 # Momentum compaction factor.
U0 = 200e3 # Energy loss per turn in [eV].
tau = np.array([1e-3, 1e-3, 2e-3]) # Horizontal, vertical and longitudinal damping times in [s].
tune = np.array([12.2, 15.3]) # Horizontal and vertical tunes.
emit = np.array([10e-9, 10e-12]) # Horizontal and vertical equilibrium emittance in [m.rad].
sigma_0 = 15e-12 # Natural bunch length in [s].
sigma_delta = 1e-3 # Equilibrium energy spread.
chro = [2.0, 3.0] # Horizontal and vertical (non-normalized) chromaticities.

local_beta = np.array([3, 2]) # Beta function at the tracking location.
local_alpha = np.array([0, 0]) # Alpha function at the tracking location.
local_dispersion = np.array([0, 0, 0, 0]) # Dispersion function and its derivative at the tracking location.
optics = Optics(local_beta=local_beta, local_alpha=local_alpha, 
                  local_dispersion=local_dispersion)

ring = Synchrotron(h=h, optics=optics, particle=particle, L=L, E0=E0, ac=ac, 
                   U0=U0, tau=tau, emit=emit, tune=tune, 
                   sigma_delta=sigma_delta, sigma_0=sigma_0, chro=chro)

LongMap = LongitudinalMap(ring)
RF = RFCavity(ring, m=1, Vc=1e6, theta=np.arccos(ring.U0/1e6))
SR = SynchrotronRadiation(ring)
TransMap = TransverseMap(ring)

ExponentialDamper

from mbtrack2 import TransverseExponentialDamper, BunchMonitor

Initialize a single particle with a 1 mm horizontal offset

mybunch = Bunch(ring, mp_number=1, current=1e-3)
mybunch["x"] += 1e-3

Initialize an TransverseExponentialDamper in the horizontal plane with a damping time of 50 turns and a phase of 90° (fully resistive damper).

damper = TransverseExponentialDamper(ring=ring, damping_time=[50, 0], phase_diff=[90, 0])

bunchmon = BunchMonitor(bunch_number=0, save_every=1, buffer_size=10, total_size=1000, file_name="damper")

turns = 1000
for i in range(turns):
  LongMap.track(mybunch)
  TransMap.track(mybunch)
  RF.track(mybunch)
  damper.track(mybunch)
  bunchmon.track(mybunch)
bunchmon.close()

Plot the particle position and Courant-Snyder invariant:

from mbtrack2 import plot_bunchdata

fig = plot_bunchdata("damper.hdf5",0,"mean","x")

fig = plot_bunchdata("damper.hdf5",0,"cs_invariant","x")

Compute the damping time from tracking data in turn number:

from scipy.optimize import curve_fit

g = hp.File("damper.hdf5","r")
xdata = np.array(g["BunchData_0"]["time"])*ring.T0
ydata = np.array(g["BunchData_0"]["cs_invariant"][0,:])
g.close()

def func(x, a, b):
    return a * np.exp(-b * x)

popt, pcov = curve_fit(func, xdata, ydata)
print(f"Fitted damping time is {1/popt[1]/ring.T0*2} turns.")
print(f"Input damping time is {damper.damping_time[0]} turns.")
plt.plot(xdata, ydata)
plt.plot(xdata, func(xdata, *popt),"--")

Fitted damping time is 49.06419715611718 turns.
Input damping time is 50 turns.

[<matplotlib.lines.Line2D at 0x7fd6a152a150>]

FIRDamper

from mbtrack2 import FIRDamper 

Stable case (low gain)

mybunch = Bunch(ring, mp_number=1, current=1e-3, track_alive=False)
mybunch["x"] += 1e-3

Initialize a FIRDamper element on the horizontal plane.

The FIR filter is computed for a reference tune of 0.2 (usually the same as ring.tune) with a gain of 0.1 and a phase of 90°. A delay of 1 turn between the measurement and the kick is also taken into account.

If the gain is low enough, the FIRDamper element should provide damping.

fir_damp = FIRDamper(ring=ring, plane="x", tune=0.2, turn_delay=1, tap_number=5, 
                     gain=0.1, phase=90, bpm_error=None, max_kick=None)

The FIR filter computed from the input can then be plotted:

fig = fir_damp.plot_fir()

bunchmon = BunchMonitor(bunch_number=0, save_every=1, buffer_size=10, total_size=100, file_name="fir")

turns = 100
for i in range(turns):
  LongMap.track(mybunch)
  TransMap.track(mybunch)
  RF.track(mybunch)
  fir_damp.track(mybunch)
  bunchmon.track(mybunch)
bunchmon.close()

fig = plot_bunchdata("fir.hdf5",0,"mean","x")

fig = plot_bunchdata("fir.hdf5",0,"cs_invariant","x")

g = hp.File("fir.hdf5","r")
xdata = np.array(g["BunchData_0"]["time"])*ring.T0
ydata = np.array(g["BunchData_0"]["cs_invariant"][0,:])
g.close()

def func(x, a, b):
    return a * np.exp(-b * x)

popt, pcov = curve_fit(func, xdata, ydata)
print(f"Fitted damping time is {1/popt[1]/ring.T0*2} turns.")
plt.plot(xdata, ydata)
plt.plot(xdata, func(xdata, *popt),"--")

Fitted damping time is 6.990584501564287 turns.

[<matplotlib.lines.Line2D at 0x7fd6a11dd220>]

Unstable case (high gain)

mybunch = Bunch(ring, mp_number=1, current=1e-3, track_alive=False)
mybunch["x"] += 1e-3

If the gain is too high (here gain = 0.5 is enough), the FIRDamper leads to unstable dynamics.

fir_damp = FIRDamper(ring=ring, plane="x", tune=0.2, turn_delay=1, tap_number=5, 
                     gain=0.5, phase=90, bpm_error=None, max_kick=None)

fig = fir_damp.plot_fir()

bunchmon = BunchMonitor(bunch_number=0, save_every=1, buffer_size=10, total_size=100, file_name="fir")

turns = 100
for i in range(turns):
  LongMap.track(mybunch)
  TransMap.track(mybunch)
  RF.track(mybunch)
  fir_damp.track(mybunch)
  bunchmon.track(mybunch)
bunchmon.close()

fig = plot_bunchdata("fir.hdf5",0,"mean","x")

fig = plot_bunchdata("fir.hdf5",0,"cs_invariant","x")

Physical limitations

Physical limitations such as the maximum kick strength or the BPM measurement errors can possibly change the damping time.

Maximum kick limit

mybunch = Bunch(ring, mp_number=1, current=1e-3, track_alive=False)
mybunch["x"] += 1e-3

Here the max_kick of the FIRDamper is set to 0.1 mrad:

fir_damp = FIRDamper(ring=ring, plane="x", tune=0.2, turn_delay=1, tap_number=5, 
                     gain=0.1, phase=90, bpm_error=None, max_kick=0.1e-3)

bunchmon = BunchMonitor(bunch_number=0, save_every=1, buffer_size=10, total_size=100, file_name="fir")

turns = 100
for i in range(turns):
  LongMap.track(mybunch)
  TransMap.track(mybunch)
  RF.track(mybunch)
  fir_damp.track(mybunch)
  bunchmon.track(mybunch)
bunchmon.close()

fig = plot_bunchdata("fir.hdf5",0,"mean","x")

fig = plot_bunchdata("fir.hdf5",0,"cs_invariant","x")

g = hp.File("fir.hdf5","r")
xdata = np.array(g["BunchData_0"]["time"])*ring.T0
ydata = np.array(g["BunchData_0"]["cs_invariant"][0,:])
g.close()

def func(x, a, b):
    return a * np.exp(-b * x)

The damping time is increased from 7 turns to 11 turns due to the kick strength limitation.

popt, pcov = curve_fit(func, xdata, ydata)
print(f"Fitted damping time is {1/popt[1]/ring.T0*2} turns.")
plt.plot(xdata, ydata)
plt.plot(xdata, func(xdata, *popt),"--")

Fitted damping time is 10.853769407288695 turns.

[<matplotlib.lines.Line2D at 0x7fd6a103be90>]

BPM measurement error

mybunch = Bunch(ring, mp_number=1, current=1e-3, track_alive=False)
mybunch["x"] += 1e-3

Here a random measurement error of 50 um RMS is added to the bunch mean position used for the feedback correction.

fir_damp = FIRDamper(ring=ring, plane="x", tune=0.2, turn_delay=1, tap_number=5, 
                     gain=0.1, phase=90, bpm_error=50e-6, max_kick=None)

bunchmon = BunchMonitor(bunch_number=0, save_every=1, buffer_size=10, total_size=100, file_name="fir")

turns = 100
for i in range(turns):
  LongMap.track(mybunch)
  TransMap.track(mybunch)
  RF.track(mybunch)
  fir_damp.track(mybunch)
  bunchmon.track(mybunch)
bunchmon.close()

fig = plot_bunchdata("fir.hdf5",0,"mean","x")

fig = plot_bunchdata("fir.hdf5",0,"cs_invariant","x")

g = hp.File("fir.hdf5","r")
xdata = np.array(g["BunchData_0"]["time"])*ring.T0
ydata = np.array(g["BunchData_0"]["cs_invariant"][0,:])
g.close()

def func(x, a, b):
    return a * np.exp(-b * x)

popt, pcov = curve_fit(func, xdata, ydata)
print(f"Fitted damping time is {1/popt[1]/ring.T0*2} turns.")
plt.plot(xdata, ydata)
plt.plot(xdata, func(xdata, *popt),"--")

Fitted damping time is 7.024288846804304 turns.

[<matplotlib.lines.Line2D at 0x7fd6a1038080>]