From 947b01b443e731b861cf33681643f5766b9325d4 Mon Sep 17 00:00:00 2001 From: Ritu Bharadwaj Date: Thu, 10 Jul 2025 13:55:08 +0200 Subject: [PATCH 1/3] These are the files to model the IntraBeam Scattering effect in pyat. --- pyat/at/collective/__init__.py | 2 + pyat/at/collective/ibs.py | 293 ++++++++++++++++++++++++++++ pyat/at/collective/ibs_element.py | 165 ++++++++++++++++ pyat/at/integrators/pyIBSRadPass.py | 241 +++++++++++++++++++++++ pyat/test/test_ibs.py | 63 ++++++ 5 files changed, 764 insertions(+) create mode 100644 pyat/at/collective/ibs.py create mode 100644 pyat/at/collective/ibs_element.py create mode 100644 pyat/at/integrators/pyIBSRadPass.py create mode 100644 pyat/test/test_ibs.py diff --git a/pyat/at/collective/__init__.py b/pyat/at/collective/__init__.py index 264492fa5f..01d680c0ec 100644 --- a/pyat/at/collective/__init__.py +++ b/pyat/at/collective/__init__.py @@ -1,7 +1,9 @@ """ Collective effects """ + from .wake_elements import * from .wake_functions import * from .wake_object import * from .beam_loading import * +from .ibs_element import * diff --git a/pyat/at/collective/ibs.py b/pyat/at/collective/ibs.py new file mode 100644 index 0000000000..8d076a2c81 --- /dev/null +++ b/pyat/at/collective/ibs.py @@ -0,0 +1,293 @@ +""" +IBS module for Acclerator toolbox(AT) inspired from mbtrack2 +from Alexin Gamelin , Vadim Gubaidulin (SOLEIL) +Github : https://gitlab.synchrotron-soleil.fr/PA/collective-effects/mbtrack2 + +""" + +import numpy as np +import at +import sys +import scipy.integrate as quad +from scipy.constants import c, elementary_charge, epsilon_0, m_e +from scipy.special import hyp2f1 +from scipy.interpolate import interp1d +from at.constants import qe +from at.constants import e_mass +from at.constants import clight + + +def initialize(bunch, ring): + """ + calculates bunch parameters at each turn. + Parameters + --------------- + bunch object that is being tracked. + """ + + d = 4 * np.std(bunch[2, :]) + sigma_s = np.std(bunch[5, :] / clight) # sigma_s in m + sigma_p = np.std(bunch[4, :]) + sigma_px = np.std(bunch[1, :]) + sigma_py = np.std(bunch[3, :]) + + return (sigma_py, sigma_px, sigma_p, sigma_s, d) + + +def paramerters( + bunch, + ring, + model, + sigma_p, + beta_x, + beta_y, + dispX_sq_over_beta_x, + dispY_sq_over_beta_y, + H_x, + H_y, + gamma, + beta, + d, + r_0, +): + """ + Calculates emittances and derived parameters (a, b, q, sigma_H) based on the selected Intrabeam Scattering (IBS) model. + These computations are performed every 500 turns by default, or at user-defined intervals specified by the update_turns variable, + depending on the lattice configuration. + + """ + + sigma = at.sigma_matrix(beam=bunch) + eigs, _ = np.linalg.eig(sigma @ at.jmat(3)) + emits = np.sort( + np.abs(eigs)[::2] + ) # order after sorting emits = [emittance_vertical , emittance_horizontal , emittance_longitudinal] + emitx = emits[1] + emity = emits[0] + + if model == "PS": + h = ( + (1 / sigma_p**2) + + (dispX_sq_over_beta_x / emitx) + + (dispY_sq_over_beta_y / emity) + ) + + sigma_h = np.sqrt(1 / h) + sigma_H = sigma_h + + elif model in ["CIMP", "PM", "Bane"]: + H = (1 / sigma_p**2) + (H_x / emitx) + (H_y / emity) + sigma_H = np.sqrt(1 / H) + + a = (sigma_H / gamma) * np.sqrt(beta_x / emitx) + b = (sigma_H / gamma) * np.sqrt(beta_y / emity) + q = sigma_H * beta * np.sqrt(2 * d / r_0) + + return (a, b, q, sigma_H, emitx, emity) + + +def scatter(model, n_points, b, a, q, C_a=None): + """ + Computes IBS scattering integrals using model-specific formulas. + This calculation is performed every 500 turns by default or at user-defined intervals specified by the update_turns variable, + depending on the lattice configuration. + + """ + + if model in ["PS", "PM"]: + vabq = np.zeros(n_points, dtype=np.float64) + v1aq = np.zeros(n_points, dtype=np.float64) + v1bq = np.zeros(n_points, dtype=np.float64) + + 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 + """ + P2 = x**2 + ((1 - x**2) * u**2) + Q2 = y**2 + ((1 - y**2) * u**2) + P = np.sqrt(P2) + Q = np.sqrt(Q2) + 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(n_points): + el_1aq, err = quad.quad( + scattering, 0, 1, args=(1 / b[i], a[i] / b[i], q[i] / b[i]) + ) + el_1bq, err = quad.quad( + scattering, 0, 1, args=(1 / a[i], b[i] / a[i], q[i] / a[i]) + ) + el_abq = -(el_1aq * (1 / b[i] ** 2)) - (el_1bq * (1 / a[i] ** 2)) + vabq[i] = el_abq + v1aq[i] = el_1aq + v1bq[i] = el_1bq + return vabq, v1aq, v1bq + + elif model == "Bane": + gval = np.zeros(n_points, dtype=np.float64) + + def g_func(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(n_points): + reslt, err = quad.quad(g_func, 0, np.inf, args=(j, C_a)) + gval[j] = reslt + return gval + + elif 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, (0.5 - (0.5 * x)) + ) + else: + val = ((1 + x) / (x - 1)) ** (u / 2) * hyp2f1( + v + 1, -v, 1 - u, (0.5 - (0.5 * x)) + ) + return val + + def g_func(u): + """ + Eq. (34) in [4]. + """ + x_arg = (1 + u**2) / (2 * u) + if u >= 1: + g_val = np.sqrt(np.pi / u) * ( + (Puv(0, -0.5, x_arg)) + ((3 / 2) * (Puv(-1, -0.5, x_arg))) + ) + else: + g_val = np.sqrt(np.pi / u) * ( + (Puv(0, -0.5, x_arg)) - ((3 / 2) * (Puv(-1, -0.5, x_arg))) + ) + return g_val + + g_ab = np.zeros(n_points) + g_ba = np.zeros(n_points) + for i in range(n_points): + val_ab = g_func(a[i] / b[i]) + val_ba = g_func(b[i] / a[i]) + g_ab[i] = val_ab + g_ba[i] = val_ba + return g_ab, g_ba + + +def get_scatter_T( + vabq=None, + v1aq=None, + v1bq=None, + g_ab=None, + g_ba=None, + gval=None, + model=None, + A=None, + sigma_H=None, + dispX=None, + dispY=None, + beta_x=None, + beta_y=None, + r_0=None, + N=None, + C_log=None, + gamma=None, + beta=None, + emitx=None, + emity=None, + sigma_s=None, + sigma_p=None, + H_x=None, + H_y=None, + a=None, + b=None, + q=None, + dispX_sq_over_beta_x=None, + dispY_sq_over_beta_y=None, +): + """ + Calculates IBS growth rates (T_x, T_y, T_p) using model-specific scattering integrals and beam parameters. + This calculation is performed every 500 turns by default or at user-defined intervals specified by the update_turns variable, + depending on the lattice configuration. + + """ + + if model == "PS": + T_p = A * (vabq * (sigma_H**2 / sigma_p**2)) + T_x = A * (v1bq + (vabq * ((dispX_sq_over_beta_x * sigma_H**2) / (emitx)))) + T_y = A * (v1aq + (vabq * ((dispY_sq_over_beta_y * sigma_H**2) / (emity)))) + + elif model == "PM": + + T_p = A * (vabq * (sigma_H**2 / sigma_p**2)) + T_x = A * (v1bq + (vabq * ((H_x * sigma_H**2) / (emitx)))) + T_y = A * (v1aq + (vabq * ((H_y * sigma_H**2) / (emity)))) + + elif model == "Bane": + T_pp = (r_0**2 * N * C_log * sigma_H * gval * (beta_x * beta_y) ** (-1 / 4)) / ( + 16 * gamma**3 * emitx ** (3 / 4) * emity ** (3 / 4) * sigma_s * sigma_p**3 + ) + T_p = np.average(T_pp) + T_x = (sigma_p**2 * H_x * T_pp) / emitx + T_y = (sigma_p**2 * H_y * T_pp) / emity + + elif model == "CIMP": + K_a = ((np.log(q**2 / a**2) * g_ba) / a) + ((np.log(q**2 / b**2) * g_ab) / b) + + T_p = 2 * np.pi ** (3 / 2) * A * ((sigma_H**2 / sigma_p**2) * K_a) + T_x = ( + 2 + * np.pi ** (3 / 2) + * A + * ((-a * np.log(q**2 / a**2) * g_ba) + (((H_x * sigma_H**2) / emitx) * K_a)) + ) + T_y = ( + 2 + * np.pi ** (3 / 2) + * A + * ((-b * np.log(q**2 / b**2) * g_ab) + (((H_y * sigma_H**2) / emity) * K_a)) + ) + + T_x = np.average(T_x) + T_y = np.average(T_y) + T_p = np.average(T_p) + + if T_p <= 0: + T_p = 0 + if T_x <= 0: + T_x = 0 + if T_y <= 0: + T_y = 0 + + return T_x, T_y, T_p diff --git a/pyat/at/collective/ibs_element.py b/pyat/at/collective/ibs_element.py new file mode 100644 index 0000000000..b0d03c201b --- /dev/null +++ b/pyat/at/collective/ibs_element.py @@ -0,0 +1,165 @@ +""" +IBS element for Acclerator toolbox(AT) inspired from mbtrack2 +from Alexin Gamelin , Vadim Gubaidulin (SOLEIL) +Github : https://gitlab.synchrotron-soleil.fr/PA/collective-effects/mbtrack2 + +""" + +import numpy as np +import at +from enum import IntEnum +from ..lattice import Lattice, AtWarning +from at.lattice import Collective +from at.lattice.elements import _array +from at.lattice.utils import Refpts, uint32_refpts, make_copy +from at.physics import get_optics, avlinopt +from at.constants import clight, qe, e_mass +from typing import Sequence, Optional, Union +import warnings + +import scipy.integrate as quad +from scipy.constants import c, elementary_charge, epsilon_0, m_e +from scipy.special import hyp2f1 +from scipy.interpolate import interp1d +from at import Element + + +class IBSElement(Collective, Element): + """Class to generate an IBS element, inherits from Element""" + + default_pass = {False: "DriftPass", True: "DriftPass"} + # _conversions = dict(RFCavity._conversions, + # ) + + def __init__(self, family_name: str, length: float, ring: Lattice, **kwargs): + """ + Parameters: + ring (Lattice): Lattice object. + + model (str, optional): + Model to compute the growth rates. Options are "PS", "PM", "Bane", or "CIMP". + Default is "CIMP". + + n_points (int, optional): + Number of points at which optics is computed. + Default is 10000. + + n_bin (int, optional): + Number of bins for capturing non-uniform beam profiles. + Default is 100. + + get_opt (str, optional): + Method to calculate the optics around the lattice object. + Options: "mbtrack2_eqv", "markers", or "average". + Default is "markers". + + - "mbtrack2_eqv": Similar to mbtrack2 implementation, uses interpolation to get optics at equidistant points. + - "markers": Places markers at equidistant points and calculates optics at their positions using get_optics. + - "average": Places markers at equidistant points and averages optics using avlinopt. + + update_turns (int, optional): + Number of turns after which growth rates are recomputed. + Default is 500. + """ + + self.Turn_counter = 0 + self.ring = ring + self.revolution_frequency = ring.revolution_frequency + self.r_0 = (qe**2) / (4 * np.pi * epsilon_0 * c**2 * m_e) + self.gamma = (ring.energy * qe) / (c**2 * m_e) + self.beta = np.sqrt(1 - (e_mass**2 / ring.energy**2)) + self.T_x = 0 + self.T_y = 0 + self.T_p = 0 + + kwargs.setdefault("PassMethod", self.default_pass[True]) + + self.model = kwargs.pop("model", "CIMP") + self.get_opt = kwargs.pop("get_opt", "markers") + self.n_points = kwargs.pop( + "n_points", 10000 if self.model in ("mbtrack2_eqv", "markers") else 40000 + ) + self.update_turns = kwargs.pop("update_turns", 500) + self.n_bin = kwargs.pop("n_b", 100) + self.circumference = ring.get_s_pos(refpts=len(ring))[0] + super().__init__(family_name, **kwargs) + + def compute_optics_params(self): + + if self.get_opt == "mbtrack2_eqv": # most close to the one from mbtrack + + self.ring = self.ring.slice(slices=self.n_points) + s_points = np.linspace(0, self.circumference, self.n_points) + _, _, ld_6d = self.ring.disable_6d(copy=True).get_optics( + refpts=np.arange(len(self.ring) + 1) + ) + s_pos = self.ring.get_s_pos(refpts=np.arange(len(self.ring) + 1)) + bx_fun = interp1d(s_pos, ld_6d.beta[:, 0]) + by_fun = interp1d(s_pos, ld_6d.beta[:, 1]) + self.beta_x = bx_fun(s_points) + self.beta_y = by_fun(s_points) + ax_fun = interp1d(s_pos, ld_6d.alpha[:, 0]) + ay_fun = interp1d(s_pos, ld_6d.alpha[:, 1]) + self.alphaX = ax_fun(s_points) + self.alphaY = ay_fun(s_points) + + dx_fun = interp1d(s_pos, ld_6d.dispersion[:, 0]) + ddx_fun = interp1d(s_pos, ld_6d.dispersion[:, 1]) + dy_fun = interp1d(s_pos, ld_6d.dispersion[:, 2]) + ddy_fun = interp1d(s_pos, ld_6d.dispersion[:, 3]) + self.dispX = dx_fun(s_points) + self.disppX = ddx_fun(s_points) + self.dispY = dy_fun(s_points) + self.disppY = ddy_fun(s_points) + + elif self.get_opt == "markers": + + marker = at.Marker("ibs_marker") + marker_list = [marker] * (self.n_points) + s_points = np.linspace(0, self.circumference, self.n_points, endpoint=False) + self.ring = self.ring.sbreak(s_points, marker_list) + refpts = self.ring.get_uint32_index("ibs*") + _, _, ld_6d = self.ring.disable_6d(copy=True).get_optics(refpts=refpts) + self.beta_x = ld_6d.beta[:, 0] + self.beta_y = ld_6d.beta[:, 1] + self.alphaX = ld_6d.alpha[:, 0] + self.alphaY = ld_6d.alpha[:, 1] + self.dispX = ld_6d.dispersion[:, 0] + self.disppX = ld_6d.dispersion[:, 1] + self.dispY = ld_6d.dispersion[:, 2] + self.disppY = ld_6d.dispersion[:, 3] + + elif self.get_opt == "average": + + marker = at.Marker("ibs_marker") + marker_list = [marker] * (self.n_points) + s_points = np.linspace(0, self.circumference, self.n_points, endpoint=False) + self.ring = self.ring.sbreak(s_points, marker_list) + refpts = self.ring.get_uint32_index("ibs*") + elemdata, avebeta, _, avedisp, avespos, tune, chrom = self.ring.disable_6d( + copy=True + ).avlinopt(dp=0.0, refpts=refpts + 1) + self.beta_x = avebeta[:, 0] + self.beta_y = avebeta[:, 1] + self.alphaX = [e.alpha[0] for e in elemdata] + self.alphaY = [e.alpha[1] for e in elemdata] + self.dispX = avedisp[:, 0] + self.disppX = avedisp[:, 1] + self.dispY = avedisp[:, 2] + self.disppY = avedisp[:, 3] + + 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.dispX_sq = self.dispX**2 + self.dispX_sq_over_beta_x = self.dispX_sq / self.beta_x + self.dispY_sq = self.dispY**2 + self.dispY_sq_over_beta_y = self.dispY_sq / self.beta_y + + def clear_history(self): + pass diff --git a/pyat/at/integrators/pyIBSRadPass.py b/pyat/at/integrators/pyIBSRadPass.py new file mode 100644 index 0000000000..fdc6d28f34 --- /dev/null +++ b/pyat/at/integrators/pyIBSRadPass.py @@ -0,0 +1,241 @@ +""" +IBS passmethod for Acclerator toolbox(AT) inspired from mbtrack2 +from Alexin Gamelin , Vadim Gubaidulin (SOLEIL) +Github : https://gitlab.synchrotron-soleil.fr/PA/collective-effects/mbtrack2 + +""" + +import numpy as np +import at +import sys +import scipy.integrate as quad +from scipy.constants import c, elementary_charge, epsilon_0, m_e +from scipy.special import hyp2f1 +from scipy.interpolate import interp1d +from at.constants import qe +from at.constants import e_mass +from at.constants import clight +from at.collective import ibs + + +def kick( + bunch, n_bin, sigma_p, sigma_px, sigma_py, revolution_frequency, T_x, T_y, T_p +): + """ + Applies IBS-induced momentum kicks to the bunch every turn, scaled by growth rates, + with optional binning to capture non-uniform bunch profiles. + """ + + if n_bin > 1: + + data = bunch[5, :] / clight + bin_min = np.min(data) + bin_min = min(bin_min * 0.99, bin_min * 1.01) + bin_max = np.max(data) + bin_max = max(bin_max * 0.99, bin_max * 1.01) + + bins = np.linspace(bin_min, bin_max, n_bin + 1) + center = (bins[1:] + bins[:-1]) / 2 + sorted_index = np.searchsorted(bins, data, side="left") + sorted_index -= 1 + profile = np.bincount(sorted_index, minlength=n_bin) + + normalized_profile = profile / max(profile) + Rho = normalized_profile[sorted_index] + + else: + + Rho = 1.0 + + N_mp = len(bunch[0, :]) + Delta_pz = ( + sigma_p + * np.sqrt(np.sqrt(2) * T_p * (1 / revolution_frequency) * Rho) + * np.random.normal(size=N_mp) + ) + Delta_px = ( + sigma_px + * np.sqrt(np.sqrt(2) * T_x * (1 / revolution_frequency) * Rho) + * np.random.normal(size=N_mp) + ) + Delta_py = ( + sigma_py + * np.sqrt(np.sqrt(2) * T_y * (1 / revolution_frequency) * Rho) + * np.random.normal(size=N_mp) + ) + + bunch[1, :] += Delta_px + bunch[3, :] += Delta_py + bunch[4, :] += Delta_pz + + +def trackFunction(rin, elem=None): + + model = elem.model + revolution_frequency = elem.revolution_frequency + n_points = int(elem.n_points) + n_bin = elem.n_bin + gamma = elem.gamma + beta = elem.beta + r_0 = elem.r_0 + update_turns = int(elem.update_turns) + + ring = elem.ring + beta_x = elem.beta_x + beta_y = elem.beta_y + alphaX = elem.alphaX + alphaY = elem.alphaY + dispX = elem.dispX + disppX = elem.disppX + dispY = elem.dispY + disppY = elem.disppY + + dispX_sq_over_beta_x = elem.dispX_sq_over_beta_x + dispY_sq_over_beta_y = elem.dispY_sq_over_beta_y + + H_x = elem.H_x + H_y = elem.H_y + + N = ring.beam_current / (qe * revolution_frequency) + + sigma_py, sigma_px, sigma_p, sigma_s, d = ibs.initialize(bunch=rin, ring=ring) + + if elem.Turn_counter % update_turns == 0: + + a, b, q, sigma_H, emitx, emity = ibs.paramerters( + bunch=rin, + ring=ring, + model=model, + beta_x=beta_x, + beta_y=beta_y, + sigma_p=sigma_p, + dispX_sq_over_beta_x=dispX_sq_over_beta_x, + dispY_sq_over_beta_y=dispY_sq_over_beta_y, + H_x=H_x, + H_y=H_y, + gamma=gamma, + beta=beta, + d=d, + r_0=r_0, + ) + + if model == "Bane": + C_log = np.log(q**2 / a**2) + C_a = a / b + + elif model in ["PM", "PS", "CIMP"]: + A = (r_0**2 * N) / ( + 64 * np.pi**2 * beta**3 * gamma**4 * emitx * emity * sigma_s * sigma_p + ) + + if model in ["PM", "PS"]: + + vabq, v1aq, v1bq = ibs.scatter( + model=model, n_points=n_points, b=b, a=a, q=q + ) + + T_x, T_y, T_p = ibs.get_scatter_T( + vabq=vabq, + v1aq=v1aq, + v1bq=v1bq, + A=A, + model=model, + sigma_H=sigma_H, + dispX=dispX, + dispY=dispY, + beta_x=beta_x, + beta_y=beta_y, + r_0=r_0, + N=N, + gamma=gamma, + beta=beta, + emitx=emitx, + emity=emity, + sigma_s=sigma_s, + sigma_p=sigma_p, + H_x=H_x, + H_y=H_y, + b=b, + a=a, + q=q, + dispX_sq_over_beta_x=dispX_sq_over_beta_x, + dispY_sq_over_beta_y=dispY_sq_over_beta_y, + ) + + elif model == "Bane": + + gval = ibs.scatter(model=model, n_points=n_points, b=b, a=a, q=q, C_a=C_a) + + T_x, T_y, T_p = ibs.get_scatter_T( + gval=gval, + model=model, + sigma_H=sigma_H, + dispX=dispX, + dispY=dispY, + beta_x=beta_x, + beta_y=beta_y, + r_0=r_0, + N=N, + C_log=C_log, + gamma=gamma, + beta=beta, + emitx=emitx, + emity=emity, + sigma_s=sigma_s, + sigma_p=sigma_p, + H_x=H_x, + H_y=H_y, + b=b, + a=a, + q=q, + dispX_sq_over_beta_x=dispX_sq_over_beta_x, + dispY_sq_over_beta_y=dispY_sq_over_beta_y, + ) + + elif model == "CIMP": + + g_ab, g_ba = ibs.scatter(model=model, n_points=n_points, b=b, a=a, q=q) + + T_x, T_y, T_p = ibs.get_scatter_T( + g_ab=g_ab, + g_ba=g_ba, + A=A, + model=model, + sigma_H=sigma_H, + dispX=dispX, + dispY=dispY, + beta_x=beta_x, + beta_y=beta_y, + r_0=r_0, + N=N, + gamma=gamma, + beta=beta, + emitx=emitx, + emity=emity, + sigma_s=sigma_s, + sigma_p=sigma_p, + H_x=H_x, + H_y=H_y, + b=b, + a=a, + q=q, + dispY_sq_over_beta_y=dispY_sq_over_beta_y, + ) + + elem.T_x = T_x + elem.T_y = T_y + elem.T_p = T_p + + elem.Turn_counter += 1 + + kick( + rin, + n_bin, + sigma_p, + sigma_px, + sigma_py, + revolution_frequency, + elem.T_x, + elem.T_y, + elem.T_p, + ) diff --git a/pyat/test/test_ibs.py b/pyat/test/test_ibs.py new file mode 100644 index 0000000000..fdcdc0dff9 --- /dev/null +++ b/pyat/test/test_ibs.py @@ -0,0 +1,63 @@ +import numpy +import at +import pytest + + +def test_two_particles_for_one_turn(hmba_lattice): + + ring = hmba_lattice.enable_6d(copy=True) + bunch_current = 10e-3 + ring.set_fillpattern(1) + ring.set_beam_current(bunch_current) + optics = at.get_optics(ring.disable_6d(copy=True), refpts=0) + ibs_element = at.IBSElement( + "Ibs_element", + 0, + ring, + model="CIMP", + get_opt="markers", + n_points=100, + n_bin=1, + PassMethod="pyIBSRadPass", + ) + ibs_element.compute_optics_params() + ring.append(ibs_element) + + numpy.random.seed(42) + sigma_matrix = at.sigma_matrix( + betax=optics[0].beta[0], + betay=optics[0].beta[1], + alphax=optics[0].alpha[0], + alphay=optics[0].alpha[1], + emitx=141e-12, + emity=1e-12, + blength=0.003, + espread=0.001, + ) + rin = at.beam(2, sigma_matrix) + + rout_particle1_expected = numpy.array( + [ + -7.25596327e-07, + -7.64355188e-06, + 1.14092197e-06, + 7.33024711e-07, + 3.94815955e-04, + 6.27494565e-04, + ] + ) + rout_particle2_expected = numpy.array( + [ + -2.80118387e-05, + -3.36321436e-06, + -1.19125693e-06, + -2.81718065e-07, + 6.31241787e-04, + -5.73637434e-03, + ] + ) + + ring.track(rin) + + numpy.testing.assert_almost_equal(rin[:, 0], rout_particle1_expected, decimal=6) + numpy.testing.assert_almost_equal(rin[:, 1], rout_particle2_expected, decimal=6) From b7be90a029f25e2eceed74039932cf06bc4ddf7b Mon Sep 17 00:00:00 2001 From: Ritu Bharadwaj Date: Wed, 16 Jul 2025 16:34:09 +0200 Subject: [PATCH 2/3] This is the test file as suggested by Lee Carver --- pyat/examples/CollectiveEffects/IBS_run.py | 200 +++++++++++++++++++++ 1 file changed, 200 insertions(+) create mode 100644 pyat/examples/CollectiveEffects/IBS_run.py diff --git a/pyat/examples/CollectiveEffects/IBS_run.py b/pyat/examples/CollectiveEffects/IBS_run.py new file mode 100644 index 0000000000..763ee5d45a --- /dev/null +++ b/pyat/examples/CollectiveEffects/IBS_run.py @@ -0,0 +1,200 @@ +""" +Example python script to test the IBS effect +""" + +import at +import numpy as np +import matplotlib.pyplot as plt +from at.constants import qe, clight +from at import Wake, WakeType, WakeComponent, WakeElement +import pickle + + +class IBS_run: + + def __init__(self, bunch_current_mA, emity_pm): + + self.bunch_current = float(bunch_current_mA) * 1e-3 # Bunch current [A] + self.emity = float(emity_pm) * 1e-12 # Vertical emittance [m] + self.emitx = 141e-12 + self.model = "CIMP" # model to compute the growth rates. + self.get_opt = "markers" # method to compute the lattice optics. + self.update_turns = ( + 500 # number of turns after which the grwoth rates for IBS are computed. + ) + self.nturns = 10000 # total number of turns for tracking + + self.ring, self.mcf = self._define_ring() + + def _define_ring(self): + """ + This method is used to load the lattice and set it by defining the number of bunches , bunch current and including + radiation damping and quantum diffusion element. + """ + + ring_hmba = at.load_mat('../../../machine_data/hmba.mat') # single cell for the EBS lattice + ring = ring_hmba.repeat(32) # cell repeated 32 times to get the complete ring + ring.disable_6d() + mcf = at.get_mcf(ring) + ring.enable_6d() + ring.set_rf_voltage(5.5e6) # voltage of the RF cavity (5.5 MV for EBS ring) + ring.set_fillpattern(1) + ring.set_beam_current(self.bunch_current) + ring.set_cavity_phase() + qdelem = at.gen_quantdiff_elem(ring) # adds the quantum diffusion element + ring.append(qdelem) + return ring, mcf + + def _build_simple_ring(self, ibs_elem, params, optics, chromaticity): + """ + This method is to define a simple ring (with an RF cavity element ,a 6x6 linear transfer map with no radiation damping, + a detuning and chromaticity element, a simple radiation damping element, and a simplified quantum diffusion element + which contains equilibrium emittance) for fast tracking. + + Further, an ibs element is also added to model the intrabeam scattering effect. + + """ + sring = at.simple_ring( + energy=self.ring.energy, + circumference=self.ring.circumference, + harmonic_number=self.ring.harmonic_number, + Qx=params.tunes6[0], + Qy=params.tunes6[1], + Vrf=self.ring.rf_voltage, + alpha=self.mcf, + betax=optics[0].beta[0], + betay=optics[0].beta[1], + alphax=optics[0].alpha[0], + alphay=optics[0].alpha[1], + dispx=optics[0].dispersion[0], + dispxp=optics[0].dispersion[1], + dispy=optics[0].dispersion[2], + dispyp=optics[0].dispersion[3], + Qpx=chromaticity[0], + Qpy=chromaticity[1], + emitx=params.emittances[0], + emity=self.emity, + espread=params.sigma_e, + taux=params.Tau[0] * self.ring.revolution_frequency, + tauy=params.Tau[1] * self.ring.revolution_frequency, + tauz=params.Tau[2] * self.ring.revolution_frequency, + U0=self.ring.energy_loss, + name="esrf", + particle="electron", + TimeLag=0.0, + ) + + sring.enable_6d() + sring.set_fillpattern(1) + sring.set_beam_current(self.bunch_current) + sring.set_cavity_phase() + + ring_ibs = sring.deepcopy() + ring_ibs.append(ibs_elem) + return ring_ibs + + def _generate_bunch(self, params, optics): + """ + Generate a particle bunch using the sigma matrix based on given optics and beam parameters. + + """ + sigma_matrix = at.sigma_matrix( + betax=optics[0].beta[0], + betay=optics[0].beta[1], + alphax=optics[0].alpha[0], + alphay=optics[0].alpha[1], + emitx=self.emitx, + emity=self.emity, + blength=params.sigma_l, + espread=params.sigma_e, + ) + return at.beam( + 1000, sigma_matrix + ) # 400000 is the number of particles in the bunch. + + def _track(self, ring_ibs, bunch): + """ + Track the bunch through the ring including radiation damping, quantum diffusion, + and intrabeam scattering effects. + + Calculates bunch length, energy spread, and emittances every turn and store the values in lists. + + """ + + emit_ibs = [] + blength_ibs = [] + espread_ibs = [] + + for _ in range(self.nturns): + ring_ibs.track(bunch, refpts=None, in_place=True) + sigma = at.sigma_matrix(beam=bunch) + eigs, _ = np.linalg.eig(sigma @ at.jmat(3)) + emits = np.sort(np.abs(eigs)[::2]) + emit_ibs.append([emits[1], emits[0], emits[2]]) + blength_ibs.append(np.std(bunch[5, :])) + espread_ibs.append(np.std(bunch[4, :])) + + return blength_ibs, espread_ibs, emit_ibs + + def run_simulation(self): + """ + This method defines the IBS element with the choice of model ('CIMP', 'Bane', 'PM', 'PS'), + the method to calculate optics ('mbtrack_eqv', 'markers', 'average'), and update_turns + (i.e., how many turns between IBS growth rate computations) etc. These parameters are defined in the __init__ method. + + It then calls the previously defined functions for ring and bunch definition, and performs tracking. + + """ + + ibs_elem = at.IBSElement( + "Ibs", + 0, + self.ring, + model=self.model, + get_opt=self.get_opt, + update_turns=self.update_turns, + PassMethod="pyIBSRadPass", + ) + ibs_elem.compute_optics_params() + params = self.ring.envelope_parameters() + _, rd, _ = at.get_optics(self.ring, get_chrom=True) + chromaticity = rd.chromaticity + optics = at.get_optics(self.ring.disable_6d(copy=True)) + + ring_ibs = self._build_simple_ring(ibs_elem, params, optics, chromaticity) + bunch = self._generate_bunch(params, optics) + + return self._track(ring_ibs, bunch) + + +############################################################################################################################ + +## Run simulation with given bunch current and horizontal emittance values. +## Adjust these parameters to study how IBS effects vary. +## +## Outputs bunch length, energy spread, and emittance, then saves results to a pickle file. + + +def launch(): + bunch_current = 10.0 # mA + emity = 1.0 # pm + + sim = IBS_run(bunch_current, emity) + blength, espread, emit = sim.run_simulation() + + outdict = { + "bunch_current_mA": bunch_current, + "emity_pm": emity, + "bunch_length": blength, + "energy_spread": espread, + "emittances": emit, + } + + with open("./IBS_run_output.pkl", "wb") as f: + pickle.dump(outdict, f) + + print("Simulation complete. Output saved to IBS_run_output.pkl") + + +if __name__ == "__main__": + launch() From 8ab0f7369f9f00a41a0a4b0aecdd272b491d0d42 Mon Sep 17 00:00:00 2001 From: Ritu Bharadwaj Date: Wed, 16 Jul 2025 16:37:42 +0200 Subject: [PATCH 3/3] I reduced the number of particles in the bunch for the faster tracking. --- pyat/examples/CollectiveEffects/IBS_run.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyat/examples/CollectiveEffects/IBS_run.py b/pyat/examples/CollectiveEffects/IBS_run.py index 763ee5d45a..67e2981afb 100644 --- a/pyat/examples/CollectiveEffects/IBS_run.py +++ b/pyat/examples/CollectiveEffects/IBS_run.py @@ -110,7 +110,7 @@ def _generate_bunch(self, params, optics): ) return at.beam( 1000, sigma_matrix - ) # 400000 is the number of particles in the bunch. + ) # 1000 is the number of particles in the bunch. Increase the number of particles for more reliable results. def _track(self, ring_ibs, bunch): """