From 081815713056d7a017d14f92e2b4a2b32ab7fe83 Mon Sep 17 00:00:00 2001
From: Alexis Gamelin <alexis.gamelin@synchrotron-soleil.fr>
Date: Sat, 6 Apr 2024 22:44:04 +0200
Subject: [PATCH] BeamLoadingEquilibrium on a grid
---
mbtrack2/utilities/beamloading.py | 116 ++++++++++++++----------------
1 file changed, 53 insertions(+), 63 deletions(-)
diff --git a/mbtrack2/utilities/beamloading.py b/mbtrack2/utilities/beamloading.py
index 19e51e9..cfd3548 100644
--- a/mbtrack2/utilities/beamloading.py
+++ b/mbtrack2/utilities/beamloading.py
@@ -6,18 +6,14 @@ Module where the BeamLoadingEquilibrium class is defined.
import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import c
-from scipy.integrate import quad
+from scipy.integrate import trapz
from scipy.optimize import root
-
+from mbtrack2.utilities.spectrum import gaussian_bunch
class BeamLoadingEquilibrium():
"""Class used to compute beam equilibrium profile for a given storage ring
and a list of RF cavities of any harmonic. The class assumes an uniform
- filling of the storage ring. Based on an extension of [1].
-
- [1] Venturini, M. (2018). Passive higher-harmonic rf cavities with general
- settings and multibunch instabilities in electron storage rings.
- Physical Review Accelerators and Beams, 21(11), 114404.
+ filling of the storage ring. Based on [1] which is an extension of [2].
Parameters
----------
@@ -30,8 +26,20 @@ class BeamLoadingEquilibrium():
PHI : list of form factor phase
B1 : lower intergration boundary
B2 : upper intergration boundary
+ N : int or float, optinal
+ Number of points used for longitudinal grid.
+ Default is 1000.
+
+ References
+ ----------
+ [1] Gamelin, A., & Yamamoto, N. (2021). Equilibrium bunch density
+ distribution with multiple active and passive RF cavities.
+ IPAC'21 (MOPAB069).
+ [2] Venturini, M. (2018). Passive higher-harmonic rf cavities with general
+ settings and multibunch instabilities in electron storage rings.
+ Physical Review Accelerators and Beams, 21(11), 114404.
+
"""
-
def __init__(self,
ring,
cavity_list,
@@ -40,7 +48,8 @@ class BeamLoadingEquilibrium():
F=None,
PHI=None,
B1=-0.2,
- B2=0.2):
+ B2=0.2,
+ N=1000):
self.ring = ring
self.cavity_list = cavity_list
self.I0 = I0
@@ -57,7 +66,10 @@ class BeamLoadingEquilibrium():
self.B1 = B1
self.B2 = B2
self.mpi = False
- self.__version__ = "1.0"
+ self.N = N
+ self.z0 = np.linspace(self.B1, self.B2, self.N)
+ self.tau0 = self.z0/c
+ self.rho0 = gaussian_bunch(self.tau0, self.ring.sigma_0)
# Define constants for scaled potential u(z)
self.u0 = self.ring.U0 / (self.ring.ac * self.ring.sigma_delta**2 *
@@ -65,6 +77,11 @@ class BeamLoadingEquilibrium():
self.ug = np.zeros((self.n_cavity, ))
self.ub = np.zeros((self.n_cavity, ))
self.update_potentials()
+
+ def update_rho(self):
+ uexp_vals = self.uexp(self.z0)
+ A = trapz(uexp_vals, x=self.z0)
+ self.rho0 = uexp_vals / A
def update_potentials(self):
"""Update potentials with cavity and ring data."""
@@ -90,9 +107,8 @@ class BeamLoadingEquilibrium():
def center_of_mass(self):
"""Return center of mass position in [s]"""
- z0 = np.linspace(self.B1, self.B2, 1000)
- rho = self.rho(z0)
- CM = np.average(z0, weights=rho)
+ self.update_rho()
+ CM = np.average(self.z0, weights=self.rho0)
return CM / c
def u(self, z):
@@ -122,12 +138,12 @@ class BeamLoadingEquilibrium():
def uexp(self, z):
return np.exp(-1 * self.u(z))
-
+
def integrate_func(self, f, g):
"""Return Integral[f*g]/Integral[f] between B1 and B2"""
- A = quad(lambda x: f(x) * g(x), self.B1, self.B2)
- B = quad(f, self.B1, self.B2)
- return A[0] / B[0]
+ A = trapz(f*g, x=self.z0)
+ B = trapz(f, x=self.z0)
+ return A / B
def to_solve(self, x, CM=True):
"""System of non-linear equation to solve to find the form factors F
@@ -158,12 +174,12 @@ class BeamLoadingEquilibrium():
cavity = self.cavity_list[i]
res[2 *
i] = self.F[i] * np.cos(self.PHI[i]) - self.integrate_func(
- lambda y: self.uexp(y),
- lambda y: np.cos(self.ring.k1 * cavity.m * y))
+ self.uexp(self.z0),
+ np.cos(self.ring.k1 * cavity.m * self.z0))
res[2*i +
1] = self.F[i] * np.sin(self.PHI[i]) - self.integrate_func(
- lambda y: self.uexp(y),
- lambda y: np.sin(self.ring.k1 * cavity.m * y))
+ self.uexp(self.z0),
+ np.sin(self.ring.k1 * cavity.m * self.z0))
# Factor 1e-8 or 1e12 for better convergence
if CM is True:
res[self.n_cavity * 2] = self.center_of_mass() * 1e12
@@ -175,27 +191,17 @@ class BeamLoadingEquilibrium():
cavity = self.cavity_list[i]
res[2 *
i] = self.F[i] * np.cos(self.PHI[i]) - self.integrate_func(
- lambda y: self.uexp(y),
- lambda y: np.cos(self.ring.k1 * cavity.m * y))
+ self.uexp(self.z0),
+ np.cos(self.ring.k1 * cavity.m * self.z0))
res[2*i +
1] = self.F[i] * np.sin(self.PHI[i]) - self.integrate_func(
- lambda y: self.uexp(y),
- lambda y: np.sin(self.ring.k1 * cavity.m * y))
+ self.uexp(self.z0),
+ np.sin(self.ring.k1 * cavity.m * self.z0))
return res
- def rho(self, z):
- """Return bunch equilibrium profile at postion z"""
- A = quad(lambda y: self.uexp(y), self.B1, self.B2)
- return self.uexp(z) / A[0]
-
- def plot_rho(self, z1=None, z2=None):
- """Plot the bunch equilibrium profile between z1 and z2"""
- if z1 is None:
- z1 = self.B1
- if z2 is None:
- z2 = self.B2
- z0 = np.linspace(z1, z2, 1000)
- plt.plot(z0 / c * 1e12, self.rho(z0))
+ def plot_rho(self):
+ """Plot the bunch equilibrium profile"""
+ plt.plot(self.tau0 * 1e12, self.rho0)
plt.xlabel(r"$\tau$ [ps]")
plt.title("Equilibrium bunch profile")
@@ -231,38 +237,22 @@ class BeamLoadingEquilibrium():
Vtot += cavity.deltaVRF(z, self.I0, self.F[i], self.PHI[i])
return Vtot
- def plot_dV(self, z1=None, z2=None):
- """Plot the derivative of RF system total voltage between z1 and z2"""
- if z1 is None:
- z1 = self.B1
- if z2 is None:
- z2 = self.B2
- z0 = np.linspace(z1, z2, 1000)
- plt.plot(z0, self.dV(z0))
+ def plot_dV(self):
+ """Plot the derivative of RF system total voltage"""
+ plt.plot(self.z0, self.dV(self.z0))
plt.xlabel("z [m]")
plt.ylabel("Total RF voltage (V)")
- def plot_voltage(self, z1=None, z2=None):
+ def plot_voltage(self):
"""Plot the RF system total voltage between z1 and z2"""
- if z1 is None:
- z1 = self.B1
- if z2 is None:
- z2 = self.B2
- z0 = np.linspace(z1, z2, 1000)
- plt.plot(z0, self.voltage(z0))
+ plt.plot(self.z0, self.voltage(self.z0))
plt.xlabel("z [m]")
plt.ylabel("Total RF voltage (V)")
- def std_rho(self, z1=None, z2=None):
+ def std_rho(self):
"""Return the rms bunch equilibrium size in [m]"""
- if z1 is None:
- z1 = self.B1
- if z2 is None:
- z2 = self.B2
- z0 = np.linspace(z1, z2, 1000)
- values = self.rho(z0)
- average = np.average(z0, weights=values)
- variance = np.average((z0 - average)**2, weights=values)
+ average = np.average(self.z0, weights=self.rho0)
+ variance = np.average((self.z0 - average)**2, weights=self.rho0)
return np.sqrt(variance)
def beam_equilibrium(self,
@@ -333,7 +323,7 @@ class BeamLoadingEquilibrium():
print("The algorithm has converged: " + str(sol.success))
if plot:
- self.plot_rho(self.B1 / 4, self.B2 / 4)
+ self.plot_rho()
return sol
--
GitLab