Skip to content
Snippets Groups Projects
Commit 7d37efa0 authored by SAURET's avatar SAURET
Browse files

Feedback class vf

parent 3e8ac74f
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
mbtrack2/tracking/FeedBack.py
+ 241
− 113
View file @ 7d37efa0
# -*- coding: utf-8 -*-
"""
Created on Tue May 3 09:36:24 2022
@author: sauret
This module defines both simple and FIR based feedback for tracking.
"""
from mbtrack2.tracking.element import Element
import numpy as np
import matplotlib.pyplot as plt
class FeedBackTransverse(Element) :
class SimpleFB(Element) :
"""
Transverse feedback system to counter transverse instabilities.
Simple feedback system proceeding bunch by bunch correction to counter beam
instabilities
Parameters
----------
ring : Synchrotron object
turn_delay : int, optional
Number of turn delay before applying kick.
tap_number : int, optional
Number of tap for the FIR filters.
ampl : float, optional
Gain of the FIR filters.
phase : float, optional
Phase of the FIR filters in deg.
damping_time : array of shape (3,)
Contains the desired damping time in seconds for each plane. Element 0 of
the array defines the horizontal plane as 1 and 2 define the vertical
and longitudinal plane respectively.
phase_diff : array of shape (3,)
Contains the desired phase in degree for each plane. Element 0 of the
array defines the horizontal plane as 1 and 2 define the vertical and
longitudinal plane respectively.
plane : bool array of shape (3,), optional
Allow to choose on which plane the feedback is active. State True to
turn on the feedback or False otherwise. Element 0 of the array defines
the horizontal plane as 1 and 2 define the vertical and longitudinal plane
respectively.
"""
def __init__(self, ring, damping_time, phase_diff, plane = np.zeros((3,), dtype= bool)):
self.ring = ring
self.damping_time = damping_time
self.phase_diff = phase_diff
self.plane = plane
@Element.parallel
def track(self, bunch):
"""
Tracking method for the feedback system
No bunch to bunch interaction, so written for Bunch object and
@Element.parallel is used to handle Beam object.
Parameters
----------
bunch : Bunch or Beam object
"""
if(self.plane[0] == True):
bunch["xp"] -= (2*self.ring.T0/self.damping_time[0])*np.sin(self.phase_diff[0])*bunch.mean[1]
if(self.plane[1] == True):
bunch["yp"] -= (2*self.ring.T0/self.damping_time[1])*np.sin(self.phase_diff[1])*bunch.mean[3]
if(self.plane[2] == True):
bunch["delta"] -= (2*self.ring.T0/self.damping_time[2])*np.sin(self.phase_diff[1])*bunch.mean[5]
class FIR_Feedback(Element) :
"""
Feedback system based on FIR filters, proceeding bunch by bunch correction
to counter beam instabilities.
Parameters
----------
ring : Synchrotron object
tune : float array of shape (3,)
Betatron tune for horizontal and vertical plane in the first and second
elements of the array. Synchrotron tune for the third element.
turn_delay : int array of shape (3,)
Number of turn delay before applying kick. Element 0 of the array
defines the horizontal plane as 1 and 2 define the vertical and
longitudinal plane respectively.
tap_number : int array of shape (3,)
Number of tap for the FIR filters. Element 0 of the array defines
the horizontal plane as 1 and 2 define the vertical and longitudinal
plane respectively.
gain : float array of shape (3,)
Gain of the FIR filters. Element 0 of the array defines
the horizontal plane as 1 and 2 define the vertical and longitudinal
plane respectively.
phase : float array of shape (3,)
Phase of the FIR filters in degree. Element 0 of the array defines
the horizontal plane as 1 and 2 define the vertical and longitudinal
plane respectively.
bpm_error : float array of shape (3,), optional
Stores the BPM error [m] for horizontal and vertical, [s] for longitudinal.
First, second and third element corespond to the horizontal, vertical
and longitudinal plane respectively.
max_kick : float array of shape (3,), optional
Stores the maximum kick limitation. First, second and third
element correspond to the horizontal, vertical and longitudinal plane
respectively.
max_kick_state : bool array of shape (3,), optional
Allow to choose on which plane the max kick limitation is active. State
True to turn on the max kick limitation or False otherwise. Element 0
of the array defines the horizontal plane as 1 and 2 define the vertical
and longitudinal plane respectively.
plane : bool array of shape (3,), optional
Allow to choose on which plane the feedback is active. State True to
turn on the feedback or False otherwise. Element 0 of the array defines
the horizontal plane as 1 and 2 define the vertical and longitudinal
plane respectively.
Attributes
----------
......@@ -31,111 +109,134 @@ class FeedBackTransverse(Element) :
y_pos : array
Retrieves and stores the position of the bunch, measured by the BPM,
in vertical plane.
tau_pos : array
Retrieves and stores the position of the bunch, measured by the BPM,
in longitudinal plane.
kick_x : array
Stores the numerical values of the kicks for the horizontal plane.
kick_y : array
Stores the numerical values of the kicks for the vertical plane.
kick_tau : array
Stores the numerical values of the kicks for the longitudinal plane.
coef_x : array
Stores the coefficients of the FIR filter for the horizontal plane.
coef_y : array
Stores the coefficients of the FIR filter for the vertical plane.
coef_tau : array
Stores the coefficients of the FIR filter for the longitudinal plane.
Methods
-------
get_fir(tap_number, tune, phase, turn_delay, gain)
Initialize the FIR filter for the desired plane and return an array
containing the FIR coefficients.
plot_fir(plane)
Plot the gain and the phase of the FIR filter selected.
References
----------
[1] T.Nakamura, S.Daté, K. Kobayashi, T. Ohshima. Proceedings of EPAC
2004. Transverse bunch by bunch feedback system for the Spring-8
storage ring.
"""
def __init__(self, ring, turn_delay=1, tap_number=5, ampl =1e-1, phase = -90):
def __init__(self, ring, tune, turn_delay, tap_number, gain, phase, bpm_error = np.zeros((3,)),
max_kick = np.zeros((3,)),plane = np.zeros((3,), dtype = bool)):
self.ring = ring
self.tune = tune
self.turn_delay = turn_delay
self.tap_number = tap_number
self.ampl = ampl
self.gain = gain
self.phase = phase
self.x_pos = np.zeros((self.tap_number,))
self.y_pos = np.zeros((self.tap_number,))
self.kick_x = np.zeros((self.turn_delay+1,))
self.kick_y = np.zeros((self.turn_delay+1,))
self.coef_x = self.get_fir_x()
self.coef_y = self.get_fir_y()
def get_fir_x(self):
self.bpm_error = bpm_error
self.max_kick = max_kick
self.plane = plane
self.x_pos = np.zeros((self.tap_number[0],))
self.y_pos = np.zeros((self.tap_number[1],))
self.tau_pos = np.zeros((self.tap_number[2],))
self.kick_x = np.zeros((self.turn_delay[0]+1,))
self.kick_y = np.zeros((self.turn_delay[1]+1,))
self.kick_tau = np.zeros((self.turn_delay[2]+1,))
self.coef_x = self.get_fir(self.tap_number[0], self.tune[0], self.phase[0], self.turn_delay[0], self.gain[0])
self.coef_y = self.get_fir(self.tap_number[1], self.tune[1], self.phase[1], self.turn_delay[1], self.gain[1])
self.coef_tau = self.get_fir(self.tap_number[2], self.tune[2], self.phase[2], self.turn_delay[2], self.gain[2])
def get_fir(self, tap_number, tune, phase, turn_delay, gain):
"""
Compute the FIR coefficients for the horizontal plane.
Compute the FIR coefficients for the selected plane(s).
This method is based on the FIR filter design algorithm developped by
T.Nakamura.
Returns
-------
coef_x : array
Array containing the FIR coefficients for the horizontal plane.
FIR_coef : array
Array containing the FIR coefficients for the selected plane(s).
"""
phi_x = 2*np.pi*self.ring.tune[0]
zeta = (self.phase*2*np.pi)/360
coef_x = np.zeros((self.tap_number,))
it = np.zeros((self.tap_number,))
for k in range(self.tap_number):
it[k] = (-k - self.turn_delay)
cs_x = np.zeros((self.tap_number,))
sn_x = np.zeros((self.tap_number,))
for i in range(self.tap_number):
cs_x[i] = (np.cos(phi_x*it[i]))
sn_x[i] = (np.sin(phi_x*it[i]))
CC_x = np.zeros((5, self.tap_number))
for j in range(self.tap_number):
CC_x[0][j] = 1
CC_x[1][j] = cs_x[j]
CC_x[2][j] = sn_x[j]
CC_x[3][j] = it[j]*sn_x[j]
CC_x[4][j] = it[j]*cs_x[j]
TCC_x = np.transpose(CC_x)
W_x = np.linalg.inv(CC_x.dot(TCC_x))
D_x = W_x.dot(CC_x)
for n in range(self.tap_number):
a_x = self.ampl*(D_x[1][n]*np.cos(zeta) + D_x[2][n])
coef_x[n] = (a_x)
return coef_x
def get_fir_y(self):
it = np.zeros((tap_number,))
CC = np.zeros((5, tap_number,))
zeta = (phase*2*np.pi)/360
for k in range(tap_number):
it[k] = (-k - turn_delay)
phi = 2*np.pi*tune
cs = np.cos(phi*it)
sn = np.sin(phi*it)
CC[0][:] = 1
CC[1][:] = cs
CC[2][:] = sn
CC[3][:] = it*sn
CC[4][:] = it*cs
TCC = np.transpose(CC)
W = np.linalg.inv(CC.dot(TCC))
D = W.dot(CC)
FIR_coef = gain*(D[1][:]*np.cos(zeta) + D[2][:]*np.sin(zeta))
return FIR_coef
def plot_fir(self, plane):
"""
Compute the FIR coefficients for the vertical plane.
Plot the gain and the phase of the FIR filter for the desired plane
Parameters
----------
plane : string
States "x" for the horizontal plane, "y" for the vertical and
"long" for the longitudinal one.
Returns
-------
coef_y : array
Array containing the FIR coefficients for the vertical plane.
None.
"""
phi_y = 2*np.pi*self.ring.tune[1]
zeta = (self.phase*2*np.pi)/360
coef_y = np.zeros((self.tap_number,))
it = np.zeros((self.tap_number,))
for k in range(0,self.tap_number):
it[k] = (-k - self.turn_delay)
cs_y = np.zeros((self.tap_number,))
sn_y = np.zeros((self.tap_number,))
for i in range(self.tap_number):
cs_y[i] = (np.cos(phi_y*it[i]))
sn_y[i] = (np.sin(phi_y*it[i]))
CC_y = np.zeros((5, self.tap_number))
for j in range(self.tap_number):
CC_y[0][j] = 1
CC_y[1][j] = cs_y[j]
CC_y[2][j] = sn_y[j]
CC_y[3][j] = it[j]*sn_y[j]
CC_y[4][j] = it[j]*cs_y[j]
TCC_y = np.transpose(CC_y)
W_y = np.linalg.inv(CC_y.dot(TCC_y))
D_y = W_y.dot(CC_y)
for n in range(self.tap_number):
a_y = self.ampl*(D_y[1][n]*np.cos(zeta) + D_y[2][n]*np.sin(zeta))
coef_y[n] = (a_y)
return coef_y
tune = np.arange(0,1,0.0001)
plane_dict = {"x":0 , "y":1 , "long":2}
index_plane = plane_dict[plane]
t = np.array([len(self.coef_x), len(self.coef_y), len(self.coef_tau)])
latency = np.exp(-1j*2*np.pi*tune*self.turn_delay[index_plane])
H_FIR = 0
liste = []
if(plane == "x"):
liste = self.coef_x
elif(plane == "y"):
liste = self.coef_y
elif(plane == "long"):
liste = self.coef_tau
for k in range(t[index_plane]):
H_FIR += liste[k]*np.exp(-1j*2*np.pi*(k)*tune)
H_tot = H_FIR * latency
gain = np.abs(H_tot)
phase = np.angle(H_tot, deg = True)
plt.figure("Gain"+str(index_plane))
plt.plot(tune, gain)
plt.title("Gain" + str(index_plane))
plt.xlabel("Tune")
plt.figure("Phase"+str(index_plane))
plt.plot(tune, phase)
plt.title("Phase in degree" + str(index_plane))
plt.xlabel("Tune")
plt.ylabel("Degree")
@Element.parallel
def track(self, bunch):
......@@ -148,23 +249,50 @@ class FeedBackTransverse(Element) :
----------
bunch : Bunch or Beam object
"""
self.x_pos[0] = bunch.mean[0]
self.y_pos[0] = bunch.mean[2]
if(self.plane[0] == True):
self.x_pos[0] = bunch.mean[0] + np.random.normal(0, self.bpm_error[0])
x_kick = 0
y_kick = 0
for k in range(self.tap_number):
for k in range(self.tap_number[0]):
x_kick += self.coef_x[k]*self.x_pos[k]
y_kick += self.coef_y[k]*self.y_pos[k]
if(self.max_kick[0] != 0 and x_kick > self.max_kick[0]):
x_kick = self.max_kick[0]
if(self.max_kick[0] != 0 and x_kick < -1*self.max_kick[0]):
x_kick = -1*self.max_kick[0]
self.kick_x[-1] = x_kick
self.kick_y[-1] = y_kick
bunch["xp"] -= self.kick_x[0]
bunch["yp"] -= self.kick_y[0]
bunch["xp"] += self.kick_x[0]
self.x_pos = np.roll(self.x_pos,1)
self.y_pos = np.roll(self.y_pos,1)
self.kick_x = np.roll(self.kick_x, -1)
if(self.plane[1] == True):
self.y_pos[0] = bunch.mean[2] + np.random.normal(0, self.bpm_error[1])
y_kick = 0
for k in range(self.tap_number[1]):
y_kick += self.coef_y[k]*self.y_pos[k]
if(self.max_kick[1] != 0 and y_kick > self.max_kick[1]):
y_kick = self.max_kick[1]
if(self.max_kick[1] != 0 and y_kick < -1*self.max_kick[1]):
y_kick = -1*self.max_kick[1]
self.kick_y[-1] = y_kick
bunch["yp"] += self.kick_y[0]
self.y_pos = np.roll(self.y_pos,1)
self.kick_y = np.roll(self.kick_y, -1)
if(self.plane[2] == True):
self.tau_pos[0] = bunch.mean[4] + np.random.normal(0, self.bpm_error[2])
tau_kick = 0
for k in range(self.tap_number[2]):
tau_kick += self.coef_tau[k]*self.tau_pos[k]
if(self.max_kick[2] != 0 and tau_kick > self.max_kick[2]):
tau_kick = self.max_kick[2]
if(self.max_kick[2] != 0 and tau_kick < -1*self.max_kick[2]):
tau_kick = -1*self.max_kick[2]
self.kick_tau[-1] = tau_kick
bunch["delta"] += self.kick_tau[0]
self.tau_pos = np.roll(self.tau_pos,1)
self.kick_tau = np.roll(self.kick_tau, -1)
\ No newline at end of file
mbtrack2/tracking/SimpleFeedBack.py
+ 10
− 9
View file @ 7d37efa0
......@@ -15,11 +15,12 @@ class SimpleFB(Element) :
Parameters
----------
ring : Synchrotron object
damping_time : array of shape (3,), optional
damping_time : array of shape (3,)
Contains the desired damping time in seconds for each plane in the same
order sepcified in plane.
phase_diff : float, optional
Phase of the damper kick in degrees.
phase_diff : array of shape (3,)
Contains the desired phase in degree for each plane in the same
order sepcified in plane.
plane : bool array of shape (3,), optional
Allow to choose on which plane the feedback is active. State True to
turn on the feedback or False otherwise. Element 0 of the array defines
......@@ -27,7 +28,7 @@ class SimpleFB(Element) :
respectively.
"""
def __init__(self, ring,phase_diff = 90, damping_time = np.array([2e-3,2e-5,2e-5]),plane = np.ones((3,), dtype= bool)):
def __init__(self, ring,phase_diff, damping_time, plane = np.zeros((3,), dtype= bool)):
self.ring = ring
self.damping_time = damping_time
self.phase_diff = phase_diff
......@@ -45,10 +46,10 @@ class SimpleFB(Element) :
bunch : Bunch or Beam object
"""
if(self.plane[0] == True):
bunch["delta"] -= (2*self.ring.T0/self.damping_time[0])*np.sin(self.phase_diff)*bunch.mean[5]
bunch["xp"] -= (2*self.ring.T0/self.damping_time[0])*np.sin(self.phase_diff[0])*bunch.mean[1]
if(self.plane[1] == True):
bunch["xp"] -= (2*self.ring.T0/self.damping_time[1])*np.sin(self.phase_diff)*bunch.mean[1]
bunch["yp"] -= (2*self.ring.T0/self.damping_time[1])*np.sin(self.phase_diff[1])*bunch.mean[3]
if(self.plane[2] == True):
bunch["yp"] -= (2*self.ring.T0/self.damping_time[2])*np.sin(self.phase_diff)*bunch.mean[3]
\ No newline at end of file
bunch["delta"] -= (2*self.ring.T0/self.damping_time[2])*np.sin(self.phase_diff[1])*bunch.mean[5]
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

The place to be for your source code