diff --git a/README.md b/README.md
deleted file mode 100644
index a33e1c283e3c898242e67be804cb883ffed2b624..0000000000000000000000000000000000000000
--- a/README.md
+++ /dev/null
@@ -1,2 +0,0 @@
-# mbtrack2
-
diff --git a/SOLEIL_U_74BA_HOA_SYM02_V0313_cleaned.mat b/SOLEIL_U_74BA_HOA_SYM02_V0313_cleaned.mat
deleted file mode 100644
index f021674c7c458deeff2a876be27f7853ba1a1b3b..0000000000000000000000000000000000000000
Binary files a/SOLEIL_U_74BA_HOA_SYM02_V0313_cleaned.mat and /dev/null differ
diff --git a/collective_effects/__init__.py b/collective_effects/__init__.py
index 27c993fd532ba8a41fdb6ab902ed422559965666..867d0d105cb7a1ffa8141c96eb0c45f72a3d7975 100644
--- a/collective_effects/__init__.py
+++ b/collective_effects/__init__.py
@@ -5,3 +5,10 @@ Created on Tue Jan 14 12:25:30 2020
@author: gamelina
"""
+from mbtrack2.collective_effects.instabilities import MBI_threshold
+from mbtrack2.collective_effects.resistive_wall import skin_depth, CircularResistiveWall
+from mbtrack2.collective_effects.tapers import StupakovRectangularTaper, StupakovCircularTaper
+from mbtrack2.collective_effects.wakefield import ComplexData, Impedance, WakeFunction, ImpedanceModel, Wakefield
+from mbtrack2.collective_effects.wakefield import read_CST, read_IW2D, spectral_density
+from mbtrack2.collective_effects.wakefield import gaussian_bunch_spectrum, beam_spectrum, effective_impedance
+from mbtrack2.collective_effects.wakefield import yokoya_elliptic, beam_loss_factor
\ No newline at end of file
diff --git a/collective_effects/resistive_wall.py b/collective_effects/resistive_wall.py
index 95ae59f68698a0d0bc8f51374bd0b33d66493714..240e30e268dfa956dd1a2dbc005e88cfe6897dca 100644
--- a/collective_effects/resistive_wall.py
+++ b/collective_effects/resistive_wall.py
@@ -8,7 +8,7 @@ Define resistive wall elements based on the Wakefield class.
import numpy as np
from scipy.constants import mu_0, epsilon_0, c
-from collective_effects.wakefield import Wakefield, Impedance
+from mbtrack2.collective_effects.wakefield import Wakefield, Impedance
def skin_depth(omega, rho, mu_r = 1, epsilon_r = 1):
""" General formula for the skin depth
diff --git a/collective_effects/tapers.py b/collective_effects/tapers.py
index 0e09b36bd476ef0d97915b6d92b861bc800375d2..34dde1855a4ab3b3b2fcb0f220da80010286ab3b 100644
--- a/collective_effects/tapers.py
+++ b/collective_effects/tapers.py
@@ -8,7 +8,7 @@ Created on Tue Mar 31 13:15:36 2020
from scipy.constants import mu_0, epsilon_0, c, pi
import numpy as np
from scipy.integrate import trapz
-from collective_effects.wakefield import Wakefield, Impedance
+from mbtrack2.collective_effects.wakefield import Wakefield, Impedance
class StupakovRectangularTaper(Wakefield):
"""
diff --git a/collective_effects/tools.py b/collective_effects/tools.py
deleted file mode 100644
index 44c0efc84bb475d8d23b34794ab52449bbcb2335..0000000000000000000000000000000000000000
--- a/collective_effects/tools.py
+++ /dev/null
@@ -1,412 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Tools and utilities functions for the Wakefield and Impedances classes.
-
-@author: Alexis Gamelin
-@date: 24/04/2020
-"""
-
-import pandas as pd
-import numpy as np
-from scipy.integrate import trapz
-from scipy.interpolate import interp1d
-from collective_effects.wakefield import Wakefield, Impedance, WakeFunction
-
-def read_CST(file, impedance_type='long', divide_by=None):
- """
- Read CST file format into an Impedance object.
-
- Parameters
- ----------
- file : str
- Path to the file to read.
- impedance_type : str, optional
- Type of the Impedance object.
- divide_by : float, optional
- Divide the impedance by a value. Mainly used to normalize transverse
- impedance by displacement.
-
- Returns
- -------
- result : Impedance object
- Data from file.
- """
- df = pd.read_csv(file, comment="#", header = None, sep = "\t",
- names = ["Frequency","Real","Imaginary"])
- df["Frequency"] = df["Frequency"]*1e9
- if divide_by is not None:
- df["Real"] = df["Real"]/divide_by
- df["Imaginary"] = df["Imaginary"]/divide_by
- if impedance_type == "long":
- df["Real"] = np.abs(df["Real"])
- df.set_index("Frequency", inplace = True)
- result = Impedance(variable = df.index,
- function = df["Real"] + 1j*df["Imaginary"],
- impedance_type=impedance_type)
- return result
-
-def read_IW2D(file, impedance_type='long'):
- """
- Read IW2D file format into an Impedance object.
-
- Parameters
- ----------
- file : str
- Path to the file to read.
- impedance_type : str, optional
- Type of the Impedance object.
-
- Returns
- -------
- result : Impedance object
- Data from file.
- """
- df = pd.read_csv(file, sep = ' ', header = None,
- names = ["Frequency","Real","Imaginary"], skiprows=1)
- df.set_index("Frequency", inplace = True)
- result = Impedance(variable = df.index,
- function = df["Real"] + 1j*df["Imaginary"],
- impedance_type=impedance_type)
- return result
-
-def loss_factor(wake, sigma):
- """
- Compute the loss factor or the kick factor of a Wakefield, Impedance or
- WakeFunction for a Gaussian bunch. Formulas from Eq. (2) p239 and
- Eq.(7) p240 of [1].
-
- Parameters
- ----------
- wake : Wakefield, Impedance or WakeFunction object.
- sigma : float
- RMS bunch length in [s]
-
- Returns
- -------
- kloss: float
- Loss factor in [V/C] or kick factor in [V/C/m] depanding on the
- impedance type.
-
- References
- ----------
- [1] : Handbook of accelerator physics and engineering, 3rd printing.
- """
-
- if isinstance(wake, Impedance):
-
- positive_index = wake.data.index > 0
- frequency = wake.data.index[positive_index]
- sd = spectral_density(frequency, sigma, m=0)
-
- if(wake.impedance_type == "long"):
- kloss = trapz(x = frequency,
- y = 2*wake.data["real"][positive_index]*sd)
- elif(wake.impedance_type == "xdip" or wake.impedance_type == "ydip"):
- kloss = trapz(x = frequency,
- y = 2*wake.data["imag"][positive_index]*sd)
- elif(wake.impedance_type == "xquad" or wake.impedance_type == "yquad"):
- kloss = trapz(x = frequency,
- y = 2*wake.data["imag"][positive_index]*sd)
- else:
- raise TypeError("Impedance type not recognized.")
-
- if isinstance(wake, WakeFunction):
- # To be implemented.
- pass
-
- if isinstance(wake, Wakefield):
- # To be implemented. Create field in wakefield with values ?
- pass
-
- return kloss
-
-def spectral_density(frequency, sigma, m = 1, 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].
-
- Parameters
- ----------
- frequency : list or numpy array
- sample points of the spectral density in [Hz]
- sigma : float
- RMS bunch length in [s]
- m : int, optional
- head-tail (or azimutal/synchrotron) mode number
- mode: str, optional
- type of the mode taken into account for the computation:
- -"Hermite" modes for Gaussian bunches
-
- Returns
- -------
- numpy array
-
- References
- ----------
- [1] : Handbook of accelerator physics and engineering, 3rd printing.
- """
-
- if mode == "Hermite":
- return (2*np.pi*frequency*sigma)**(2*m)*np.exp(
- -1*(2*np.pi*frequency*sigma)**2)
- else:
- raise NotImplementedError("Not implemanted yet.")
-
-def Gaussian_bunch_spectrum(frequency, sigma):
- """
- Compute a Gaussian bunch spectrum [1].
-
- Parameters
- ----------
- frequency : array
- sample points of the beam spectrum in [Hz].
- sigma : float
- RMS bunch length in [s].
-
- Returns
- -------
- bunch_spectrum : array
- Bunch spectrum sampled at points given in frequency.
-
- References
- ----------
- [1] : Gamelin, A. (2018). Collective effects in a transient microbunching
- regime and ion cloud mitigation in ThomX. p86, Eq. 4.19
- """
- return np.exp(-1/2*(2*np.pi*frequency)**2*sigma**2)
-
-
-def beam_spectrum(frequency, M, bunch_spacing, sigma=None,
- bunch_spectrum=None):
- """
- Compute the beam spectrum assuming constant spacing between bunches [1].
-
- Parameters
- ----------
- frequency : list or numpy array
- sample points of the beam spectrum in [Hz].
- M : int
- Number of bunches.
- bunch_spacing : float
- Time between two bunches in [s].
- sigma : float, optional
- If bunch_spectrum is None then a Gaussian bunch with sigma RMS bunch
- length in [s] is assumed.
- bunch_spectrum : array, optional
- Bunch spectrum sampled at points given in frequency.
-
- Returns
- -------
- beam_spectrum : array
-
- References
- ----------
- [1] Rumolo, G - Beam Instabilities - CAS - CERN Accelerator School:
- Advanced Accelerator Physics Course - 2014, Eq. 9
- """
-
- if bunch_spectrum is None:
- bunch_spectrum = Gaussian_bunch_spectrum(frequency, sigma)
-
- beam_spectrum = (bunch_spectrum * np.exp(1j*np.pi*frequency *
- bunch_spacing*(M-1)) *
- np.sin(M*np.pi*frequency*bunch_spacing) /
- np.sin(np.pi*frequency*bunch_spacing))
-
- return beam_spectrum
-
-
-def effective_impedance(ring, imp, m, mu, sigma, xi=None, mode="Hermite"):
- """
- Compute the effective (longitudinal or transverse) impedance.
- Formulas from Eq. (1) and (2) p238 of [1].
-
- Parameters
- ----------
- ring : Synchrotron object
- imp : Impedance object
- mu : int
- coupled bunch mode number, goes from 0 to (M-1) where M is the
- number of bunches
- m : int
- head-tail (or azimutal/synchrotron) mode number
- sigma : float
- RMS bunch length in [s]
- xi : float, optional
- (non-normalized) chromaticity
- mode: str, optional
- type of the mode taken into account for the computation:
- -"Hermite" modes for Gaussian bunches
-
- Returns
- -------
- Zeff : float
- effective impedance in [ohm] or in [ohm/m] depanding on the impedance
- type.
-
- [1] Handbook of accelerator physics and engineering, 3rd printing.
- """
-
- if not isinstance(imp, Impedance):
- raise TypeError("{} should be an Impedance object.".format(imp))
-
- fmin = imp.data.index.min()
- fmax = imp.data.index.max()
- if fmin > 0:
- raise ValueError("A double sided impedance spectrum, with positive and"
- " negative frequencies, is needed to compute the "
- "effective impedance.")
-
- if mode == "Hermite":
- def h(f):
- return spectral_density(frequency=f, sigma=sigma, m=m,
- mode="Hermite")
- else:
- raise NotImplementedError("Not implemanted yet.")
-
- M = 416
-
- pmax = fmax/(ring.f0 * M) - 50
- pmin = fmin/(ring.f0 * M) + 50
-
- p = np.arange(pmin,pmax+1)
-
- if imp.impedance_type == "long":
- fp = ring.f0*(p*M + mu + m*ring.tune[2])
- fp = fp[np.nonzero(fp)] # Avoid division by 0
- num = np.sum( imp(fp) * h(fp) / (fp*2*np.pi) )
- den = np.sum( h(fp) )
- Zeff = num/den
-
- elif imp.impedance_type == "xdip" or imp.impedance_type == "ydip":
- if imp.impedance_type == "xdip":
- tuneXY = ring.tune[0]
- xi = ring.chro[0]
- elif imp.impedance_type == "ydip":
- tuneXY = ring.tune[1]
- xi = ring.chro[1]
- fp = ring.f0*(p*M + mu + tuneXY + m*ring.tuneS)
- f_xi = xi/ring.eta*ring.f0
- num = np.sum( imp(fp) * h(fp - f_xi) )
- den = np.sum( h(fp - f_xi) )
- Zeff = num/den
- else:
- raise TypeError("Effective impedance is only defined for long, xdip"
- " and ydip impedance type.")
-
- return Zeff
-
-
-def Yokoya_elliptic(x_radius , y_radius):
- """
- Compute Yokoya factors for an elliptic beam pipe.
- Function adapted from N. Mounet IW2D.
-
- Parameters
- ----------
- x_radius : float
- Horizontal semi-axis of the ellipse in [m].
- y_radius : float
- Vertical semi-axis of the ellipse in [m].
-
- Returns
- -------
- yoklong : float
- Yokoya factor for the longitudinal impedance.
- yokxdip : float
- Yokoya factor for the dipolar horizontal impedance.
- yokydip : float
- Yokoya factor for the dipolar vertical impedance.
- yokxquad : float
- Yokoya factor for the quadrupolar horizontal impedance.
- yokyquad : float
- Yokoya factor for the quadrupolar vertical impedance.
- """
- if y_radius < x_radius:
- small_semiaxis = y_radius
- large_semiaxis = x_radius
- else:
- small_semiaxis = x_radius
- large_semiaxis = y_radius
-
- # read Yokoya factors interpolation file
- # BEWARE: columns are ratio, dipy, dipx, quady, quadx
- yokoya_file = pd.read_csv("collective_effects/data/" +
- "Yokoya_elliptic_from_Elias_USPAS.csv")
- ratio_col = yokoya_file["x"]
- # compute semi-axes ratio (first column of this file)
- ratio = (large_semiaxis - small_semiaxis)/(large_semiaxis + small_semiaxis)
-
- # interpolate Yokoya file at the correct ratio
- yoklong = 1
-
- if y_radius < x_radius:
- yokydip = np.interp(ratio, ratio_col, yokoya_file["dipy"])
- yokxdip = np.interp(ratio, ratio_col, yokoya_file["dipx"])
- yokyquad = np.interp(ratio, ratio_col, yokoya_file["quady"])
- yokxquad = np.interp(ratio, ratio_col, yokoya_file["quadx"])
- else:
- yokxdip = np.interp(ratio, ratio_col, yokoya_file["dipy"])
- yokydip = np.interp(ratio, ratio_col, yokoya_file["dipx"])
- yokxquad = np.interp(ratio, ratio_col, yokoya_file["quady"])
- yokyquad = np.interp(ratio, ratio_col, yokoya_file["quadx"])
-
- return (yoklong, yokxdip, yokydip, yokxquad, yokyquad)
-
-def beam_loss_factor(impedance, frequency, spectrum, ring):
- """
- Compute "beam" loss factor using the beam spectrum, uses a sum instead of
- integral compared to loss_factor [1].
-
- Parameters
- ----------
- impedance : Impedance of type "long"
- frequency : array
- Sample points of spectrum.
- spectrum : array
- Beam spectrum to consider.
- ring : Synchrotron object
-
- Returns
- -------
- kloss_beam : float
- Beam loss factor in [V/C].
-
- References
- ----------
- [1] : Handbook of accelerator physics and engineering, 3rd printing.
- Eq (3) p239.
- """
- pmax = np.floor(impedance.data.index.max()/ring.f0)
- pmin = np.floor(impedance.data.index.min()/ring.f0)
-
- if pmin >= 0:
- # Add negative part
- negative_index = impedance.data.index*-1
- negative_data = impedance.data.set_index(negative_index)
- negative_data["imag"] = -1*negative_data["imag"]
- try:
- negative_data = negative_data.drop(0)
- except KeyError:
- pass
-
- all_data = impedance.data.append(negative_data)
- all_data = all_data.sort_index()
- impedance.data = all_data
-
- pmin = -1*pmax
-
- p = np.arange(pmin+1,pmax)
- pf0 = p*ring.f0
- ReZ = np.real(impedance(pf0))
- spectral_density = np.abs(spectrum)**2
- # interpolation of the spectrum is needed to avoid problems liked to
- # division by 0
- # computing the spectrum directly to the frequency points gives
- # wrong results
- spect = interp1d(frequency, spectral_density)
- kloss_beam = ring.f0 * np.sum(ReZ*spect(pf0))
-
- return kloss_beam
-
-
diff --git a/collective_effects/wakefield.py b/collective_effects/wakefield.py
index 4994c4a7a07e53e3660fa20e0c76ed0163d018b6..095090b7c6509cf443b943494c1efd738086d0b0 100644
--- a/collective_effects/wakefield.py
+++ b/collective_effects/wakefield.py
@@ -12,11 +12,11 @@ import re
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
-import collective_effects.tools as tools
from scipy.interpolate import interp1d
-from tracking.element import Element
+from mbtrack2.tracking.element import Element
from copy import deepcopy
from scipy.integrate import trapz
+from scipy.interpolate import interp1d
class ComplexData:
"""
@@ -327,6 +327,45 @@ class Impedance(ComplexData):
power_y += 1./2.
return power_y
+ def loss_factor(self, sigma):
+ """
+ Compute the loss factor or the kick factor for a Gaussian bunch.
+ Formulas from Eq. (2) p239 and Eq.(7) p240 of [1].
+
+ Parameters
+ ----------
+ sigma : float
+ RMS bunch length in [s]
+
+ Returns
+ -------
+ kloss: float
+ Loss factor in [V/C] or kick factor in [V/C/m] depanding on the
+ impedance type.
+
+ References
+ ----------
+ [1] : Handbook of accelerator physics and engineering, 3rd printing.
+ """
+
+ positive_index = self.data.index > 0
+ frequency = self.data.index[positive_index]
+ sd = spectral_density(frequency, sigma, m=0)
+
+ if(self.impedance_type == "long"):
+ kloss = trapz(x = frequency,
+ y = 2*self.data["real"][positive_index]*sd)
+ elif(self.impedance_type == "xdip" or self.impedance_type == "ydip"):
+ kloss = trapz(x = frequency,
+ y = 2*self.data["imag"][positive_index]*sd)
+ elif(self.impedance_type == "xquad" or self.impedance_type == "yquad"):
+ kloss = trapz(x = frequency,
+ y = 2*self.data["imag"][positive_index]*sd)
+ else:
+ raise TypeError("Impedance type not recognized.")
+
+ return kloss
+
class Wakefield:
"""
Defines a Wakefield which corresponds to a single physical element which
@@ -703,7 +742,7 @@ class ImpedanceModel(Element):
legend.append(attr)
if sigma is not None:
- spect = tools.spectral_density(zero_impedance.data.index, sigma)
+ spect = spectral_density(zero_impedance.data.index, sigma)
spect = spect/spect.max()*total_imp.max()
ax.plot(sum_imp.data.index*1e-9, spect, 'r', linewidth=2.5)
@@ -731,7 +770,7 @@ class ImpedanceModel(Element):
for j, component_name in enumerate(list_components):
try:
impedance = getattr(getattr(self, attr), component_name)
- loss_array[j,i] = tools.loss_factor(impedance, sigma)
+ loss_array[j,i] = impedance.loss_factor(sigma)
except AttributeError:
pass
@@ -800,5 +839,351 @@ class ImpedanceModel(Element):
return summary
+
+def read_CST(file, impedance_type='long', divide_by=None):
+ """
+ Read CST file format into an Impedance object.
+
+ Parameters
+ ----------
+ file : str
+ Path to the file to read.
+ impedance_type : str, optional
+ Type of the Impedance object.
+ divide_by : float, optional
+ Divide the impedance by a value. Mainly used to normalize transverse
+ impedance by displacement.
+
+ Returns
+ -------
+ result : Impedance object
+ Data from file.
+ """
+ df = pd.read_csv(file, comment="#", header = None, sep = "\t",
+ names = ["Frequency","Real","Imaginary"])
+ df["Frequency"] = df["Frequency"]*1e9
+ if divide_by is not None:
+ df["Real"] = df["Real"]/divide_by
+ df["Imaginary"] = df["Imaginary"]/divide_by
+ if impedance_type == "long":
+ df["Real"] = np.abs(df["Real"])
+ df.set_index("Frequency", inplace = True)
+ result = Impedance(variable = df.index,
+ function = df["Real"] + 1j*df["Imaginary"],
+ impedance_type=impedance_type)
+ return result
+
+def read_IW2D(file, impedance_type='long'):
+ """
+ Read IW2D file format into an Impedance object.
+
+ Parameters
+ ----------
+ file : str
+ Path to the file to read.
+ impedance_type : str, optional
+ Type of the Impedance object.
+
+ Returns
+ -------
+ result : Impedance object
+ Data from file.
+ """
+ df = pd.read_csv(file, sep = ' ', header = None,
+ names = ["Frequency","Real","Imaginary"], skiprows=1)
+ df.set_index("Frequency", inplace = True)
+ result = Impedance(variable = df.index,
+ function = df["Real"] + 1j*df["Imaginary"],
+ impedance_type=impedance_type)
+ return result
+
+def spectral_density(frequency, sigma, m = 1, 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].
+
+ Parameters
+ ----------
+ frequency : list or numpy array
+ sample points of the spectral density in [Hz]
+ sigma : float
+ RMS bunch length in [s]
+ m : int, optional
+ head-tail (or azimutal/synchrotron) mode number
+ mode: str, optional
+ type of the mode taken into account for the computation:
+ -"Hermite" modes for Gaussian bunches
+
+ Returns
+ -------
+ numpy array
+
+ References
+ ----------
+ [1] : Handbook of accelerator physics and engineering, 3rd printing.
+ """
+
+ if mode == "Hermite":
+ return (2*np.pi*frequency*sigma)**(2*m)*np.exp(
+ -1*(2*np.pi*frequency*sigma)**2)
+ else:
+ raise NotImplementedError("Not implemanted yet.")
+
+def gaussian_bunch_spectrum(frequency, sigma):
+ """
+ Compute a Gaussian bunch spectrum [1].
+
+ Parameters
+ ----------
+ frequency : array
+ sample points of the beam spectrum in [Hz].
+ sigma : float
+ RMS bunch length in [s].
+
+ Returns
+ -------
+ bunch_spectrum : array
+ Bunch spectrum sampled at points given in frequency.
+
+ References
+ ----------
+ [1] : Gamelin, A. (2018). Collective effects in a transient microbunching
+ regime and ion cloud mitigation in ThomX. p86, Eq. 4.19
+ """
+ return np.exp(-1/2*(2*np.pi*frequency)**2*sigma**2)
+
+
+def beam_spectrum(frequency, M, bunch_spacing, sigma=None,
+ bunch_spectrum=None):
+ """
+ Compute the beam spectrum assuming constant spacing between bunches [1].
+
+ Parameters
+ ----------
+ frequency : list or numpy array
+ sample points of the beam spectrum in [Hz].
+ M : int
+ Number of bunches.
+ bunch_spacing : float
+ Time between two bunches in [s].
+ sigma : float, optional
+ If bunch_spectrum is None then a Gaussian bunch with sigma RMS bunch
+ length in [s] is assumed.
+ bunch_spectrum : array, optional
+ Bunch spectrum sampled at points given in frequency.
+
+ Returns
+ -------
+ beam_spectrum : array
+
+ References
+ ----------
+ [1] Rumolo, G - Beam Instabilities - CAS - CERN Accelerator School:
+ Advanced Accelerator Physics Course - 2014, Eq. 9
+ """
+
+ if bunch_spectrum is None:
+ bunch_spectrum = Gaussian_bunch_spectrum(frequency, sigma)
+
+ beam_spectrum = (bunch_spectrum * np.exp(1j*np.pi*frequency *
+ bunch_spacing*(M-1)) *
+ np.sin(M*np.pi*frequency*bunch_spacing) /
+ np.sin(np.pi*frequency*bunch_spacing))
+
+ return beam_spectrum
+
+
+def effective_impedance(ring, imp, m, mu, sigma, xi=None, mode="Hermite"):
+ """
+ Compute the effective (longitudinal or transverse) impedance.
+ Formulas from Eq. (1) and (2) p238 of [1].
+
+ Parameters
+ ----------
+ ring : Synchrotron object
+ imp : Impedance object
+ mu : int
+ coupled bunch mode number, goes from 0 to (M-1) where M is the
+ number of bunches
+ m : int
+ head-tail (or azimutal/synchrotron) mode number
+ sigma : float
+ RMS bunch length in [s]
+ xi : float, optional
+ (non-normalized) chromaticity
+ mode: str, optional
+ type of the mode taken into account for the computation:
+ -"Hermite" modes for Gaussian bunches
+
+ Returns
+ -------
+ Zeff : float
+ effective impedance in [ohm] or in [ohm/m] depanding on the impedance
+ type.
+
+ [1] Handbook of accelerator physics and engineering, 3rd printing.
+ """
+
+ if not isinstance(imp, Impedance):
+ raise TypeError("{} should be an Impedance object.".format(imp))
+
+ fmin = imp.data.index.min()
+ fmax = imp.data.index.max()
+ if fmin > 0:
+ raise ValueError("A double sided impedance spectrum, with positive and"
+ " negative frequencies, is needed to compute the "
+ "effective impedance.")
+
+ if mode == "Hermite":
+ def h(f):
+ return spectral_density(frequency=f, sigma=sigma, m=m,
+ mode="Hermite")
+ else:
+ raise NotImplementedError("Not implemanted yet.")
+
+ M = 416
+
+ pmax = fmax/(ring.f0 * M) - 50
+ pmin = fmin/(ring.f0 * M) + 50
+
+ p = np.arange(pmin,pmax+1)
+
+ if imp.impedance_type == "long":
+ fp = ring.f0*(p*M + mu + m*ring.tune[2])
+ fp = fp[np.nonzero(fp)] # Avoid division by 0
+ num = np.sum( imp(fp) * h(fp) / (fp*2*np.pi) )
+ den = np.sum( h(fp) )
+ Zeff = num/den
+
+ elif imp.impedance_type == "xdip" or imp.impedance_type == "ydip":
+ if imp.impedance_type == "xdip":
+ tuneXY = ring.tune[0]
+ xi = ring.chro[0]
+ elif imp.impedance_type == "ydip":
+ tuneXY = ring.tune[1]
+ xi = ring.chro[1]
+ fp = ring.f0*(p*M + mu + tuneXY + m*ring.tuneS)
+ f_xi = xi/ring.eta*ring.f0
+ num = np.sum( imp(fp) * h(fp - f_xi) )
+ den = np.sum( h(fp - f_xi) )
+ Zeff = num/den
+ else:
+ raise TypeError("Effective impedance is only defined for long, xdip"
+ " and ydip impedance type.")
+
+ return Zeff
+
+
+def yokoya_elliptic(x_radius , y_radius):
+ """
+ Compute Yokoya factors for an elliptic beam pipe.
+ Function adapted from N. Mounet IW2D.
+
+ Parameters
+ ----------
+ x_radius : float
+ Horizontal semi-axis of the ellipse in [m].
+ y_radius : float
+ Vertical semi-axis of the ellipse in [m].
+
+ Returns
+ -------
+ yoklong : float
+ Yokoya factor for the longitudinal impedance.
+ yokxdip : float
+ Yokoya factor for the dipolar horizontal impedance.
+ yokydip : float
+ Yokoya factor for the dipolar vertical impedance.
+ yokxquad : float
+ Yokoya factor for the quadrupolar horizontal impedance.
+ yokyquad : float
+ Yokoya factor for the quadrupolar vertical impedance.
+ """
+ if y_radius < x_radius:
+ small_semiaxis = y_radius
+ large_semiaxis = x_radius
+ else:
+ small_semiaxis = x_radius
+ large_semiaxis = y_radius
+
+ # read Yokoya factors interpolation file
+ # BEWARE: columns are ratio, dipy, dipx, quady, quadx
+ yokoya_file = pd.read_csv("collective_effects/data/" +
+ "Yokoya_elliptic_from_Elias_USPAS.csv")
+ ratio_col = yokoya_file["x"]
+ # compute semi-axes ratio (first column of this file)
+ ratio = (large_semiaxis - small_semiaxis)/(large_semiaxis + small_semiaxis)
+
+ # interpolate Yokoya file at the correct ratio
+ yoklong = 1
+
+ if y_radius < x_radius:
+ yokydip = np.interp(ratio, ratio_col, yokoya_file["dipy"])
+ yokxdip = np.interp(ratio, ratio_col, yokoya_file["dipx"])
+ yokyquad = np.interp(ratio, ratio_col, yokoya_file["quady"])
+ yokxquad = np.interp(ratio, ratio_col, yokoya_file["quadx"])
+ else:
+ yokxdip = np.interp(ratio, ratio_col, yokoya_file["dipy"])
+ yokydip = np.interp(ratio, ratio_col, yokoya_file["dipx"])
+ yokxquad = np.interp(ratio, ratio_col, yokoya_file["quady"])
+ yokyquad = np.interp(ratio, ratio_col, yokoya_file["quadx"])
+
+ return (yoklong, yokxdip, yokydip, yokxquad, yokyquad)
+
+def beam_loss_factor(impedance, frequency, spectrum, ring):
+ """
+ Compute "beam" loss factor using the beam spectrum, uses a sum instead of
+ integral compared to loss_factor [1].
+
+ Parameters
+ ----------
+ impedance : Impedance of type "long"
+ frequency : array
+ Sample points of spectrum.
+ spectrum : array
+ Beam spectrum to consider.
+ ring : Synchrotron object
+
+ Returns
+ -------
+ kloss_beam : float
+ Beam loss factor in [V/C].
+
+ References
+ ----------
+ [1] : Handbook of accelerator physics and engineering, 3rd printing.
+ Eq (3) p239.
+ """
+ pmax = np.floor(impedance.data.index.max()/ring.f0)
+ pmin = np.floor(impedance.data.index.min()/ring.f0)
+
+ if pmin >= 0:
+ # Add negative part
+ negative_index = impedance.data.index*-1
+ negative_data = impedance.data.set_index(negative_index)
+ negative_data["imag"] = -1*negative_data["imag"]
+ try:
+ negative_data = negative_data.drop(0)
+ except KeyError:
+ pass
+
+ all_data = impedance.data.append(negative_data)
+ all_data = all_data.sort_index()
+ impedance.data = all_data
+
+ pmin = -1*pmax
+
+ p = np.arange(pmin+1,pmax)
+ pf0 = p*ring.f0
+ ReZ = np.real(impedance(pf0))
+ spectral_density = np.abs(spectrum)**2
+ # interpolation of the spectrum is needed to avoid problems liked to
+ # division by 0
+ # computing the spectrum directly to the frequency points gives
+ # wrong results
+ spect = interp1d(frequency, spectral_density)
+ kloss_beam = ring.f0 * np.sum(ReZ*spect(pf0))
+
+ return kloss_beam
\ No newline at end of file
diff --git a/machines/__init__.py b/machines/__init__.py
deleted file mode 100644
index c67145e06e038f4a48ceda33f888dc60ab733eea..0000000000000000000000000000000000000000
--- a/machines/__init__.py
+++ /dev/null
@@ -1,7 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Tue Jan 14 18:11:33 2020
-
-@author: gamelina
-"""
-
diff --git a/machines/soleil.py b/machines/soleil.py
deleted file mode 100644
index c8ac723c058f7c2d1860dc7859939cd7b0adbf55..0000000000000000000000000000000000000000
--- a/machines/soleil.py
+++ /dev/null
@@ -1,63 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-SOLEIL synchrotron parameters script.
-
-@author: Alexis Gamelin
-@date: 14/01/2020
-"""
-
-import numpy as np
-from tracking.synchrotron import Synchrotron
-from tracking.optics import Optics
-from tracking.particles import Electron, Beam
-
-def soleil(mode = 'Uniform'):
- """
-
- """
-
- h = 416
- L = 3.540969742590899e+02
- E0 = 2.75e9
- particle = Electron()
- ac = 1.47e-4
- U0 = 0.310e6
- tau = np.array([6.56e-3, 6.56e-3, 3.27e-3])
- tune = np.array([18.15687, 10.22824, 0.00502])
- emit = np.array([4.5e-9, 4.5e-9*0.01])
- sigma_0 = 8e-12
- sigma_delta = 8.6e-4
- chro = [2,3]
-
- # mean values
- beta = np.array([3, 1.3])
- alpha = np.array([0, 0])
- dispersion = np.array([0, 0, 0, 0])
- optics = Optics(local_beta=beta, local_alpha=alpha,
- local_dispersion=dispersion)
-
- ring = Synchrotron(h, optics, particle, L=L, E0=E0, ac=ac, U0=U0, tau=tau,
- emit=emit, tune=tune, sigma_delta=sigma_delta,
- sigma_0=sigma_0, chro=chro)
-
- return ring
-
-def v0313():
- """
-
- """
-
- h = 416
- particle = Electron()
- tau = np.array([7.3e-3, 13.1e-3, 11.7e-3])
- emit = np.array([53.47e-12, 53.47e-12])
- sigma_0 = 9.2e-12
- sigma_delta = 9e-4
-
- lattice_file = "SOLEIL_U_74BA_HOA_SYM02_V0313_cleaned.mat"
- optics = Optics(lattice_file=lattice_file)
-
- ring = Synchrotron(h, optics, particle, tau=tau, emit=emit,
- sigma_0=sigma_0, sigma_delta=sigma_delta)
-
- return ring
\ No newline at end of file
diff --git a/plotting.py b/plotting.py
deleted file mode 100644
index f03c40eefa605069541e65ebb0ad359a86036018..0000000000000000000000000000000000000000
--- a/plotting.py
+++ /dev/null
@@ -1,74 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Running script to test getplot method in Bunch class, particles module
-
-@author: Watanyu Foosang
-@date: 27/03/2020
-"""
-
-import numpy as np
-import matplotlib.pyplot as plt
-import seaborn as sns
-import pandas as pd
-from tracking.particles import Bunch, Beam
-from machines.soleil import soleil
-from tracking.rf import *
-from scipy.constants import e
-
-ring = soleil() # Define the properties of the storage ring
-mybunch = Bunch(ring,mp_number=1000,current=1e-3) # Define bunch parameters
-
-mybunch.init_gaussian() # Initialize 6D bunch parameters
-
-
-v_rf = 1.1e6 # RF peak voltage [v]
-sync_phase = np.arccos(ring.U0/v_rf) # Determine the synchronous phase
-
-rf = RFCavity(ring, m=1, Vc=v_rf, theta=sync_phase)
-
-rev = []
-mean_tau = []
-mean_delta = []
-def track_long(bunch,ring,rf,turn):
- """
- Track longitudinal beam parameters for one turn
- """
-
- rf.track(bunch) # update "delta" for one turn
- bunch.particles["tau"] = bunch.particles["tau"] - bunch.particles["delta"]\
- *ring.T0*ring.ac
- turn += 1
- rev.append(turn)
- mean_tau.append(bunch.particles["tau"].mean())
- mean_delta.append(bunch.particles["delta"].mean())
- return turn
-
-def get_meanplot(ini,fin,step,val):
- """
- Plot the evolution of a mean value over turns.
- (For now, this function only works subject to plotting module.)
-
- Parameters
- ----------
- ini: int, starting revolution
- fin: int, final revolution
- step: int, step size on the revolution axis
- val: str, "tau" for maen values of tau or "del" for mean values of delta
-
- """
-
- x = rev[ini-1:fin:step]
-
- if val == "tau":
- val = mean_tau[ini-1:fin:step]
- else: val = mean_delta[ini-1:fin:step]
-
- plt.plot(x,val)
-
-
-
-
-
-
-
-
diff --git a/tracking/__init__.py b/tracking/__init__.py
index c67145e06e038f4a48ceda33f888dc60ab733eea..45db00132ecab8571c819849c041e992ca8826fe 100644
--- a/tracking/__init__.py
+++ b/tracking/__init__.py
@@ -5,3 +5,14 @@ Created on Tue Jan 14 18:11:33 2020
@author: gamelina
"""
+from mbtrack2.tracking.particles import (Electron, Proton, Bunch, Beam,
+ Particle)
+from mbtrack2.tracking.synchrotron import Synchrotron
+from mbtrack2.tracking.rf import RFCavity
+from mbtrack2.tracking.parallel import Mpi
+from mbtrack2.tracking.optics import Optics, PhyisicalModel
+from mbtrack2.tracking.element import (Element, LongitudinalMap,
+ TransverseMap, SynchrotronRadiation)
+from mbtrack2.tracking.aperture import (CircularAperture, ElipticalAperture,
+ RectangularAperture)
+
diff --git a/tracking/aperture.py b/tracking/aperture.py
index 4cba4b3de7d395777e340a22104b3cd8d9308498..177184563831d4ef5a93bce825127240c6d7a6c0 100644
--- a/tracking/aperture.py
+++ b/tracking/aperture.py
@@ -7,7 +7,7 @@ This module defines aperture elements for tracking.
"""
import numpy as np
-from tracking.element import Element
+from mbtrack2.tracking.element import Element
class CircularAperture(Element):
"""
diff --git a/tracking/element.py b/tracking/element.py
index bee031ce02323679e53832f4745196ef33dc1666..8692856c1c54e9f17f7a40fbed11bc693b0273f9 100644
--- a/tracking/element.py
+++ b/tracking/element.py
@@ -11,7 +11,7 @@ included in the tracking.
import numpy as np
from abc import ABCMeta, abstractmethod
from functools import wraps
-from tracking.particles import Beam
+from mbtrack2.tracking.particles import Beam
class Element(metaclass=ABCMeta):
"""
diff --git a/tracking/monitors/__init__.py b/tracking/monitors/__init__.py
index c67145e06e038f4a48ceda33f888dc60ab733eea..028472c4c7fe270d4434391af894d9fcc56a35eb 100644
--- a/tracking/monitors/__init__.py
+++ b/tracking/monitors/__init__.py
@@ -5,3 +5,9 @@ Created on Tue Jan 14 18:11:33 2020
@author: gamelina
"""
+from mbtrack2.tracking.monitors.monitors import (Monitor, BunchMonitor,
+ PhaseSpaceMonitor,
+ BeamMonitor,
+ ProfileMonitor)
+from mbtrack2.tracking.monitors.plotting import (plot_bunchdata,
+ plot_phasespacedata)
\ No newline at end of file
diff --git a/tracking/monitors/monitors.py b/tracking/monitors/monitors.py
index e6234651dc30eed6693b8200f86ded7b96d7fbe9..c7c6ac86173bda7087ac8e529879671b9ff32544 100644
--- a/tracking/monitors/monitors.py
+++ b/tracking/monitors/monitors.py
@@ -10,8 +10,8 @@ during tracking.
import numpy as np
import h5py as hp
-from tracking.element import Element
-from tracking.particles import Bunch, Beam
+from mbtrack2.tracking.element import Element
+from mbtrack2.tracking.particles import Bunch, Beam
from abc import ABCMeta
from mpi4py import MPI
diff --git a/tracking/particles.py b/tracking/particles.py
index cb4c0468306f63b81fb9d45fc56f0127713b6f83..d3e7f171ab827c2aeec740a19c7cc17c940ee47d 100644
--- a/tracking/particles.py
+++ b/tracking/particles.py
@@ -10,7 +10,7 @@ Module where particles, bunches and beams are described as objects.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
-from tracking.parallel import Mpi
+from mbtrack2.tracking.parallel import Mpi
import seaborn as sns
from scipy.constants import c, m_e, m_p, e
diff --git a/tracking/rf.py b/tracking/rf.py
index 24928d33841f52697dccbf31a5ccf970ea8d889a..451d220c88b6bda6ea0aaafa7f00775a29699909 100644
--- a/tracking/rf.py
+++ b/tracking/rf.py
@@ -7,7 +7,7 @@ This module handles radio-frequency (RF) cavitiy elements.
"""
import numpy as np
-from tracking.element import Element
+from mbtrack2.tracking.element import Element
class RFCavity(Element):
"""