# -*- coding: utf-8 -*-
"""
Periodic transient beam loading (PTBL) (also called coupled bunch mode l=±1 instability).
Three backends:
- ptbl_he: He's perturbation criterion [1].
- ptbl_alves: Gaussian longitudinal mode-coupling instabilty method [2].
- coupled_bunch_mode1: Bosch [1] l=±1 coupled-bunch dipole mode [3].
References
----------
[1] He, Tianlong. "Novel perturbation method for judging the stability of the
equilibrium solution in the presence of passive harmonic cavities." Physical
Review Accelerators and Beams 25.9 (2022): 094402.
[2] Alves, Murilo B. "Theoretical models for longitudinal coupled-bunch instabilities
driven by harmonic cavities in electron storage rings." Physical Review Accelerators
and Beams 28.3 (2025): 034401.
[3] He, Tianlong. "Novel perturbation method for judging the stability of the
equilibrium solution in the presence of passive harmonic cavities." Physical
Review Accelerators and Beams 25.9 (2022): 094402.
"""
from __future__ import annotations
from math import factorial
from typing import TYPE_CHECKING
import numpy as np
from scipy.constants import c
from scipy.optimize import root_scalar
from albums.physics.robinson import phi_pm
if TYPE_CHECKING: # pragma: no cover - hints only
from mbtrack2 import CavityResonator, Synchrotron
from pycolleff.longitudinal_equilibrium import LongitudinalEquilibrium
[docs]def ptbl_he(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
PHI: np.ndarray,
) -> bool:
"""PTBL stability test via He's perturbation criterion (eta > 1).
Parameters
----------
ring : Synchrotron
Ring parameters (uses h, T0).
cavity_list : list of CavityResonator
Solved cavities (main then harmonic).
I0 : float
Beam current in A.
F : numpy.ndarray
Amplitude form factor per cavity (F[1] is the harmonic cavity).
PHI : numpy.ndarray
Phase form factor per cavity (PHI[1] is the harmonic cavity).
Returns
-------
bool
True if PTBL-unstable (eta > 1).
References
----------
He [1].
"""
m = ring.h
MC = cavity_list[0]
HC = cavity_list[1]
a = np.exp(-HC.wr * ring.T0 / (2 * HC.Q))
theta = HC.detune * 2 * np.pi * ring.T0
dtheta = np.arcsin(
(1 - a) * np.cos(theta / 2) / (np.sqrt(1 + a**2 - 2 * a * np.cos(theta)))
)
k = np.arange(1, m)
d_k = np.exp(-1 * HC.wr * ring.T0 * (k - 1) / (2 * HC.Q * m))
theta_k = theta / 2 + dtheta - (k - 1) / m * theta
eps_k = 1 - np.sin(np.pi / 2 - k * 2 * np.pi / m)
num = np.sum(eps_k * d_k * np.cos(theta_k))
f = num / (m * np.sqrt(1 + a**2 - 2 * a * np.cos(theta)))
eta = (
2
* np.pi
* HC.m**2
* F[1]
* ring.h
* I0
* HC.Rs
/ HC.Q
* f
/ (MC.Vc * np.sin(MC.theta - PHI[1] / HC.m))
)
return bool(eta > 1)
[docs]def coupled_bunch_mode1_damping_rate(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
Omega: float,
mu: int,
bunch_length: float,
) -> float:
"""Damping rate of the l=±1 coupled-bunch mode. Eq. (15) in [3]."""
coef1 = ring.ac * ring.omega1 * I0 / (Omega * ring.E0 * ring.T0)
coef2 = (
mu * (ring.omega1 * bunch_length) ** (2 * mu - 2) / (2**mu * factorial(mu - 1))
)
sum_val = 0
for i, cavity in enumerate(cavity_list):
diff = (
np.cos(phi_pm(ring, cavity, Omega, "-")) ** 2
- np.cos(phi_pm(ring, cavity, Omega, "+")) ** 2
)
sum_val += cavity.m ** (2 * mu - 1) * F[i] ** 2 * cavity.RL * diff
return coef1 * coef2 * sum_val
[docs]def coupled_bunch_mode1(
ring: Synchrotron,
cavity_list: list[CavityResonator],
I0: float,
F: np.ndarray,
omega_r: float,
bunch_length: float,
landau_threshold: np.ndarray,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Coupled-bunch dipole mode l=±1 (Bosch mode-1 criterion).
Parameters
----------
ring : Synchrotron
Ring parameters.
cavity_list : list of CavityResonator
Solved cavities.
I0 : float
Beam current in A.
F : numpy.ndarray
Amplitude form factor per cavity.
omega_r : float
Synchrotron angular frequency.
bunch_length : float
RMS bunch length in [s].
landau_threshold : numpy.ndarray
Per-mode thresholds from
landau_threshold.
Returns
-------
coupled_mode1 : numpy.ndarray of bool, shape (2,)
Instability flag for the +1 and -1 modes.
modes : numpy.ndarray, shape (2,)
Coherent angular frequency for each mode.
converged : numpy.ndarray of bool, shape (2,)
Root-finder convergence per mode.
References
----------
Eq. (15) in [3].
"""
modes = ["+1", "-1"]
mu = 1
coupled_mode1 = np.zeros(len(modes), dtype=bool)
mode1_modes = np.zeros(len(modes), dtype=float)
converged = np.zeros(len(modes), dtype=bool)
for j, mode in enumerate(modes):
if mode == "+1":
added = ring.omega0
elif mode == "-1":
added = -ring.omega0
Omega0 = mu * omega_r
def f(Omega):
coef1 = ring.ac * ring.omega1 * I0 / (ring.E0 * ring.T0)
coef2 = (
mu
* (ring.omega1 * bunch_length) ** (2 * mu - 2)
/ (2**mu * factorial(mu - 1))
)
sum_val = 0
for i, cavity in enumerate(cavity_list):
brackets = np.sin(
2 * phi_pm(ring, cavity, Omega + added, "-")
) + np.sin(2 * phi_pm(ring, cavity, Omega + added, "+"))
sum_val += cavity.m ** (2 * mu - 1) * cavity.RL * F[i] ** 2 * brackets
inner = (mu * omega_r) ** 2 - coef1 * coef2 * sum_val
return Omega - np.sqrt(inner)
sol = root_scalar(f, x0=Omega0, xtol=1e-6)
if sol.converged:
Omega = sol.root
mode1_modes[j] = Omega
converged[j] = True
alpha_r = coupled_bunch_mode1_damping_rate(
ring, cavity_list, I0, F, Omega + added, mu, bunch_length
)
alpha_r_incl_rad_damp = alpha_r + mu / ring.tau[2]
delta_Omega = Omega - 1j * alpha_r - mu * omega_r
if alpha_r_incl_rad_damp < 0:
if np.abs(delta_Omega) > landau_threshold[mu - 1]:
coupled_mode1[j] = True
else:
converged[j] = False
return coupled_mode1, mode1_modes, converged
[docs]def prepare_alves_le(ring: Synchrotron, LE: LongitudinalEquilibrium) -> None:
"""Set the synchrotron tune from the bunch length and select ImpedanceDFT."""
from pycolleff.longitudinal_equilibrium import ImpedanceSource
sync_freq = ring.sigma_delta * ring.eta() * c
sync_freq /= LE.ring.bunlen * 2 * np.pi
LE.ring.sync_tune = sync_freq / ring.f0
for idx, _ in enumerate(LE.impedance_sources):
LE.impedance_sources[idx].calc_method = ImpedanceSource.Methods.ImpedanceDFT
# Workaround: some ImpedanceSource objects carry loop_ctrl_transfer as an int
# flag instead of a callable; make it a function so get_impedance won't crash.
for imp in LE.impedance_sources:
if hasattr(imp, "loop_ctrl_transfer") and not callable(imp.loop_ctrl_transfer):
imp.loop_ctrl_transfer = lambda w, wc: 0 * np.asarray(w)
[docs]def ptbl_alves(
ring: Synchrotron, LE: LongitudinalEquilibrium, max_azi: int = 10, max_rad: int = 10
) -> tuple[bool, float]:
"""PTBL via Gaussian longitudinal mode-coupling instability method (coupled-bunch mode 1).
Also flag the usual coupled-bunch instability of mode 1 as PTBL.
Parameters
----------
ring : Synchrotron
Ring parameters (uses omega1, tau).
LE : LongitudinalEquilibrium
Solved pycolleff equilibrium handle (Equilibrium.backend_LE).
max_azi : int, optional
Maximum azimuthal mode number (default 10).
max_rad : int, optional
Maximum radial mode number (default 10).
Returns
-------
unstable : bool
True if the largest net growth rate is positive.
growth : float
Net growth rate (eigenvalue imaginary part minus radiation damping).
References
----------
Alves [2].
"""
prepare_alves_le(ring, LE)
eigenfreq, *_ = LE.calc_mode_coupling(
w=[-10 * ring.omega1, +10 * ring.omega1],
cbmode=1,
max_azi=max_azi,
max_rad=max_rad,
use_fokker=False,
delete_m0=True,
delete_m0k0=True,
reduced=True,
)
idx = np.argmax(eigenfreq.imag)
grate = eigenfreq.imag[idx] - 1 / ring.tau[2]
return bool(grate > 0), grate