Add ion critical mass and fixes

Fix mbi_threshold to have sigma in [s}
Add plot_critical_mass to instabilities (ions)
Add sigma(position) to synchrotron to compute the beam size at any ring location

collective_effects/instabilities.py

@@ -6,6 +6,8 @@ General calculations about instabilities
 @date: 15/01/2020
 """
 
+import numpy as np
+import matplotlib.pyplot as plt
 from scipy.constants import c, m_e, e, pi, epsilon_0
 
 def mbi_threshold(ring, sigma, R, b):
@@ -26,14 +28,15 @@ def mbi_threshold(ring, sigma, R, b):
     Returns
     -------
     I : float
-        MBI current threshold
+        MBI current threshold in [A]
         
     [1] : Y. Cai, "Theory of microwave instability and coherent synchrotron 
     radiation in electron storage rings", SLAC-PUB-14561
     [2] : D. Zhou, "Coherent synchrotron radiation and microwave instability in 
     electron storage rings", PhD thesis, p112
     """
-
+    
+    sigma = sigma * c
     Ia = 4*pi*epsilon_0*m_e*c**3/e # Alfven current
     chi = sigma*(R/b**3)**(1/2) # Shielding paramter
     xi = 0.5 + 0.34*chi
@@ -84,4 +87,57 @@ def cbi_threshold(ring, I, Vrf, f, beta, Ncav=1):
     Zydip = 1/(ring.tau[1]*beta[1]) * (2*ring.E0) / (ring.f0 * I * Ncav)
     
     return (Zlong, Zxdip, Zydip)
+
+def plot_critical_mass(ring, bunch_charge, bunch_spacing, n_points=1e4):
+    """
+    Plot ion critical mass, using Eq. (7.70) p147 of [1]
+
+    Parameters
+    ----------
+    ring : Synchrotron object
+    bunch_charge : float
+        Bunch charge in [C].
+    bunch_spacing : float
+        Time in between two adjacent bunches in [s].
+    n_points : float or int, optional
+        Number of point used in the plot. The default is 1e4.
+
+    Returns
+    -------
+    fig : figure
+        
+    References
+    ----------
+    [1] : Gamelin, A. (2018). Collective effects in a transient microbunching 
+    regime and ion cloud mitigation in ThomX (Doctoral dissertation, 
+    Université Paris-Saclay).
+
+    """
+    
+    n_points = int(n_points)
+    s = np.linspace(0, ring.L, n_points)
+    sigma = ring.sigma(s)
+    rp = 1.534698250004804e-18 # Proton classical radius, m
+    N = np.abs(bunch_charge/e)
+    
+    Ay = N*rp*bunch_spacing*c/(2*sigma[2,:]*(sigma[2,:] + sigma[0,:]))
+    Ax = N*rp*bunch_spacing*c/(2*sigma[0,:]*(sigma[2,:] + sigma[0,:]))
+    
+    fig = plt.figure()
+    ax = plt.gca()
+    ax.plot(s, Ax, label=r"$A_x^c$")
+    ax.plot(s, Ay, label=r"$A_y^c$")
+    ax.set_yscale("log")
+    ax.plot(s, np.ones_like(s)*2, label=r"$H_2^+$")
+    ax.plot(s, np.ones_like(s)*16, label=r"$H_2O^+$")
+    ax.plot(s, np.ones_like(s)*18, label=r"$CH_4^+$")
+    ax.plot(s, np.ones_like(s)*28, label=r"$CO^+$")
+    ax.plot(s, np.ones_like(s)*44, label=r"$CO_2^+$")
+    ax.legend()
+    ax.set_ylabel("Critical mass")
+    ax.set_xlabel("Longitudinal position [m]")
+    
+    return fig
+    
+    
     
\ No newline at end of file

tracking/synchrotron.py

@@ -223,18 +223,47 @@ class Synchrotron:
         """Momentum compaction"""
         return self.ac - 1/(self.gamma**2)
     
-    @property
-    def sigma(self):
-        """RMS beam size at equilibrium"""
-        sigma = np.zeros((4,))
-        sigma[0] = (self.emit[0]*self.optics.local_beta[0] +
-                    self.optics.local_dispersion[0]**2*self.sigma_delta)**0.5
-        sigma[1] = (self.emit[0]*self.optics.local_alpha[0] +
-                    self.optics.local_dispersion[1]**2*self.sigma_delta)**0.5
-        sigma[2] = (self.emit[1]*self.optics.local_beta[1] +
-                    self.optics.local_dispersion[2]**2*self.sigma_delta)**0.5
-        sigma[3] = (self.emit[1]*self.optics.local_alpha[1] +
-                    self.optics.local_dispersion[3]**2*self.sigma_delta)**0.5
+    def sigma(self, position=None):
+        """
+        RMS beam size at equilibrium
+
+        Parameters
+        ----------
+        position : float or array, optional
+            Longitudinal position where the beam size is computed. If None, 
+            the local values are used.
+
+        Returns
+        -------
+        sigma : array
+            RMS beam size at position location or at local positon if position
+            is None.
+
+        """
+        if position is None:
+            sigma = np.zeros((4,))
+            sigma[0] = (self.emit[0]*self.optics.local_beta[0] +
+                        self.optics.local_dispersion[0]**2*self.sigma_delta)**0.5
+            sigma[1] = (self.emit[0]*self.optics.local_alpha[0] +
+                        self.optics.local_dispersion[1]**2*self.sigma_delta)**0.5
+            sigma[2] = (self.emit[1]*self.optics.local_beta[1] +
+                        self.optics.local_dispersion[2]**2*self.sigma_delta)**0.5
+            sigma[3] = (self.emit[1]*self.optics.local_alpha[1] +
+                        self.optics.local_dispersion[3]**2*self.sigma_delta)**0.5
+        else:
+            if isinstance(position, (float, int)):
+                n = 1
+            else:
+                n = len(position)
+            sigma = np.zeros((4, n))
+            sigma[0,:] = (self.emit[0]*self.optics.beta(position)[0] +
+                        self.optics.dispersion(position)[0]**2*self.sigma_delta)**0.5
+            sigma[1,:] = (self.emit[0]*self.optics.alpha(position)[0] +
+                        self.optics.dispersion(position)[1]**2*self.sigma_delta)**0.5
+            sigma[2,:] = (self.emit[1]*self.optics.beta(position)[1] +
+                        self.optics.dispersion(position)[2]**2*self.sigma_delta)**0.5
+            sigma[3,:] = (self.emit[1]*self.optics.alpha(position)[1] +
+                        self.optics.dispersion(position)[3]**2*self.sigma_delta)**0.5
         return sigma
     
     def synchrotron_tune(self, Vrf):