Change beam phasor evolution + phasor init

Add phasor initialization (on the "rf" reference)
Change beam phasor evolution to allow for "beam" or "rf" reference

tracking/rf.py

@@ -188,7 +188,7 @@ class CavityResonator():
         self.valid_bunch_index = np.where(beam.filling_pattern == True)[0]
         self.tracking = True
     
-    def track(self, beam):
+    def track(self, beam, ref_frame="beam"):
         """version type mbtrack + formule analytique"""
         
         if self.tracking is False:
@@ -210,7 +210,7 @@ class CavityResonator():
             
             # phasor is given at t=0 synchronous particle of the bunch
             # remove part of beam phasor decay to be at the start of the binning
-            self.phasor_decay(center[0])
+            self.phasor_decay(center[0], ref_frame=ref_frame)
             
             if index != self.bunch_index:
                 self.phasor_evol(profile, bin_length, charge_per_mp)
@@ -229,17 +229,43 @@ class CavityResonator():
                     Vtot = Vgene + Vbeam - charge_per_mp*self.loss_factor*mp_per_bin
                     energy_change[ind] = Vtot / self.ring.E0
 
-                    self.phasor_decay(bin_length)
-                    
                     self.beam_phasor -= 2*charge_per_mp*self.loss_factor*mp_per_bin
+                    self.phasor_decay(bin_length, ref_frame=ref_frame)
                 
-            self.phasor_decay( (self.distance[index] * self.ring.T1) - center[-1])
+            self.phasor_decay( (self.distance[index] * self.ring.T1) - center[-1], ref_frame=ref_frame)
             
             if index == self.bunch_index:
                 # apply kick
                 bunch["delta"] += energy_change
+                bunch.energy_change = energy_change
+                
+        
+    def init_phasor(self, beam):
+        """version type mbtrack + formule analytique"""
+        
+        if self.tracking is False:
+            self.init_tracking(beam)
+            
+        n_turn = int(self.filling_time/self.ring.T0*5)
         
-    def phasor_decay(self, time):
+        for i in range(n_turn):
+            for j, bunch in enumerate(beam.not_empty):
+                
+                index = self.valid_bunch_index[j]
+                
+                # get shared bunch profile for current bunch
+                center = beam.mpi.center_all[j]
+                profile = beam.mpi.profile_all[j]
+                bin_length = center[1]-center[0]
+                charge_per_mp = beam.mpi.charge_per_mp_all[j]
+                
+                # phasor is given at t=0 synchronous particle of the bunch
+                # remove part of beam phasor decay to be at the start of the binning
+                self.phasor_decay(center[0], ref_frame="rf")
+                self.phasor_evol(profile, bin_length, charge_per_mp, ref_frame="rf")
+                self.phasor_decay( (self.distance[index] * self.ring.T1) - center[-1], ref_frame="rf")
+            
+    def phasor_decay(self, time, ref_frame="beam"):
         """
         Compute the beam phasor decay during a given time span, assuming that 
         no particles are crossing the cavity during the time span.
@@ -250,10 +276,14 @@ class CavityResonator():
             Time span in [s], can be positive or negative.
 
         """
+        if ref_frame == "beam":
+            delta = self.wr
+        elif ref_frame == "rf":
+            delta = (self.wr - self.ring.omega1)
         self.beam_phasor = self.beam_phasor * np.exp((-1/self.filling_time +
-                                  1j*2*np.pi*self.detune)*time)
+                                  1j*delta)*time)
         
-    def phasor_evol(self, profile, bin_length, charge_per_mp):
+    def phasor_evol(self, profile, bin_length, charge_per_mp, ref_frame="beam"):
         """
         Compute the beam phasor evolution during the crossing of a bunch.
 
@@ -266,15 +296,20 @@ class CavityResonator():
         charge_per_mp : float
             Charge per macro-particle in [C].
         """
+        if ref_frame == "beam":
+            delta = self.wr
+        elif ref_frame == "rf":
+            delta = (self.wr - self.ring.omega1)
+            
         n_bin = len(profile)
         
         # Phasor decay during crossing time
         deltaT = n_bin*bin_length
-        self.phasor_decay(deltaT)
+        self.phasor_decay(deltaT, ref_frame)
         
         # Phasor evolution due to induced voltage by marco-particles
         k = np.arange(1, n_bin + 1)
-        var = np.exp( (-1/self.filling_time + 1j*2*np.pi*self.detune) * 
+        var = np.exp( (-1/self.filling_time + 1j*delta) * 
                      (n_bin-k) * deltaT )
         sum_tot = np.sum(profile * var)
         sum_val = -2 * sum_tot * charge_per_mp * self.loss_factor
@@ -297,11 +332,28 @@ class CavityResonator():
         """Cavity total voltage in [V]"""
         return np.abs(self.cavity_phasor)
     
+    def cavity_real_voltage(self, t):
+        """Cavity total voltage in [V]"""
+        Vg = self.Vg*np.exp(1j*(self.m*self.ring.omega1*t + self.theta_g))
+        Vb = self.beam_phasor*np.exp((-1/self.filling_time +
+                                  1j*2*np.pi*self.detune)*t)
+        return np.real(Vg+Vb)
+    
     @property
     def cavity_phase(self):
         """Cavity total phase in [rad]"""
         return np.angle(self.cavity_phasor)
     
+    @property
+    def beam_voltage(self):
+        """Beam loading voltage in [V]"""
+        return np.abs(self.beam_phasor)
+    
+    @property
+    def beam_phase(self):
+        """Beam loading phase in [rad]"""
+        return np.angle(self.beam_phasor)
+    
     @property
     def loss_factor(self):
         """Cavity loss factor in [V/C]"""