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¶

pip install mbtrack2

Tracking set-up¶

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

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")

../_images/2e1b2fcbeb001884f9f28a45e7b091c85f636a2d73f6a859de7d692c272bbd8d.png

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

../_images/775b4958d7952d2cd56ba76949766f91650318329cf8652cec2d3dac63b023a0.png

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 0x7658d5422300>]

../_images/02b14caea0d9675c92e82178d5802edbc1b903a53bd7e1f10a815e3c19f8233d.png

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()

../_images/cd26b53d9d1b9518fe1d894d9d4900cf4db468a56e5b0ccf376b530cb951fc09.png

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")

../_images/58904201affa26255f2a13b1834339d9eebcb626f551cf2c2c6d8d5d6cb315b3.png

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

../_images/16e3ff3ab2eeed1536b390c034d5ae5449341219f735aa0488229b3e98f05950.png

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 0x7658d31f7f20>]

../_images/2e57a4c4458c2f6b421dc7fbba54157fbd268b128884eaf1c00286aacd2ebe6a.png

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()

../_images/85ef66e903b269417f86c0da70e6ca24131c0c593918668b8f19636ab2b6ecf5.png

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")

../_images/3c6625318ca1a37b90c250b0f493ff299c4a43aab97a471a8a5c838a9a53a60f.png

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

../_images/4382521a1e5c26faf4b553f828ea841e4dbbf4702c77342f5b9e89472bb36c4b.png

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")

../_images/522c067a743bce3e0e400e266300de625584b0b80ff87390270577b6f8ef57d6.png

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

../_images/db8ec76ad417173cdfae96c55649db49c3e07a1e01a48bf9d330deb5d433e4ae.png

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 0x7658d2f46d50>]

../_images/121060b38dbae1d8f3077df405b61b47849815ec1b1268e74fb9a0be73ad7ad0.png

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")

../_images/fdf54eab568deff48df981c80efb92472846b8519d7652c7c82dfe74196a00d0.png

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

../_images/d504c4a2149405fecf84d25e6f055e31ccfa55ddd47849c41b3c37bae6f18d70.png

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.956954332114212 turns.

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

../_images/e4e16a9eac6cdb74ae3b6325f40fe6e0cb2d55c6f17d4139d76aa3f3c2617174.png