# -*- coding: utf-8 -*-
"""
Module for intrabeam scattering computations.
"""
import warnings
import numpy as np
import scipy.integrate as quad
from numpy.typing import NDArray
from scipy.constants import c, elementary_charge, epsilon_0
from scipy.special import hyp2f1
from mbtrack2.tracking.element import Element
from mbtrack2.tracking.particles import Beam, Bunch
from mbtrack2.utilities.synchrotron import Synchrotron
[docs]
class IntrabeamScattering(Element):
"""
Intrabeam scattering class computes the IBS growth rate analytically each
turn and apply corresponding kicks to each particle according to the
following models:
Piwinski Modified (PM):
Modification of the original Piwinski model [1] in order to take into account the
dispersion and the variation of the optics functions around the ring [2].
Bane (Bane):
Approximation for high-energy beams which allows computing a single
scattering integral instead of two [3].
Completely integrated modified Piwinski model (CIMP):
Approximating elliptic integral as a special Legendre function [4].
Parameters:
-----------
ring: Synchrotron object
Ring to consider.
bunch: Bunch or Beam object
Bunch or Beam object which will be tracked.
model: {'CIMP','Bane','PM','PS'}
Implemented computational models for computing the growth rate T_(x,y,p):
CIMP: the CIMP model [4].
Bane: the Bane model [3].
PM: the modified model of Piwinski [2].
n_points: int, optional
Number of partitions (sampling) of the optics functions (depends on the
lattice complexity), can use a lower value (50) if the beta function is
symmetric, thus gaining computation speed.
The Default is 1000
n_bin: int, optional
Number of bins for the bunch.binning profile function. Will determine
the number of slices for the profile function.
If the value of n_bin is greater than 1 the kick method will compute
the kick[5], applying momentum change according to the positions
of each macroparticle with respect to the density of macroparticles in
that position.
If n_bin is set to 1, the kick will be computed assuming Rho(z) as a
uniform distribution.
The Default is 100.
Methods:
--------
initialize(bunch):
Initializes the dynamic parameters at each turn, modifies the class
variables according to the selected model.
get_scattering_integrals():
Computes the scattering integrals according to the selected model.
get_ibs_growthrate():
Computes the growth rate from the scattering integrals.
kick(bunch, T_x, T_y, T_p):
Tracking method of IntrabeamScattering takes T_(x,y,p) and apply
momentum change to the coordinates of macroparticles inside the bunch.
References:
-----------
[1] A. Piwinski, Intra-Beam-Scattering, (1974).
http://dx.doi.org/10.5170/CERN-1992-001.226
[2] K. L. F. Bane, A Simplified Model of Intrabeam Scattering, in 8th
European Particle Accelerator Conference (Paris, France, 2002), p. 1443.
https://doi.org/10.48550/arXiv.physics/0206002
[3] K. L. F. Bane, H. Hayano, K. Kubo, T. Naito, T. Okugi, and J. Urakawa,
Intrabeam Scattering Analysis of Measurements at KEK’s Accelerator Test
Facility Damping Ring, Phys. Rev. ST Accel. Beams 5, 084403 (2002).
https://doi.org/10.1103/PhysRevSTAB.5.084403
[4] K. Kubo, S. K. Mtingwa, and A. Wolski, Intrabeam Scattering Formulas
for High Energy Beams, Phys. Rev. ST Accel. Beams 8, 081001 (2005).
https://doi.org/10.1103/PhysRevSTAB.8.081001
[5] R. Bruce, J. M. Jowett, M. Blaskiewicz, and W. Fischer, Time
Evolution of the Luminosity of Colliding Heavy-Ion Beams in BNL
Relativistic Heavy Ion Collider and CERN Large Hadron Collider,
Phys. Rev. ST Accel. Beams 13, 091001 (2010).
https://doi.org/10.1103/PhysRevSTAB.13.091001
"""
[docs]
def __init__(
self,
ring: Synchrotron,
model: str,
n_points: int = 1000,
n_bin: int = 100,
):
self.ring = ring
self.n_points = int(n_points)
self.n_bin = int(n_bin)
self.model = model
self.s = np.linspace(0, ring.L, n_points)
self.beta_x, self.beta_y = self.ring.optics.beta(self.s)
self.dispX, self.disppX, self.dispY, self.disppY = self.ring.optics.dispersion(
self.s)
self.alphaX, self.alphaY = self.ring.optics.alpha(self.s)
self.H_x = (1 / self.beta_x) * (self.dispX**2 +
((self.beta_x * self.disppX) +
(self.alphaX * self.dispX))**2)
self.H_y = (1 / self.beta_y) * (self.dispY**2 +
((self.beta_y * self.disppY) +
(self.alphaY * self.dispY))**2)
self.local_dispersion_x = self.ring.optics.dispersion(0)[1]
self.local_dispersion_y = self.ring.optics.dispersion(0)[3]
if self.ring.optics.use_local_values:
warnings.warn(
"Lattice file not loaded, intiating optics with approximated values",
UserWarning)
self.n_points = 1
if self.model not in ["Bane", "CIMP", "PM"]:
raise ValueError(
"Incorrect model name. Allowed model names are: PM, Bane, CIMP."
)
self.r_0 = (ring.particle.charge**
2) / (4 * np.pi * epsilon_0 * c**2 * ring.particle.mass)
def __repr__(self):
return f"IntrabeamScattering(model={self.model}, " \
"n_points={self.n_points}, n_bin={self.n_bin})"
[docs]
def initialize(self, bunch: Bunch):
"""
Initializes the dynamic parameters at each turn, modifies the class
variables according to the selected model.
Parameters
----------
bunch: Bunch or Beam object
Bunch or Beam object which will be tracked.
"""
self.N = bunch.current * self.ring.T0 / elementary_charge
self.sigma_s = np.std(bunch['tau'])
self.sigma_p = np.std(bunch['delta'])
self.emit_x, self.emit_y, _ = bunch.emit
self.d = 4 * np.min((np.std(bunch['y']), np.std(bunch['x'])))
self.sigma_px = np.std(bunch['xp'] -
bunch['delta'] * self.local_dispersion_x)
self.sigma_py = np.std(bunch['yp'] -
bunch['delta'] * self.local_dispersion_y)
self.H = (1 / self.sigma_p**2) + (self.H_x / self.emit_x) + \
(self.H_y / self.emit_y)
sigma_H = np.sqrt(1 / self.H)
self.sigma_H = sigma_H
self.a = (self.sigma_H / self.ring.gamma) * np.sqrt(
self.beta_x / self.emit_x)
self.b = (self.sigma_H / self.ring.gamma) * np.sqrt(
self.beta_y / self.emit_y)
self.q = self.sigma_H * self.ring.beta * np.sqrt(2 * self.d / self.r_0)
self.C_log = np.log(self.q**2 / self.a**2)
self.C_a = self.a / self.b
self.A = (self.r_0**2 * self.N) / (
64 * np.pi**2 * self.ring.beta**3 * self.ring.gamma**4 *
self.emit_x * self.emit_y * self.sigma_s * self.sigma_p)
[docs]
def get_scattering_integrals(self, bunch):
"""
Computes the scattering integrals according to the selected model.
Parameters
----------
bunch: Bunch or Beam object
Bunch or Beam object which will be tracked.
Returns
-------
If self.model == "PM":
vabq, v1aq, v1bq: arrays
Array of scattering integral values at each point around the
ring.
If self.model == "Bane":
gval: array
Array of scattering integral values at each point around the
ring.
If self.model == "CIMP":
g_ab, g_ba: arrays
Array of analytical function g(a/b) and g(b/a) respectively,
that simulate the integrals at each point around the ring.
"""
self.initialize(bunch)
v_abq = np.zeros(self.n_points, dtype=np.float64)
v_1aq = np.zeros(self.n_points, dtype=np.float64)
v_1bq = np.zeros(self.n_points, dtype=np.float64)
if self.model == "PM":
def scattering(u, x, y, z):
"""
Eq. (17) in:
L. R. Evans and B. W. Zotter, Intrabeam Scattering in the SPS.
https://cds.cern.ch/record/126036
"""
P = np.sqrt(x**2 + (1 - x**2) * u**2)
Q = np.sqrt(y**2 + (1 - y**2) * u**2)
f_abq = 8 * np.pi * (1 - 3 * u**2) / (P * Q) * \
(2 * np.log(z / 2 * (1/P + 1/Q)) - 0.5777777777)
return f_abq
for i in range(self.n_points):
el_1aq, _ = quad.quad(scattering,
0,
1,
args=(1 / self.b[i],
self.a[i] / self.b[i],
self.q[i] / self.b[i]))
el_1bq, _ = quad.quad(scattering,
0,
1,
args=(1 / self.a[i],
self.b[i] / self.a[i],
self.q[i] / self.a[i]))
el_abq = -el_1aq / self.b[i]**2 - el_1bq / self.a[i]**2
v_abq[i] = el_abq
v_1aq[i] = el_1aq
v_1bq[i] = el_1bq
elif self.model == "Bane":
def _g_bane(u, j, C_a):
"""
Eq. (12) in [2].
Parameters
----------
u : float
integration variable.
j : int
index.
C_a : float
result of a/b
Returns
-------
g_val : array
Scattering integral value at a given point.
"""
g_val = 2 * np.sqrt(C_a[j])/ np.pi * \
(1 / (np.sqrt(1 + u**2) * np.sqrt(C_a[j]**2 + u**2)))
return g_val
for j in range(self.n_points):
reslt, _ = quad.quad(_g_bane, 0, np.inf, args=(j, self.C_a))
v_abq[j] = 4 * np.pi**2 * reslt * self.C_log[j] / np.sqrt(
self.a[j] * self.b[j])
elif self.model == "CIMP":
def _Puv(u, v, x):
"""
https://dlmf.nist.gov/14.3
"""
if x < 1:
val = ((1 + x) / (1 - x))**(u/2) * \
hyp2f1(v+1, -v, 1-u, .5 - .5 * x)
else:
val = ((1 + x) / (x - 1))**(u/2) * \
hyp2f1(v+1, -v, 1-u, .5 - .5 * x)
return val
@np.vectorize
def _g_cimp(u):
"""
Eq. (34) in [4].
"""
x_arg = (1 + u**2) / (2*u)
if u >= 1:
g_val = np.sqrt(np.pi / u) * (_Puv(0, -.5, x_arg) +
3 / 2 * _Puv(-1, -.5, x_arg))
else:
g_val = np.sqrt(np.pi / u) * (_Puv(0, -.5, x_arg) -
3 / 2 * _Puv(-1, -.5, x_arg))
return g_val
g_ab = _g_cimp(self.a / self.b)
g_ba = _g_cimp(self.b / self.a)
v_abq = 4 * np.pi**(3 /
2) * (g_ab / self.b * np.log(self.q / self.b) +
g_ba / self.a * np.log(self.q / self.a))
v_1aq = self.b * g_ab * np.log(self.q / self.b)
v_1bq = self.a * g_ba * np.log(self.q / self.a)
else:
raise ValueError("Wrong model, must be 'PM', 'Bane' or 'CIMP'.")
return v_abq, v_1aq, v_1bq
[docs]
def get_ibs_growthrate(self, bunch):
"""
Computes the growth rate from the scattering integrals:
Piwinski: Eq. (16-18) in [2].
Bane: Eq (14) and (20) in [3].
CIMP: Eq (35-37) in [4].
The growth rate can be defined as:
1/T_ibs = 1/sigma * d(sigma)/dt
Parameters
----------
vabq, v1aq, v1bq: arrays
Takes integral values from scatter() method.
Returns
-------
r_x: Array
Array of IBS growth rates on the horizontal plane over the entire
ring in [1/s].
r_y: Array
Array of IBS growth rates on the vertical plane over the entire
ring in [1/s].
r_p: Array
Array of IBS growth rates on the longitudinal plane over the entire
ring in [1/s].
"""
vabq, v1aq, v1bq = self.get_scattering_integrals(bunch)
r_p = self.A * vabq * self.sigma_H**2 / self.sigma_p**2
r_x = self.A * (v1bq + vabq * self.H_x * self.sigma_H**2 / self.emit_x)
r_y = self.A * (v1aq + vabq * self.H_y * self.sigma_H**2 / self.emit_y)
return r_x, r_y, r_p
[docs]
def kick(self, bunch: Bunch, r_x: float, r_y: float, r_p: float):
"""
Apply kick to the bunch coordinates by converting growth rate to
momentum change [1].
Parameters
----------
bunch : Object
Bunch or Beam object.
r_p : Array
Growth rate of 'delta' in [1/s].
r_x : Array
Growth rate of 'xp' in [1/s].
r_y : Array
Growth rate of 'yp' in [1/s].
Reference:
----------
[1] R. Bruce, J. M. Jowett, M. Blaskiewicz, and W. Fischer, Time
Evolution of the Luminosity of Colliding Heavy-Ion Beams in BNL
Relativistic Heavy Ion Collider and CERN Large Hadron Collider,
Phys. Rev. ST Accel. Beams 13, 091001 (2010).
https://doi.org/10.1103/PhysRevSTAB.13.091001
"""
# Growth rate operations are moved here to get real values incase of benchmarking with other codes
r_p[r_p <= 0] = 0
r_x[r_x <= 0] = 0
r_y[r_y <= 0] = 0
r_p = np.average(r_p)
r_x = np.average(r_x)
r_y = np.average(r_y)
if self.n_bin > 1:
_, sorted_index, profile, _ = bunch.binning(n_bin=self.n_bin)
normalized_profile = profile / max(profile)
Rho = normalized_profile[sorted_index]
else:
Rho = 1.0
N_mp = len(bunch)
bunch['delta'] += self.sigma_p * np.sqrt(
np.sqrt(2) * r_p * self.ring.T0 *
Rho) * np.random.normal(size=N_mp)
bunch['xp'] += self.sigma_px * np.sqrt(
np.sqrt(2) * r_x * self.ring.T0 *
Rho) * np.random.normal(size=N_mp)
bunch['yp'] += self.sigma_py * np.sqrt(
np.sqrt(2) * r_y * self.ring.T0 *
Rho) * np.random.normal(size=N_mp)
[docs]
@Element.parallel
@Element.track_bunch_if_non_empty
def track(self, bunch: Bunch | Beam):
"""
Tracking method of IntrabeamScattering takes T_(x,y,p) and apply
momentum change to the coordinates of macroparticles inside the bunch.
parameters:
-----------
bunch: Object
Bunch or Beam object.
"""
rx, ry, rp = self.get_ibs_growthrate(bunch)
self.kick(bunch, rx, ry, rp)