diff --git a/collective_effects/instabilities.py b/collective_effects/instabilities.py
index d33e3047858fe89c979638859e95c15413154f75..12598cae3576dde196f51d92b812735be6285bc8 100644
--- a/collective_effects/instabilities.py
+++ b/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
diff --git a/tracking/synchrotron.py b/tracking/synchrotron.py
index e3df7cd4fac3d68571d1f57c6aa0fee073ae8a6e..d5f3107b2cdf0aa1097bf60d868911c7a360f413 100644
--- a/tracking/synchrotron.py
+++ b/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):