"<h2>Simulation of Full CERN PS Acceleration Ramp</h2>\n",
"<h3>CERN PS Machine parameters</h3>"
]
},
{
"cell_type": "markdown",
"id": "taken-reminder",
"metadata": {},
"source": [
"The CERN Proton Synchrotron\n",
"- has a circumference of 2π·100m\n",
"- takes protons from the PS Booster at a kinetic energy of 2GeV corresponding to a $\\gamma=3.13$\n",
"- injects with 50kV of rf voltage, up to 200kV for ramp\n",
"- runs at harmonic $h=7$\n",
"- has a momentum compaction factor of $\\alpha_c=0.027$\n",
"- typical acceleration rate of (up to) $\\dot{B}=2$ T/s, the bending radius is $\\rho=70.08$ m"
]
},
{
"cell_type": "markdown",
"id": "outdoor-police",
"metadata": {},
"source": [
"We start by instantiating the PS, `Machine(...)`:\n",
"\n",
"(<i>hint: hit both `Shift+Tab` keys inside the parentheses `()` below to get info about the possible arguments to `Machine` and the initial values</i>)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "developed-canadian",
"metadata": {},
"outputs": [],
"source": [
"m = Machine()\n",
"\n",
"assert m.phi_s > 0, \"machine is not accelerating...?\""
]
},
{
"cell_type": "markdown",
"id": "accepted-table",
"metadata": {},
"source": [
"You can check the rf systems setup by the following plot command from the previous lecture:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "diagnostic-gibson",
"metadata": {},
"outputs": [],
"source": [
"plot_rf_overview(m);"
]
},
{
"cell_type": "markdown",
"id": "configured-terminology",
"metadata": {},
"source": [
"We initialise a Gaussian bunch distribution with very small rms bunch length $\\sigma_z=1$ m (so that the small-amplitude approximation holds):"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "figured-evening",
"metadata": {},
"outputs": [],
"source": [
"sigma_z=1"
]
},
{
"cell_type": "markdown",
"id": "reduced-penalty",
"metadata": {},
"source": [
"The \"matched\" rms momentum spread $\\sigma_{\\Delta p}$ (remember, $\\sigma_{\\Delta p}$ and $\\sigma_z$ are linked via equal Hamiltonian values $\\mathcal{H}_0$, the equilibrium condition):"
"<h3>Generating Macro-particles via Box-Muller</h3>"
]
},
{
"cell_type": "markdown",
"id": "realistic-marsh",
"metadata": {},
"source": [
"Limit by machine precision: the smallest number at FP64 is $\\varepsilon\\approx2^{-53}$, therefore one can only generate pseudo-random numbers from the Gaussian distribution up to an amplitude of $x$ where $\\exp(-x^2/2)=2^{-53}$, i.e."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "reduced-aside",
"metadata": {},
"outputs": [],
"source": [
"np.sqrt(2*-np.log(2**-53))"
]
},
{
"cell_type": "markdown",
"id": "homeless-universe",
"metadata": {},
"source": [
"$\\implies$ given $\\sigma_z=1\\,$</i>m<i>, no particles can be generated outside of $z=8.6\\,$m and the equivalent Hamiltonian contour in phase space)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "brief-presence",
"metadata": {},
"outputs": [],
"source": [
"N = 1000 # Number of macro-particles"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "classified-dietary",
"metadata": {},
"outputs": [],
"source": [
"np.random.seed(12345) # set seed to get always the same \"random\" distribution\n",
"\n",
"z = np.random.normal(loc=0, scale=sigma_z, size=N) # generate N particle positions with sigma_z around z=0\n",
"deltap = np.random.normal(loc=0, scale=sigma_deltap, size=N) # generate N particle momenta with sigma_deltap around deltap=0"
"$\\implies$ Something went wrong, since emittance has increased significantly after a few thousand turns!\n",
"\n",
"You can use the following phase-space plots as diagnostics. For this, you can stop the tracking after a certain turn to investigate, e.g. by changing the value of `n_turns` inside the `trange` counter in the `for` loop."
"$\\implies$ Emittance conserved up to about 2% if phase jump at transition is included. An additional jump in transition energy can be introduced to further reduce the emittance growth."
]
},
{
"cell_type": "markdown",
"id": "ae20d4a8-a9b7-44a3-9425-c6b2d893c141",
"metadata": {},
"source": [
"<h2>PyHEADTAIL</h2>\n",
"<h3>RF Buckets in PyHEADTAIL</h3>\n",
"\n",
"`PyHEADTAIL` provides a class to represent rf buckets (for plotting as well as for matching)"
"Calling the `RFBucketMatcher.generate` method will \n",
"\n",
"(1.) iterate on $\\mathcal{H}_0$ until the numerical integration of $\\psi(\\mathcal{H})$ for the bunch length converges to $\\hat{\\sigma}_z$, and then \n",
"\n",
"(2.) sample this $\\psi(\\mathcal{H})$ by <b>rejection sampling</b> (see previous lecture) to generate the macro-particle phase-space coordinates $(z,\\delta)$"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0446ef82-1470-437b-9abc-31d17cf3ec64",
"metadata": {},
"outputs": [],
"source": [
"z, delta, _, _ = rfb_matcher.generate(int(1e5))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b27ac4a3-d224-4e1a-9c86-73946a9e583a",
"metadata": {},
"outputs": [],
"source": [
"print(\"Converged value of H0: \" + str(rfb_matcher.psi_object.H0))\n",
"(The `scipy` implementation of the secant method in `newton()` [computes](https://github.com/scipy/scipy/blob/v1.9.3/scipy/optimize/_zeros_py.py#L330) the second guess $x_1$ as $x_1=x_0\\cdot 1.0001\\pm 0.0001$.)"
]
},
{
"cell_type": "markdown",
"id": "b35e9444-6f32-453f-90d3-306636bd7f87",
"metadata": {},
"source": [
"<h2>Numerical emittance growth</h2>"
]
},
{
"cell_type": "markdown",
"id": "aba93342-c608-47d7-b263-b3e4f97c2ec8",
"metadata": {},
"source": [
"Start with the bi-Gaussian (simulation as previous lecture):"
"$\\leadsto$ the Gaussian particle distribution is <b>not exactly</b> in equilibrium for sufficiently large rms values in the nonlinear potential, the particles <b>filament</b> and the rms emittance grows (a little)! \n",
"\n",
"$\\implies$ compare to using full nonlinear Hamiltonian to construct PDF $\\psi(\\mathcal{H})\\propto\\exp\\left(\\cfrac{\\mathcal{H}}{\\mathcal{H}_0}\\right)$"
"epsn_z = np.array([emittance(z_i, deltap_i) for z_i, deltap_i in zip(z, deltap)])\n",
"\n",
"ylim_m = np.median(4 * np.pi * epsn_z / e)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "943c4c16-c27c-43ef-90ba-d76b5a96d4d4",
"metadata": {},
"outputs": [],
"source": [
"plt.plot(4 * np.pi * epsn_z / e)\n",
"\n",
"plt.ylim(ylim_m - ylim_d, ylim_m + ylim_d)\n",
"\n",
"plt.xlabel('Turns')\n",
"plt.ylabel('$4\\pi\\epsilon_z$ [eV.s]');"
]
},
{
"cell_type": "markdown",
"id": "f04d5b2a-564b-4a15-a113-486f46748eae",
"metadata": {},
"source": [
"$\\implies$ this result shows that the nonlinearly matched thermal distribution is in equilibrium from the start (up to macro-particle noise, the fluctuations reduce with $1/\\sqrt{N}$)!"
"$\\implies$ note the offset towards negative $z$, the contours of the macro-particle density are no longer matched to the Hamiltonian contours."
]
},
{
"cell_type": "markdown",
"id": "c6800012-8075-4472-b721-c47a7a1de54f",
"metadata": {},
"source": [
"The safety `margin` inside the separatrix (where no particles are generated in `generate_gaussian_in_rfbucket`) should be chosen large enough such that no particles are located outside the rf bucket after the mismatch:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1d86c381-37a4-404c-aa63-987be0aa19a2",
"metadata": {},
"outputs": [],
"source": [
"assert all(hamiltonian(z_ini, deltap_ini, m) < 0), 'particles have been generated outside the rf bucket!'"
]
},
{
"cell_type": "markdown",
"id": "c41a169f-41ae-4e3c-b2e5-6aa2d8e30674",
"metadata": {},
"source": [
"Tracking the mismatched distribution of macro-particles:"
"$\\implies$ only residual centroid fluctuations (due to macro-particle noise), note the amplitude of the oscillation in comparison to the rf bucket length!"
]
},
{
"cell_type": "markdown",
"id": "c81c20cb-5118-4129-8b73-11b256e31ada",
"metadata": {},
"source": [
"<h3>RMS bunch length results</h3>"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e1c65718-496f-4882-8cb5-5cc6bea08bd7",
"metadata": {},
"outputs": [],
"source": [
"plt.plot(np.std(z, axis=1))\n",
"\n",
"plt.xlabel('Turns')\n",
"plt.ylabel(r'$\\sigma_z$ [m]');"
]
},
{
"cell_type": "markdown",
"id": "069a9bd6-bd05-4554-a7e7-6e7c9eeb696d",
"metadata": {},
"source": [
"$\\implies$ exponential decay of the initial momentum mismatch (due to the non-linearity of the rf bucket)"
]
},
{
"cell_type": "markdown",
"id": "d3204663-1fbe-49d9-af46-eef592089372",
"metadata": {},
"source": [
"<h3>RMS emittance results</h3>"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1bae04d7-57e7-47d5-bd16-9960bf9d8141",
"metadata": {},
"outputs": [],
"source": [
"epsn_z = np.array([emittance(z_i, deltap_i) for z_i, deltap_i in zip(z, deltap)])"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9bc47eca-c715-4352-9e9e-6f18112c9d57",
"metadata": {},
"outputs": [],
"source": [
"plt.plot(epsn_z / e)\n",
"\n",
"plt.xlabel('Turns')\n",
"plt.ylabel('$\\epsilon_z$ [eV.s]');"
]
},
{
"cell_type": "markdown",
"id": "bc13e46e-0d5f-4a08-bd2c-adff48bf2380",
"metadata": {},
"source": [
"$\\implies$ in this example, 20% emittance growth as a result of the 50% momentum spread mismatch."
<h2>Simulation of Full CERN PS Acceleration Ramp</h2>
<h3>CERN PS Machine parameters</h3>
%% Cell type:markdown id:taken-reminder tags:
The CERN Proton Synchrotron
- has a circumference of 2π·100m
- takes protons from the PS Booster at a kinetic energy of 2GeV corresponding to a $\gamma=3.13$
- injects with 50kV of rf voltage, up to 200kV for ramp
- runs at harmonic $h=7$
- has a momentum compaction factor of $\alpha_c=0.027$
- typical acceleration rate of (up to) $\dot{B}=2$ T/s, the bending radius is $\rho=70.08$ m
%% Cell type:markdown id:outdoor-police tags:
We start by instantiating the PS, `Machine(...)`:
(<i>hint: hit both `Shift+Tab` keys inside the parentheses `()` below to get info about the possible arguments to `Machine` and the initial values</i>)
%% Cell type:code id:developed-canadian tags:
``` python
m=Machine()
assertm.phi_s>0,"machine is not accelerating...?"
```
%% Cell type:markdown id:accepted-table tags:
You can check the rf systems setup by the following plot command from the previous lecture:
We initialise a Gaussian bunch distribution with very small rms bunch length $\sigma_z=1$ m (so that the small-amplitude approximation holds):
%% Cell type:code id:figured-evening tags:
``` python
sigma_z=1
```
%% Cell type:markdown id:reduced-penalty tags:
The "matched" rms momentum spread $\sigma_{\Delta p}$ (remember, $\sigma_{\Delta p}$ and $\sigma_z$ are linked via equal Hamiltonian values $\mathcal{H}_0$, the equilibrium condition):
<h3>Generating Macro-particles via Box-Muller</h3>
%% Cell type:markdown id:realistic-marsh tags:
Limit by machine precision: the smallest number at FP64 is $\varepsilon\approx2^{-53}$, therefore one can only generate pseudo-random numbers from the Gaussian distribution up to an amplitude of $x$ where $\exp(-x^2/2)=2^{-53}$, i.e.
%% Cell type:code id:reduced-aside tags:
``` python
np.sqrt(2*-np.log(2**-53))
```
%% Cell type:markdown id:homeless-universe tags:
$\implies$ given $\sigma_z=1\,$</i>m<i>, no particles can be generated outside of $z=8.6\,$m and the equivalent Hamiltonian contour in phase space)
%% Cell type:code id:brief-presence tags:
``` python
N=1000# Number of macro-particles
```
%% Cell type:code id:classified-dietary tags:
``` python
np.random.seed(12345)# set seed to get always the same "random" distribution
z=np.random.normal(loc=0,scale=sigma_z,size=N)# generate N particle positions with sigma_z around z=0
deltap=np.random.normal(loc=0,scale=sigma_deltap,size=N)# generate N particle momenta with sigma_deltap around deltap=0
We have reached the extraction energy of $\gamma=27.7$:
%% Cell type:code id:front-brunei tags:
``` python
m.gamma_ref
```
%% Cell type:markdown id:given-absence tags:
<h3>How did we do?</h3>
Check the rms emittance:
%% Cell type:code id:sharp-rescue tags:
``` python
plt.plot(np.arange(n_turns)/100000,epsn_z/e)
plt.xlabel('Turns [100k]')
plt.ylabel('$\epsilon_z$ [eV.s]');
```
%% Cell type:markdown id:typical-collection tags:
$\implies$ Something went wrong, since emittance has increased significantly after a few thousand turns!
You can use the following phase-space plots as diagnostics. For this, you can stop the tracking after a certain turn to investigate, e.g. by changing the value of `n_turns` inside the `trange` counter in the `for` loop.
%% Cell type:code id:furnished-links tags:
``` python
plot_hamiltonian(m);
plt.scatter(z%600,deltap/m.p0(),marker='.',s=1)
```
%% Cell type:code id:electrical-curve tags:
``` python
np.random.seed(12345)# set seed to get always the same "random" distribution
z=np.random.normal(loc=0,scale=sigma_z,size=N)# generate N particle positions with sigma_z around z=0
deltap=np.random.normal(loc=0,scale=sigma_deltap,size=N)# generate N particle momenta with sigma_deltap around deltap=0
m=Machine()
```
%% Cell type:code id:related-privacy tags:
``` python
fori_turnintrange(1,int(n_turns/100)):
z,deltap=track_one_turn(z,deltap,m)
epsn_z[i_turn]=emittance(z,deltap)
```
%% Cell type:code id:boring-comedy tags:
``` python
print(m.gamma_ref)
plot_hamiltonian(m);
plt.scatter(z,deltap/m.p0(),marker='.',s=1)
```
%% Cell type:code id:faced-princess tags:
``` python
m.gamma_ref
```
%% Cell type:code id:constant-newark tags:
``` python
np.sqrt(1/m.alpha_c)
```
%% Cell type:markdown id:amber-retreat tags:
<h3>Crossing transition</h3>
%% Cell type:code id:limited-reset tags:
``` python
np.random.seed(12345)# set seed to get always the same "random" distribution
z=np.random.normal(loc=0,scale=sigma_z,size=N)# generate N particle positions with sigma_z around z=0
deltap=np.random.normal(loc=0,scale=sigma_deltap,size=N)# generate N particle momenta with sigma_deltap around deltap=0
m=Machine()
```
%% Cell type:code id:visible-school tags:
``` python
phi_s_1=np.pi-m.phi_s# synchronous phase after transition calculated from synchronous phase before transition
fori_turnintrange(1,n_turns):
z,deltap=track_one_turn(z,deltap,m)
epsn_z[i_turn]=emittance(z,deltap)
condition=m.gamma_ref>np.sqrt(1/m.alpha_c)# condition to change synchronous phase: gamma above gamma_transition
ifcondition:
m.phi_s=phi_s_1
```
%% Cell type:code id:capable-ladder tags:
``` python
plt.plot(np.arange(n_turns)/100000,epsn_z/e)
plt.xlabel('Turns [100k]')
plt.ylabel('$\epsilon_z$ [eV.s]');
```
%% Cell type:markdown id:velvet-technician tags:
$\implies$ Emittance conserved up to about 2% if phase jump at transition is included. An additional jump in transition energy can be introduced to further reduce the emittance growth.
Calling the `RFBucketMatcher.generate` method will
(1.) iterate on $\mathcal{H}_0$ until the numerical integration of $\psi(\mathcal{H})$ for the bunch length converges to $\hat{\sigma}_z$, and then
(2.) sample this $\psi(\mathcal{H})$ by <b>rejection sampling</b> (see previous lecture) to generate the macro-particle phase-space coordinates $(z,\delta)$
(The `scipy` implementation of the secant method in `newton()` [computes](https://github.com/scipy/scipy/blob/v1.9.3/scipy/optimize/_zeros_py.py#L330) the second guess $x_1$ as $x_1=x_0\cdot 1.0001\pm 0.0001$.)
$\leadsto$ the Gaussian particle distribution is <b>not exactly</b> in equilibrium for sufficiently large rms values in the nonlinear potential, the particles <b>filament</b> and the rms emittance grows (a little)!
$\implies$ compare to using full nonlinear Hamiltonian to construct PDF $\psi(\mathcal{H})\propto\exp\left(\cfrac{\mathcal{H}}{\mathcal{H}_0}\right)$
$\implies$ this result shows that the nonlinearly matched thermal distribution is in equilibrium from the start (up to macro-particle noise, the fluctuations reduce with $1/\sqrt{N}$)!
The safety `margin` inside the separatrix (where no particles are generated in `generate_gaussian_in_rfbucket`) should be chosen large enough such that no particles are located outside the rf bucket after the mismatch:
$\implies$ only residual centroid fluctuations (due to macro-particle noise), note the amplitude of the oscillation in comparison to the rf bucket length!
