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

@@ -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):