From d5848c30a44da2a621682fb53cb2588888e6488a Mon Sep 17 00:00:00 2001
From: Gamelin Alexis <gamelin@synchrotron-soleil.fr>
Date: Wed, 28 Apr 2021 17:38:31 +0200
Subject: [PATCH] Add phasor initialization using analytic formula and fix
phase slipping factor for init
CavityResonator.init_phasor is now CavityResonator.init_phasor_track
New CavityResonator.init_phasor uses analytic formula
Fix phase slipping factor for initialization case (rf frame should be delta = (self.wr - self.m*self.ring.omega1))
fix init
---
tracking/rf.py | 90 ++++++++++++++++++++++++++++++++++++++++++++++----
1 file changed, 84 insertions(+), 6 deletions(-)
diff --git a/tracking/rf.py b/tracking/rf.py
index 8b553c5..49a743b 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):
--
GitLab