Skip to content

Commit 6d2a2ff

Browse files
giadarolgiadarol
authored
Merge pull request #157 from SebastienJoly/feature/equilibrium_emittance
Equilibrium emittance in presence of SR and IBS
2 parents a477233 + 71ddab6 commit 6d2a2ff

File tree

8 files changed

+1123
-10
lines changed

8 files changed

+1123
-10
lines changed

‎contributors_and_credits.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,4 +13,8 @@ CERN contributors:
1313
- F. Soubelet
1414
- E. Waagaard
1515

16+
External contributors:
17+
- S. Joly (Helmholtz-Zentrum, Berlin) - Development of IBS + SR equilibrium calculation
18+
19+
1620
The Particle In Cell implementation is largely based on the PyPIC package (https://github.com/pycomplete/pypic)

‎examples/005_ibs/005_steady_state_emittances_coupling.py‎

Lines changed: 130 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,130 @@
1+
# copyright ################################# #
2+
# This file is part of the Xfields Package. #
3+
# Copyright (c) CERN, 2021. #
4+
# ########################################### #
5+
import matplotlib.pyplot as plt
6+
import xtrack as xt
7+
from scipy.constants import e
8+
9+
##########################
10+
# Load xt.Line from file #
11+
##########################
12+
13+
fname_line_particles = "../../../xtrack/test_data/bessy3/bessy3.json"
14+
line = xt.Line.from_json(fname_line_particles) # has particle_ref
15+
line.build_tracker()
16+
17+
########################
18+
# Twiss with Radiation #
19+
########################
20+
21+
# We need to Twiss with Synchrotron Radiation enabled to obtain
22+
# the SR equilibrium emittances and damping constants
23+
24+
line.matrix_stability_tol = 1e-2
25+
line.configure_radiation(model="mean")
26+
line.compensate_radiation_energy_loss()
27+
tw = line.twiss(eneloss_and_damping=True)
28+
29+
######################################
30+
# Steady-State Emittance Calculation #
31+
######################################
32+
33+
bunch_intensity = 1e-9 / e # 1C bunch charge
34+
emittance_coupling_factor = 1 # round beam
35+
36+
# If not providing starting emittances, the function will
37+
# default to the SR equilibrium emittances from the TwissTable
38+
39+
result = tw.get_ibs_and_synrad_emittance_evolution(
40+
formalism="nagaitsev", # can also be "bjorken-mtingwa"
41+
total_beam_intensity=bunch_intensity,
42+
# gemitt_x=..., # defaults to tw.eq_gemitt_x
43+
# gemitt_y=..., # defaults to tw.eq_gemitt_x
44+
# gemitt_zeta=..., # defaults to tw.eq_gemitt_zeta
45+
emittance_coupling_factor=emittance_coupling_factor,
46+
emittance_constraint="coupling",
47+
)
48+
49+
# The returned object is an xtrack Table
50+
print(result)
51+
52+
# Table: 1089 rows, 9 cols
53+
# time gemitt_x nemitt_x gemitt_y nemitt_y gemitt_zeta nemitt_zeta sigma_zeta sigma_delta Kx Ky Kz
54+
# 1.2104374139405318e-06 6.8355e-11 3.34418e-07 6.8355e-11 3.34418e-07 3.38167e-06 0.0581296 0.00344698 0.000981052 58.4838 -0.0165806 17.0318
55+
# 9.306647341196751e-05 6.87219e-11 3.36214e-07 6.87219e-11 3.36214e-07 3.39225e-06 0.0583114 0.00345237 0.000982586 58.4736 -0.0165775 17.0292
56+
# 0.0001849225094099945 6.90808e-11 3.37969e-07 6.90808e-11 3.37969e-07 3.40267e-06 0.0584907 0.00345767 0.000984095 57.7058 -0.0163509 16.8361
57+
# ...
58+
# 0.0997568655312697 8.4492e-11 4.13367e-07 8.4492e-11 4.13367e-07 4.51774e-06 0.0776582 0.00398413 0.00113393 29.799 -0.00861232 8.14787
59+
# 0.09984872156726772 8.4492e-11 4.13367e-07 8.4492e-11 4.13367e-07 4.51774e-06 0.0776583 0.00398413 0.00113393 29.799 -0.00861232 8.14786
60+
# 0.09994057760326575 8.44919e-11 4.13367e-07 8.44919e-11 4.13367e-07 4.51775e-06 0.0776584 0.00398414 0.00113393 29.799 -0.00861232 8.14784
61+
62+
######################################
63+
# Comparison with analytical results #
64+
######################################
65+
66+
# These are analytical estimate (from the last step's IBS growth rates)
67+
# The factor below is to be respected with the coupling constraint
68+
factor = 1 + emittance_coupling_factor * (tw.partition_numbers[1] / tw.partition_numbers[0])
69+
analytical_x = result.gemitt_x[0] / (1 - result.Kx[-1] / (tw.damping_constants_s[0] * factor))
70+
analytical_y = result.gemitt_y[0] / (1 - result.Kx[-1] / (tw.damping_constants_s[0] * factor))
71+
analytical_z = result.gemitt_zeta[0] / (1 - result.Kz[-1] / (tw.damping_constants_s[2]))
72+
73+
print()
74+
print("Emittance Constraint: Coupling")
75+
print("Horizontal steady-state emittance:")
76+
print("---------------------------------")
77+
print(f"Analytical: {analytical_x}")
78+
print(f"ODE: {result.eq_sr_ibs_gemitt_x}")
79+
print("Vertical steady-state emittance:")
80+
print("-------------------------------")
81+
print(f"Analytical: {analytical_y}")
82+
print(f"ODE: {result.eq_sr_ibs_gemitt_y}")
83+
print("Longitudinal steady-state emittance:")
84+
print("-----------------------------------")
85+
print(f"Analytical: {analytical_z}")
86+
print(f"ODE: {result.eq_sr_ibs_gemitt_zeta}")
87+
88+
# Emittance Constraint: Coupling
89+
# Horizontal steady-state emittance:
90+
# ---------------------------------
91+
# Analytical: 8.450195693021898e-11
92+
# ODE: 8.449194372183254e-11
93+
# Vertical steady-state emittance:
94+
# -------------------------------
95+
# Analytical: 8.450195693021898e-11
96+
# ODE: 8.449194372183254e-11
97+
# Longitudinal steady-state emittance:
98+
# -----------------------------------
99+
# Analytical: 4.518939710729963e-06
100+
# ODE: 4.5177482475985255e-06
101+
102+
# The results from the table can easily be plotted to view
103+
# at the evolution of various parameters across time steps
104+
105+
#!end-doc-part
106+
# fmt: off
107+
108+
fig, (ax0, ax1) = plt.subplots(2, 1, sharex=True, layout="constrained")
109+
110+
(l1,) = ax0.plot(result.time * 1e3, result.gemitt_x * 1e12, ls="-", label=r"$\tilde{\varepsilon}_x$")
111+
(l2,) = ax0.plot(result.time * 1e3, result.gemitt_y * 1e12, ls="--", label=r"$\tilde{\varepsilon}_y$")
112+
l4 = ax0.axhline(analytical_x * 1e12, color="C0", ls="-.", label=r"Analytical $\varepsilon_{x}^{eq}$")
113+
l5 = ax0.axhline(analytical_y * 1e12, color="C1", ls="-.", label=r"Analytical $\varepsilon_{y}^{eq}$")
114+
ax0b = ax0.twinx()
115+
(l3,) = ax0b.plot(result.time * 1e3, result.gemitt_zeta * 1e6, color="C2", label=r"$\varepsilon_z$")
116+
l6 = ax0b.axhline(analytical_z * 1e6, color="C2", ls="-.", label=r"Analytical $\varepsilon_{\zeta}^{eq}$")
117+
ax0.legend(handles=[l1, l2, l3, l4, l5], ncols=2)
118+
119+
ax1.plot(result.time * 1e3, result.Kx, label=r"$K_{x}^{IBS}$")
120+
ax1.plot(result.time * 1e3, result.Ky, label=r"$K_{y}^{IBS}$")
121+
ax1.plot(result.time * 1e3, result.Kz, label=r"$K_{z}^{IBS}$")
122+
ax1.legend()
123+
124+
ax1.set_xlabel("Time [ms]")
125+
ax0.set_ylabel(r"$\tilde{\varepsilon}_{x,y}$ [pm.rad]")
126+
ax0b.set_ylabel(r"$\varepsilon_{\zeta}$ [m]")
127+
ax1.set_ylabel(r"$K_{x,y,z}^{IBS}$ [$s^{-1}$]")
128+
fig.align_ylabels((ax0, ax1))
129+
plt.tight_layout()
130+
plt.show()

‎examples/005_ibs/006_steady_state_emittances_excitation.py‎

Lines changed: 99 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,99 @@
1+
# copyright ################################# #
2+
# This file is part of the Xfields Package. #
3+
# Copyright (c) CERN, 2021. #
4+
# ########################################### #
5+
import matplotlib.pyplot as plt
6+
import xtrack as xt
7+
from scipy.constants import e
8+
9+
##########################
10+
# Load xt.Line from file #
11+
##########################
12+
13+
fname_line_particles = "../../../xtrack/test_data/bessy3/bessy3.json"
14+
line = xt.Line.from_json(fname_line_particles) # has particle_ref
15+
line.build_tracker()
16+
17+
########################
18+
# Twiss with Radiation #
19+
########################
20+
21+
# We need to Twiss with Synchrotron Radiation enabled to obtain
22+
# the SR equilibrium emittances and damping constants
23+
24+
line.matrix_stability_tol = 1e-2
25+
line.configure_radiation(model="mean")
26+
line.compensate_radiation_energy_loss()
27+
tw = line.twiss(eneloss_and_damping=True)
28+
29+
######################################
30+
# Steady-State Emittance Calculation #
31+
######################################
32+
33+
bunch_intensity = 1e-9 / e # 1C bunch charge
34+
emittance_coupling_factor = 0.5 # for excitation this time
35+
36+
# One can provide specific values for starting emittances,
37+
# but we need to ensure they respect the emittance coupling
38+
# contraint we want to enforce
39+
gemitt_x = 1.1 * tw.eq_gemitt_x # larger horizontal emittance
40+
gemitt_y = emittance_coupling_factor * gemitt_x # enforce the constraint
41+
42+
# One can overwrite sigma_zeta / sigma_delta (larger
43+
# values from potential well distortion for example)
44+
overwrite_sigma_zeta = 1.2 * (tw.eq_gemitt_zeta * tw.bets0) ** 0.5 # larger sigma_zeta
45+
overwrite_sigma_delta = 1.2 * (tw.eq_gemitt_zeta / tw.bets0) ** 0.5 # larger sigma_delta
46+
47+
# A specific time step or relative tolerance for convergence can also be provided.
48+
result = tw.get_ibs_and_synrad_emittance_evolution(
49+
formalism="nagaitsev", # can also be "bjorken-mtingwa"
50+
total_beam_intensity=bunch_intensity,
51+
gemitt_x=gemitt_x, # provided explicitely
52+
gemitt_y=gemitt_y, # provided explicitely
53+
overwrite_sigma_zeta=overwrite_sigma_zeta, # will recompute gemitt_zeta
54+
overwrite_sigma_delta=overwrite_sigma_delta, # will recompute gemitt_zeta
55+
emittance_coupling_factor=emittance_coupling_factor,
56+
emittance_constraint="excitation",
57+
)
58+
59+
# The returned object is a Table
60+
print(result)
61+
62+
# Table: 989 rows, 9 cols
63+
# time gemitt_x nemitt_x gemitt_y nemitt_y gemitt_zeta nemitt_zeta sigma_zeta sigma_delta Kx Ky Kz
64+
# 1.2104374139405318e-06 1.07706e-10 5.26938e-07 5.3853e-11 2.63469e-07 4.05791e-06 0.069754 0.00413633 0.000981042 32.1477 -0.00949839 13.5932
65+
# 9.306647341196751e-05 1.08146e-10 5.29091e-07 5.4073e-11 2.64546e-07 4.06402e-06 0.069859 0.00413944 0.00098178 32.1437 -0.00949715 13.5919
66+
# 0.0001849225094099945 1.08574e-10 5.31185e-07 5.4287e-11 2.65592e-07 4.07004e-06 0.0699624 0.00414251 0.000982507 31.844 -0.00940374 13.4942
67+
# ...
68+
# 0.09057126193146717 1.22535e-10 5.99488e-07 6.12675e-11 2.99744e-07 4.70398e-06 0.0808596 0.00445345 0.00105626 21.8776 -0.00648841 9.10733
69+
# 0.09066311796746519 1.22535e-10 5.99487e-07 6.12675e-11 2.99744e-07 4.70398e-06 0.0808597 0.00445345 0.00105626 21.8776 -0.00648841 9.10732
70+
# 0.09075497400346322 1.22535e-10 5.99487e-07 6.12675e-11 2.99744e-07 4.70399e-06 0.0808598 0.00445345 0.00105626 21.8776 -0.00648841 9.10731
71+
72+
73+
# The results from the table can easily be plotted . One notices
74+
# the transverse coupling factor is respected at all steps.
75+
76+
#!end-doc-part
77+
# fmt: off
78+
79+
fig, (ax0, ax1) = plt.subplots(2, 1, sharex=True, layout="constrained")
80+
81+
(l1,) = ax0.plot(result.time * 1e3, result.gemitt_x * 1e12, ls="-", label=r"$\tilde{\varepsilon}_x$")
82+
(l2,) = ax0.plot(result.time * 1e3, result.gemitt_y * 1e12, ls="-", label=r"$\tilde{\varepsilon}_y$")
83+
(l3,) = ax0.plot(result.time * 1e3, result.gemitt_y / emittance_coupling_factor * 1e12, ls=":", c="C1", label=r"$\tilde{\varepsilon}_y$ / factor")
84+
ax0b = ax0.twinx()
85+
(l4,) = ax0b.plot(result.time * 1e3, result.gemitt_zeta * 1e6, color="C2", label=r"$\varepsilon_z$")
86+
ax0.legend(handles=[l1, l2, l3, l4], ncols=2)
87+
88+
ax1.plot(result.time * 1e3, result.Kx, label=r"$K_{x}^{IBS}$")
89+
ax1.plot(result.time * 1e3, result.Ky, label=r"$K_{y}^{IBS}$")
90+
ax1.plot(result.time * 1e3, result.Kz, label=r"$K_{z}^{IBS}$")
91+
ax1.legend()
92+
93+
ax1.set_xlabel("Time [ms]")
94+
ax0.set_ylabel(r"$\tilde{\varepsilon}_{x,y}$ [pm.rad]")
95+
ax0b.set_ylabel(r"$\varepsilon_{\zeta}$ [m]")
96+
ax1.set_ylabel(r"$K_{x,y,z}^{IBS}$ [$s^{-1}$]")
97+
fig.align_ylabels((ax0, ax1))
98+
plt.tight_layout()
99+
plt.show()

‎tests/ibs_conftest.py‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,23 +1,23 @@
11
from __future__ import annotations
22

33
from pathlib import Path
4-
from typing import Tuple
54

65
import xtrack as xt
76
from cpymad.madx import Madx
87

98
# /!\ This assumes xtrack repo is sitting next to xfields repo
109
XTRACK_TEST_DATA = Path(__file__).parent.parent.parent / "xtrack" / "test_data/"
1110

11+
1212
def set_madx_beam_parameters(
1313
madx: Madx,
1414
total_beam_intensity: int,
15-
gemitt_x: float = None,
16-
nemitt_x: float = None,
17-
gemitt_y: float = None,
18-
nemitt_y: float = None,
19-
sigma_delta: float = None,
20-
bunch_length: float = None,
15+
gemitt_x: float | None = None,
16+
nemitt_x: float | None = None,
17+
gemitt_y: float | None = None,
18+
nemitt_y: float | None = None,
19+
sigma_delta: float | None = None,
20+
bunch_length: float | None = None,
2121
bunched: bool = True,
2222
) -> None:
2323
"""
@@ -82,7 +82,7 @@ def set_madx_beam_parameters(
8282
madx.sequence[seq_name].beam.bunched = bunched # set if the beam is bunched
8383

8484

85-
def get_madx_ibs_growth_rates(madx: Madx) -> Tuple[float, float, float]:
85+
def get_madx_ibs_growth_rates(madx: Madx) -> tuple[float, float, float]:
8686
"""
8787
Calls the IBS module then return horizontal, vertical and longitudinal
8888
emittance growth rates. A Twiss call is done to make sure it is centered.
@@ -138,7 +138,7 @@ def get_ref_particle_from_madx_beam(madx: Madx) -> xt.Particles:
138138
return xt.Particles(q0=q0, mass0=mass0, gamma0=gamma0)
139139

140140

141-
def get_parameters_from_madx_beam(madx: Madx) -> Tuple[float, float, float, float, float]:
141+
def get_parameters_from_madx_beam(madx: Madx) -> tuple[float, float, float, float, float]:
142142
"""
143143
Get beam intensity, transverse emittances, momentum spread and
144144
bunch length from the MAD-X's beam parameters. A Twiss call is

0 commit comments

Comments
 (0)