diff --git a/blond/impedances/impedance.py b/blond/impedances/impedance.py
index ff3fe1091b7980cda7600a11f1328e7423b9bd6a..05e4172065fd34a0e109596e9cd73059b6fffd50 100644
--- a/blond/impedances/impedance.py
+++ b/blond/impedances/impedance.py
@@ -357,7 +357,6 @@ class _InducedVoltage:
if self.multi_turn_wake:
# Number of points of the memory array for multi-turn wake
self.n_mtw_memory = self.n_induced_voltage
-
self.front_wake_buffer = 0
if self.mtw_mode == 'freq':
@@ -476,7 +475,7 @@ class _InducedVoltage:
# self.mtw_memory = bm.interp_const_space(self.time_mtw + t_rev,
self.mtw_memory = bm.interp(self.time_mtw + t_rev,
self.time_mtw, self.mtw_memory,
- left=0, right=0)
+ left=0, right=0) # todo
def _track(self):
"""
@@ -551,12 +550,15 @@ class InducedVoltageTime(_InducedVoltage):
####################################
# Call the __init__ method of the parent class [calls process()]
- _InducedVoltage.__init__(self, beam, profile,
- frequency_resolution=None,
- wake_length=wake_length,
- multi_turn_wake=multi_turn_wake,
- rf_station=rf_station, mtw_mode=mtw_mode,
- use_regular_fft=use_regular_fft)
+ super().__init__(beam=beam,
+ profile=profile,
+ frequency_resolution=None,
+ wake_length=wake_length,
+ multi_turn_wake=multi_turn_wake,
+ rf_station=rf_station,
+ mtw_mode=mtw_mode,
+ use_regular_fft=use_regular_fft,
+ )
def process(self) -> None:
"""
@@ -732,11 +734,11 @@ class InducedVoltageFreq(_InducedVoltage):
###############
# Call the __init__ method of the parent class
- _InducedVoltage.__init__(self, beam, profile, wake_length=None,
- frequency_resolution=frequency_resolution,
- multi_turn_wake=multi_turn_wake,
- rf_station=rf_station, mtw_mode=mtw_mode,
- use_regular_fft=use_regular_fft)
+ super().__init__(beam, profile, wake_length=None,
+ frequency_resolution=frequency_resolution,
+ multi_turn_wake=multi_turn_wake,
+ rf_station=rf_station, mtw_mode=mtw_mode,
+ use_regular_fft=use_regular_fft)
def process(self) -> None:
"""
@@ -868,7 +870,7 @@ class InductiveImpedance(_InducedVoltage):
self.deriv_mode = deriv_mode
# Call the __init__ method of the parent class
- _InducedVoltage.__init__(self, beam, profile, rf_station=rf_station)
+ super().__init__(beam, profile, rf_station=rf_station)
def induced_voltage_1turn(self, beam_spectrum_dict={}):
"""
@@ -924,7 +926,7 @@ class InducedVoltageResonator(_InducedVoltage):
solving the convolution integral with the resonator impedance analytically.
The line density need NOT be sampled at equidistant points. The times when
the induced voltage is calculated need to be the same where the line
- density is sampled. If no timeArray is passed, the induced voltage is
+ density is sampled. If no time_array is passed, the induced voltage is
evaluated at the points of the line density. This is necessary of
compatibility with other functions that calculate the induced voltage.
Currently, it requires the all quality factors :math:`Q>0.5`
@@ -942,6 +944,9 @@ class InducedVoltageResonator(_InducedVoltage):
Array of time values where the induced voltage is calculated.
If left out, the induced voltage is calculated at the times of the line
density.
+ array_length : int, optional
+ length of an array of one turn
+
Attributes
----------
@@ -949,15 +954,15 @@ class InducedVoltageResonator(_InducedVoltage):
Copy of the Beam object in order to access the beam info.
profile : Profile
Copy of the Profile object in order to access the line density.
- tArray : float array
+ timr_array : float array
Array of time values where the induced voltage is calculated.
If left out, the induced voltage is calculated at the times of the
line density
atLineDensityTimes : boolean
- flag indicating if the induced voltage has to be computed for timeArray
+ flag indicating if the induced voltage has to be computed for time_array
or for the line density
n_time : int
- length of tArray
+ length of time_array
R, omega_r, Q : lists of float
Resonators parameters
n_resonators : int
@@ -967,8 +972,16 @@ class InducedVoltageResonator(_InducedVoltage):
"""
@handle_legacy_kwargs
- def __init__(self, beam: Beam, profile: Profile, resonators: Resonators,
- time_array: Optional[NDArray] = None):
+ def __init__(self, beam: Beam,
+ profile: Profile,
+ resonators: Optional[Resonators] = None,
+ frequency_resolution: Optional[float] = None,
+ wake_length: Optional[float] = None,
+ multi_turn_wake: bool = False,
+ time_array: Optional[NDArray] = None,
+ rf_station: Optional[RFStation] = None,
+ array_length: Optional[int] = None,
+ use_regular_fft: bool = True) -> None:
# Test if one or more quality factors is smaller than 0.5.
if sum(resonators.Q < 0.5) > 0:
@@ -984,14 +997,15 @@ class InducedVoltageResonator(_InducedVoltage):
# If left out, the induced voltage is calculated at the times of the
# line density.
if time_array is None:
- self.tArray = self.profile.bin_centers
+ self.time_array = self.profile.bin_centers
self.atLineDensityTimes = True
else:
- self.tArray = time_array
+ self.time_array = time_array
self.atLineDensityTimes = False
- # Length of timeArray
- self.n_time = len(self.tArray)
+ self.array_length = array_length
+
+ self.n_time = len(self.time_array)
# Copy of the shunt impedances of the Resonators in* :math:`\Omega`
self.R = resonators.R_S
@@ -1001,7 +1015,6 @@ class InducedVoltageResonator(_InducedVoltage):
self.Q = resonators.Q
# Number of resonators
self.n_resonators = len(self.R)
-
# For internal use
self._Qtilde = self.Q * np.sqrt(1. - 1. / (4. * self.Q ** 2.))
self._reOmegaP = self.omega_r * self._Qtilde / self.Q
@@ -1016,15 +1029,17 @@ class InducedVoltageResonator(_InducedVoltage):
self._kappa1 = np.zeros(
int(self.profile.n_slices - 1), dtype=bm.precision.real_t, order='C')
- # Matrix to hold n_times many tArray[t]-bin_centers arrays.
+ # Matrix to hold n_times many time_array[t]-bin_centers arrays.
self._deltaT = np.zeros(
(self.n_time, self.profile.n_slices), dtype=bm.precision.real_t, order='C')
# Call the __init__ method of the parent class [calls process()]
- _InducedVoltage.__init__(self, beam, profile, wake_length=None,
- frequency_resolution=None,
- multi_turn_wake=False, rf_station=None,
- mtw_mode='time')
+ super().__init__(beam=beam, profile=profile,
+ frequency_resolution=frequency_resolution,
+ wake_length=wake_length,
+ multi_turn_wake=multi_turn_wake,
+ rf_station=rf_station, mtw_mode='time',
+ use_regular_fft=use_regular_fft)
def process(self):
r"""
@@ -1052,7 +1067,7 @@ class InducedVoltageResonator(_InducedVoltage):
self.profile.bin_centers,
self.profile.bin_size,
self.n_time, self._deltaT,
- self.tArray, self._reOmegaP,
+ self.time_array, self._reOmegaP,
self._imOmegaP, self._Qtilde,
self.n_resonators,
self.omega_r,
@@ -1063,6 +1078,23 @@ class InducedVoltageResonator(_InducedVoltage):
self.induced_voltage,
bm.precision.real_t))
+ def induced_voltage_mtw(self, beam_spectrum_dict={}):
+ r"""
+ Induced voltage method for InducedVoltageResonator.
+ mtw_memory is shifted by one turn, setting the values to 0.
+ The current turn's induced voltage is added to the memory of the previous turn.
+ Implementation by F. Batsch.
+ """
+ # shift the entries in array by 1 t_rev and set to 0
+ self.mtw_memory = np.append(self.mtw_memory, np.zeros(self.array_length))
+ # remove one turn length of memory
+ self.mtw_memory = self.mtw_memory[self.array_length:]
+ # Induced voltage of the current turn
+ self.induced_voltage_1turn(beam_spectrum_dict)
+ # Add induced voltage of the current turn to the memory from previous
+ self.mtw_memory[:int(self.n_time)] += self.induced_voltage
+ self.induced_voltage = self.mtw_memory[:self.n_time]
+
def to_gpu(self, recursive=True):
"""
Transfer all necessary arrays to the GPU
@@ -1075,7 +1107,7 @@ class InducedVoltageResonator(_InducedVoltage):
self.induced_voltage = cp.array(self.induced_voltage)
self._kappa1 = cp.array(self._kappa1)
self._deltaT = cp.array(self._deltaT)
- self.tArray = cp.array(self.tArray)
+ self.time_array = cp.array(self.time_array)
self._tmp_matrix = cp.array(self._tmp_matrix)
# to make sure it will not be called again
self._device: DeviceType = 'GPU'
@@ -1092,7 +1124,7 @@ class InducedVoltageResonator(_InducedVoltage):
self.induced_voltage = cp.asnumpy(self.induced_voltage)
self._kappa1 = cp.asnumpy(self._kappa1)
self._deltaT = cp.asnumpy(self._deltaT)
- self.tArray = cp.asnumpy(self.tArray)
+ self.time_array = cp.asnumpy(self.time_array)
self._tmp_matrix = cp.asnumpy(self._tmp_matrix)
# to make sure it will not be called again
diff --git a/blond/input_parameters/ring.py b/blond/input_parameters/ring.py
index 2322027e1941bb541153c8d4c1f66ab6a5ae90e5..f7137f0f7e70a5b40b85957247da423fbb1e8fd8 100644
--- a/blond/input_parameters/ring.py
+++ b/blond/input_parameters/ring.py
@@ -279,13 +279,20 @@ class Ring:
self.gamma: NDArray = np.sqrt(1 + (self.momentum / self.particle.mass) ** 2)
self.energy: NDArray = np.sqrt(self.momentum ** 2 + self.particle.mass ** 2)
self.kin_energy: NDArray = (np.sqrt(self.momentum ** 2 + self.particle.mass ** 2) -
- self.particle.mass)
- self.delta_E: NDArray = np.diff(self.energy, axis=1)
+ self.particle.mass)
self.t_rev: NDArray = np.dot(self.ring_length, 1 / (self.beta * c))
self.cycle_time: NDArray = np.cumsum(self.t_rev) # Always starts with zero
self.f_rev: NDArray = 1 / self.t_rev
self.omega_rev: NDArray = 2 * np.pi * self.f_rev
+ if self.n_sections == 1:
+ self.delta_E = np.diff(self.energy, axis=1)
+ else:
+ # when there is more than 1 RF station, self.energy has shape (n_sections, n_turns+1)
+ self.delta_E = np.zeros((n_sections, n_turns))
+ self.delta_E[0, :] = self.energy[0, 1:n_turns + 1] - self.energy[-1, 0:n_turns]
+ self.delta_E[1:, :] = self.energy[1:, 1:n_turns + 1] - self.energy[:-1, 1:n_turns + 1]
+
# Momentum compaction, checks, and derived slippage factors
if ring_options.t_start is None:
interp_time = self.cycle_time
diff --git a/unittests/impedances/data/mtw_resonator_induced_voltage.npy b/unittests/impedances/data/mtw_resonator_induced_voltage.npy
new file mode 100644
index 0000000000000000000000000000000000000000..f175f950e39d39e9b0ec0ef253f842baaaabbecb
Binary files /dev/null and b/unittests/impedances/data/mtw_resonator_induced_voltage.npy differ
diff --git a/unittests/impedances/test_impedance.py b/unittests/impedances/test_impedance.py
index ce73d65332f410d6425b4ab63bc6da281dd62821..e00e8c40abd94c8987531d2de5c273e91c2c2d85 100644
--- a/unittests/impedances/test_impedance.py
+++ b/unittests/impedances/test_impedance.py
@@ -16,16 +16,20 @@ Unittest for impedances.impedance
import unittest
import numpy as np
-
-from blond.beam.profile import CutOptions, Profile
-from blond.impedances.impedance import InducedVoltageFreq, InducedVoltageTime
+import os
+from blond.beam.profile import CutOptions, Profile, FitOptions
+from blond.impedances.impedance import InducedVoltageFreq, InducedVoltageTime, InducedVoltageResonator, \
+ TotalInducedVoltage
from blond.impedances.impedance_sources import Resonators
+from blond.beam.beam import Beam, Proton
+from blond.input_parameters.rf_parameters import RFStation
+from blond.input_parameters.ring import Ring
+from blond.trackers.tracker import RingAndRFTracker
class TestInducedVoltageFreq(unittest.TestCase):
def setUp(self):
-
# Nyquist frequency 1.6 GHz; frequency spacing 50 MHz
self.profile = Profile(None,
cut_options=CutOptions(cut_left=0, cut_right=5e-9, n_slices=16))
@@ -86,7 +90,6 @@ class TestInducedVoltageFreq(unittest.TestCase):
class TestInducedVoltageTime(unittest.TestCase):
def setUp(self):
-
# Nyquist frequency 1.6 GHz; frequency spacing 50 MHz
self.profile = Profile(None,
cut_options=CutOptions(cut_left=0, cut_right=5e-9, n_slices=16))
@@ -121,6 +124,73 @@ class TestInducedVoltageTime(unittest.TestCase):
np.testing.assert_allclose(test_object.wake_length_input, 11e-9)
-if __name__ == '__main__':
+class TestInducedVoltageResonatorMultiTurnWake(unittest.TestCase):
+ def setUp(self):
+
+ mom_compaction = 4.68e-4
+ circ = 5989.95
+ self.n_stations = 1
+ l_per_section = circ / self.n_stations
+ self.n_turns = 3
+ sectionlengths = np.full(self.n_stations, l_per_section)
+ alpha_c = np.full(self.n_stations, mom_compaction)
+ # program should be shape (n_stations, n_turns+1)
+ energy_program = np.array([1 * 25e9, 1 * 25e9, 1.04 * 25e9, 1.05e9])
+
+ self.ring = Ring(ring_length=sectionlengths, alpha_0=alpha_c,
+ particle=Proton(), n_turns=self.n_turns,
+ n_sections=self.n_stations, synchronous_data=energy_program)
+ self.R_shunt = 11897424000
+ self.Q_factor = 0.696e6
+ fdet = -1320
+ self.f_res = 1297263703
+ self.harmonic = int(self.f_res * self.ring.t_rev[0])
+ self.rf_station = RFStation(ring=self.ring, harmonic=[self.harmonic], voltage=[300e6], phi_rf_d=[0.0], n_rf=1)
+ self.n_slices = 2 ** 9
+ self.beam = Beam(ring=self.ring,
+ n_macroparticles=1001,
+ intensity=int(1e10))
+ cut_options = CutOptions(cut_left=0, cut_right=2 * np.pi, n_slices=self.n_slices,
+ cuts_unit='rad', rf_station=self.rf_station)
+ self.profile = Profile(self.beam, cut_options,
+ FitOptions(fit_option='gaussian'))
+
+ self.resonator = Resonators(self.R_shunt, self.f_res + fdet, self.Q_factor)
+
+ def test_total_induced_voltage(self):
+ potential_min_cav = self.rf_station.phi_s[0] / self.rf_station.omega_rf[0, 0]
+ min_index = np.abs(self.profile.bin_centers - potential_min_cav).argmin()
+ timeArray = []
+ for turn_ind in range(self.n_turns):
+ timeArray = np.append(timeArray,
+ self.rf_station.t_rev[turn_ind] * turn_ind +
+ np.linspace(self.profile.bin_centers[0],
+ self.profile.bin_centers[-1] + 2 * (
+ self.profile.bin_centers[min_index] -
+ self.profile.bin_centers[0]),
+ self.n_slices + 2 * min_index))
+ ivr = InducedVoltageResonator(beam=self.beam, profile=self.profile, time_array=timeArray,
+ multi_turn_wake=True,
+ resonators=self.resonator,
+ rf_station=self.rf_station,
+ wake_length=len(timeArray) * self.profile.bin_size,
+ array_length=int(len(timeArray) / self.n_turns))
+ totalvolt = TotalInducedVoltage(beam=self.beam, profile=self.profile, induced_voltage_list=[ivr])
+ longitudinal_tracker = RingAndRFTracker(rf_station=self.rf_station, beam=self.beam, profile=self.profile,
+ total_induced_voltage=totalvolt,
+ interpolation=True)
+ induced_voltages = []
+ for i in range(self.n_turns):
+ totalvolt.induced_voltage_sum()
+ longitudinal_tracker.track()
+ induced_voltages.append(longitudinal_tracker.totalInducedVoltage.induced_voltage)
+ induced_voltages = np.array(induced_voltages)
+
+ folder = os.path.abspath(os.path.dirname(__file__)) + "/data/"
+ expected_voltages = np.load(folder + "mtw_resonator_induced_voltage.npy")
+ np.testing.assert_allclose(induced_voltages, expected_voltages, rtol=1e-5, atol=1e-1)
+
+
+if __name__ == '__main__':
unittest.main()
diff --git a/unittests/input_parameters/test_ring.py b/unittests/input_parameters/test_ring.py
index 71863da90f1443bf3ad2dc177cd8f22127355a52..6289c21fc29ce4c382872911ed7d82296d3ceb11 100644
--- a/unittests/input_parameters/test_ring.py
+++ b/unittests/input_parameters/test_ring.py
@@ -191,6 +191,21 @@ class TestGeneralParameters(unittest.TestCase):
[np.float64(1.0), np.float64(2.0), np.float64(3.0)]) # might crash
Ring(C, alpha, momentum, Proton(), n_turns) # might crash
+ def test_multi_RF_dE(self):
+ # test if the delta_E is correctly calculated with multi-RF stations
+ C = 2 * np.pi * 1100.009 # Ring circumference [m]
+ l_per_section = C / self.num_sections
+ section_lengths = np.full(self.num_sections, l_per_section)
+ # the shape of momentum should be (n_sections, n_turns+1)
+ momentum = np.arange(0, (self.n_turns + 1) * self.num_sections).reshape(
+ (self.n_turns + 1, self.num_sections)).transpose() * 1e10
+ ring = Ring(ring_length=section_lengths, alpha_0=self.alpha_0,
+ synchronous_data=momentum, particle=self.particle,
+ n_turns=self.n_turns, n_sections=self.num_sections,
+ synchronous_data_type='momentum')
+ delta_E_values = ring.delta_E.flatten()
+ np.testing.assert_allclose(delta_E_values, delta_E_values[0], rtol=1e-3)
+
if __name__ == '__main__':
unittest.main()