diff --git a/tracking/particles.py b/tracking/particles.py
index 1165ee0ca5024b94cfbcd78d0b73c7852d6404ac..feea6f4d95f08d07495e68cae2273328645275c6 100644
--- a/tracking/particles.py
+++ b/tracking/particles.py
@@ -692,7 +692,7 @@ class Beam:
self.mpi_switch = False
self.mpi = None
- def mpi_share_distributions(self, dimensions="tau", n_bins=75):
+ def mpi_share_distributions(self, dimensions="tau", n_bin=75):
"""
Compute the bunch profile and share it between the different bunches.
@@ -711,8 +711,8 @@ class Beam:
if isinstance(dimensions, str):
dimensions = [dimensions]
- if isinstance(n_bins, int):
- n_bins = np.ones((len(dimensions),), dtype=int)*n_bins
+ if isinstance(n_bin, int):
+ n_bin = np.ones((len(dimensions),), dtype=int)*n_bin
bunch = self[self.mpi.bunch_num]
@@ -722,8 +722,8 @@ class Beam:
for i in range(len(dimensions)):
dim = dimensions[i]
- n_bin = n_bins[i]
- bins, sorted_index, profile, center = bunch.binning(dimension=dim, n_bin=n_bin)
+ n = n_bin[i]
+ bins, sorted_index, profile, center = bunch.binning(dimension=dim, n_bin=n)
self.mpi.__setattr__(dim + "_center", np.empty((len(self), len(center)), dtype=np.float64))
self.mpi.comm.Allgather([center, MPI.DOUBLE], [self.mpi.__getattribute__(dim + "_center"), MPI.DOUBLE])
diff --git a/tracking/rf.py b/tracking/rf.py
index 8b553c5b223ff5abae612c83b1ad37cae5ba6a8c..49a743b1b821c65f60c228bbce1c6b98c199d274 100644
--- a/tracking/rf.py
+++ b/tracking/rf.py
@@ -196,6 +196,7 @@ class CavityResonator():
Track a Beam object through the CavityResonator object.
Can be used with or without mpi.
+ If used with mpi, beam.mpi_share_distributions must be called before.
The beam phasor is given at t=0 (synchronous particle) of the first
non empty bunch.
@@ -270,9 +271,10 @@ class CavityResonator():
self.nturn += 1
- def init_phasor(self, beam):
+ def init_phasor_track(self, beam):
"""
- Initialize the beam phasor for a given beam distribution.
+ Initialize the beam phasor for a given beam distribution using a
+ tracking like method.
Follow the same steps as the track method but in the "rf" reference
frame and without any modifications on the beam.
@@ -335,13 +337,14 @@ class CavityResonator():
if ref_frame == "beam":
delta = self.wr
elif ref_frame == "rf":
- delta = (self.wr - self.ring.omega1)
+ delta = (self.wr - self.m*self.ring.omega1)
self.beam_phasor = self.beam_phasor * np.exp((-1/self.filling_time +
1j*delta)*time)
def phasor_evol(self, profile, bin_length, charge_per_mp, ref_frame="beam"):
"""
- Compute the beam phasor evolution during the crossing of a bunch.
+ Compute the beam phasor evolution during the crossing of a bunch using
+ an analytic formula [1].
Assume that the phasor decay happens before the beam loading.
@@ -356,11 +359,15 @@ class CavityResonator():
ref_frame : string, optional
Reference frame to be used, can be "beam" or "rf".
+ References
+ ----------
+ [1] mbtrack2 manual.
+
"""
if ref_frame == "beam":
delta = self.wr
elif ref_frame == "rf":
- delta = (self.wr - self.ring.omega1)
+ delta = (self.wr - self.m*self.ring.omega1)
n_bin = len(profile)
@@ -376,7 +383,78 @@ class CavityResonator():
sum_val = -2 * sum_tot * charge_per_mp * self.loss_factor
self.beam_phasor += sum_val
-
+ def init_phasor(self, beam, n_bin=75):
+ """
+ Initialize the beam phasor for a given beam distribution using an
+ analytic formula [1].
+
+ No modifications on the Beam object.
+
+ Parameters
+ ----------
+ beam : Beam object
+ n_bin : int, optional
+ Number of bins to use. The default is 75.
+
+ References
+ ----------
+ [1] mbtrack2 manual.
+
+ """
+
+ # Initialization
+ if self.tracking is False:
+ self.init_tracking(beam)
+
+ N = n_bin-1
+ delta = (self.wr - self.m*self.ring.omega1)
+ n_turn = int(self.filling_time/self.ring.T0*10)
+
+ T = np.ones(self.ring.h)*self.ring.T1
+ bin_length = np.zeros(self.ring.h)
+ charge_per_mp = np.zeros(self.ring.h)
+ profile = np.zeros((N, self.ring.h))
+ center = np.zeros((N, self.ring.h))
+
+ # Gather beam distribution data
+ for j, bunch in enumerate(beam.not_empty):
+ index = self.valid_bunch_index[j]
+ if beam.mpi_switch:
+ beam.mpi_share_distributions(n_bin=n_bin)
+ center[:,index] = beam.mpi.tau_center[j]
+ profile[:,index] = beam.mpi.tau_profile[j]
+ bin_length[index] = center[1, index]-center[0, index]
+ charge_per_mp[index] = beam.mpi.charge_per_mp_all[j]
+ else:
+ (bins, sorted_index, profile[:, index], center[:, index]) = bunch.binning(n_bin=n_bin)
+ bin_length[index] = center[1, index]-center[0, index]
+ charge_per_mp[index] = bunch.charge_per_mp
+ T[index] -= (center[-1, index] + bin_length[index]/2)
+ if index != 0:
+ T[index - 1] += (center[0, index] - bin_length[index]/2)
+ T[self.ring.h - 1] += (center[0, 0] - bin_length[0]/2)
+
+ # Compute matrix coefficients
+ k = np.arange(0, N)
+ Tkj = np.zeros((N, self.ring.h))
+ for j in range(self.ring.h):
+ sum_t = np.array([T[n] + N*bin_length[n] for n in range(j+1,self.ring.h)])
+ Tkj[:,j] = (N-k)*bin_length[j] + T[j] + np.sum(sum_t)
+
+ var = np.exp( (-1/self.filling_time + 1j*delta) * Tkj )
+ sum_tot = np.sum((profile*charge_per_mp) * var)
+
+ # Use the formula n_turn times
+ for i in range(n_turn):
+ # Phasor decay during one turn
+ self.phasor_decay(self.ring.T0, ref_frame="rf")
+ # Phasor evolution due to induced voltage by marco-particles during one turn
+ sum_val = -2 * sum_tot * self.loss_factor
+ self.beam_phasor += sum_val
+
+ # Replace phasor at t=0 (synchronous particle) of the first non empty bunch.
+ idx0 = self.valid_bunch_index[0]
+ self.phasor_decay(center[-1,idx0] + bin_length[idx0]/2, ref_frame="rf")
@property
def generator_phasor(self):