- Added Chebyshev, Legenre and Sacherer modes

for spectral_density()
- Corrected Hermite mode to include normalisation by the mode number
- some authorefactoring

mbtrack2/utilities/spectrum.py

@@ -4,8 +4,10 @@ Module where bunch and beam spectrums and profile are defined.
 """
 
 import numpy as np
+from scipy.special import jv, spherical_jn
 
-def spectral_density(frequency, sigma, m = 1, mode="Hermite"):
+
+def spectral_density(frequency, sigma, m=1, k=0, mode="Hermite"):
     """
     Compute the spectral density of different modes for various values of the
     head-tail mode number, based on Table 1 p238 of [1].
@@ -18,9 +20,14 @@ def spectral_density(frequency, sigma, m = 1, mode="Hermite"):
         RMS bunch length in [s]
     m : int, optional
         head-tail (or azimutal/synchrotron) mode number
+    k : int, optional
+        radial mode number (such that |q|=m+2k, where |q| is the head-tail mode number)
     mode: str, optional
         type of the mode taken into account for the computation:
-        -"Hermite" modes for Gaussian bunches
+        -"Hermite" modes for Gaussian bunches (typical for electrons)
+        -"Chebyshev" for airbag bunches
+        -"Legendre" for parabolic bunches (typical for protons)
+        -"Sacherer" or "Sinusoidal" simplifying approximation of Legendre modes from [3]
 
     Returns
     -------
@@ -29,14 +36,26 @@ def spectral_density(frequency, sigma, m = 1, mode="Hermite"):
     References
     ----------
     [1] : Handbook of accelerator physics and engineering, 3rd printing.
+    [2] : Ng, K. Y. (2005). Physics of intensity dependent beam instabilities. WORLD SCIENTIFIC. https://doi.org/10.1142/5835
+    [3] : Sacherer, F. J. (1972). Methods for computing bunched beam instabilities. CERN Tech. rep. CERN/SI-BR/72-5 https://cds.cern.ch/record/2291670?ln=en
     """
 
     if mode == "Hermite":
-        return (2*np.pi*frequency*sigma)**(2*m)*np.exp(
-                -1*(2*np.pi*frequency*sigma)**2)
+        return 1/(np.math.factorial(m)*2**m)*(2*np.pi*frequency*sigma)**(2*m)*np.exp(
+            -(2*np.pi*frequency*sigma)**2)
+    elif mode == "Chebyshev":
+        tau_l = 4*sigma
+        return (jv(m, 2*np.pi*frequency*tau_l))**2
+    elif mode == "Legendre":
+        tau_l = 4*sigma
+        return (spherical_jn(m,  np.abs(2*np.pi*frequency*tau_l)))**2
+    elif mode == "Sacherer" or mode == "Sinusoidal":
+        y = 4*2*np.pi*frequency*sigma/np.pi
+        return (2*(m+1)/np.pi*1/np.abs(y**2-(m+1)**2)*np.sqrt(1+(-1)**m*np.cos(np.pi*y)))**2
     else:
         raise NotImplementedError("Not implemanted yet.")
 
+
 def gaussian_bunch_spectrum(frequency, sigma):
     """
     Compute a Gaussian bunch spectrum [1].
@@ -60,6 +79,7 @@ def gaussian_bunch_spectrum(frequency, sigma):
     """
     return np.exp(-1/2*(2*np.pi*frequency)**2*sigma**2)
 
+
 def gaussian_bunch(time, sigma):
     """
     Compute a Gaussian bunch profile.