diff --git a/week5/nmap5.ipynb b/week5/nmap5.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..a251577116f8c163b6a55c03322f98eb1e4e23b4
--- /dev/null
+++ b/week5/nmap5.ipynb
@@ -0,0 +1,2013 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "incomplete-medline",
+ "metadata": {},
+ "source": [
+ "<h1 style=\"text-align: center; vertical-align: middle;\">Numerical Methods in Accelerator Physics</h1>\n",
+ "<h2 style=\"text-align: center; vertical-align: middle;\">Python examples -- Week 5</h2>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "executive-television",
+ "metadata": {},
+ "source": [
+ "<h2>Run this first!</h2>\n",
+ "\n",
+ "Imports and modules:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "biological-product",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from config5 import np, plt, plot_rfwave, tqdm, trange\n",
+ "from scipy.constants import m_p, e, c\n",
+ "import PyNAFF\n",
+ "%matplotlib inline"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "transsexual-renewal",
+ "metadata": {},
+ "source": [
+ "If the progress bar by `tqdm` (`trange`) later in this document does not work, run this:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "maritime-warner",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "!jupyter nbextension enable --py widgetsnbextension"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "present-mother",
+ "metadata": {},
+ "source": [
+ "<h3>Machine parameters for exercises</h3>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "taken-reminder",
+ "metadata": {},
+ "source": [
+ "The CERN Proton Synchrotron (PS):\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": "primary-sullivan",
+ "metadata": {},
+ "source": [
+ "<h3 style=\"color: #e6541a;\">Exercises for topics of week 4 (Solutions below)</h3>\n",
+ "\n",
+ "<div style=\"color: #e6541a; margin-top: 2em;\">Based on the CERN Proton Synchrotron (PS) machine parameters: <br />\n",
+ "\n",
+ "1. Track particles given a stationary synchronous particle ($\\varphi_s=0$).\n",
+ "2. Determine the frequency of the synchrotron motion (via NAFF).\n",
+ "3. Track particles given an accelerating synchronous particle ($\\varphi_s>0$).\n",
+ "4. Compute the transition energy $\\gamma_{\\text{t}}$ for CERN PS.\n",
+ "4. Set machine energy ($\\gamma$) above transition / relativistic regime, track like (4) again.\n",
+ "\n",
+ "$\\implies$ take note of the observed trajectories in phase space (for the following)!\n",
+ "</div>\n",
+ "<!-- \n",
+ "3. Derivation of Hamiltonian $\\mathcal{H}(z,\\delta)=T(\\delta) + V(z)$ for stationary synchronous particle ($\\varphi_s=0$) -->"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "sensitive-nature",
+ "metadata": {},
+ "source": [
+ "Some convenience functions to compute the speed β and the relativistic Lorentz factor γ:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "fifteen-utility",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def beta(gamma):\n",
+ " '''Speed β in units of c from relativistic Lorentz factor γ.'''\n",
+ " return np.sqrt(1 - gamma**-2)\n",
+ "\n",
+ "def gamma(p):\n",
+ " '''Relativistic Lorentz factor γ from total momentum p.'''\n",
+ " return np.sqrt(1 + (p / (mass * c))**2)\n",
+ "\n",
+ "def p_from_gamma(gamma):\n",
+ " '''Total momentum p from relativistic Lorentz factor γ.'''\n",
+ " return mass * c * np.sqrt(gamma**2 - 1)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "outdoor-police",
+ "metadata": {},
+ "source": [
+ "We gather all the machine parameters in a class named `Machine`.\n",
+ "\n",
+ "A `Machine` instance knows\n",
+ "- at which energy the synchronous particle (reference γ, `gamma_ref`, or alternatively the momentum `p0()`) currently runs,\n",
+ "- what the acceleration rate is in terms of the synchronous phase φ_s, `phi_s`\n",
+ "- how to compute the phase-slip factor $\\eta$, `eta`, for a particle at a certain momentum p0 + Δp\n",
+ "- how to update the energy of the synchronous particle via `update_gamma_ref()` when a turn has passed"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "developed-canadian",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "charge = e\n",
+ "mass = m_p\n",
+ "\n",
+ "class Machine(object):\n",
+ " gamma_ref = 3.13\n",
+ " circumference = 2*np.pi*100\n",
+ " voltage = 50e3\n",
+ " harmonic = 7\n",
+ " alpha_c = 0.027\n",
+ " phi_s = 0\n",
+ " \n",
+ " def __init__(self, gamma_ref=gamma_ref, circumference=circumference,\n",
+ " voltage=voltage, harmonic=harmonic, \n",
+ " alpha_c=alpha_c, phi_s=phi_s):\n",
+ " '''Override default settings by giving explicit arguments.'''\n",
+ " self.gamma_ref = gamma_ref\n",
+ " self.circumference = circumference\n",
+ " self.voltage = voltage\n",
+ " self.harmonic = harmonic\n",
+ " self.alpha_c = alpha_c\n",
+ " self.phi_s = phi_s\n",
+ " \n",
+ " def eta(self, deltap):\n",
+ " '''Phase-slip factor for a particle.'''\n",
+ " p = self.p0() + deltap\n",
+ " return self.alpha_c - gamma(p)**-2\n",
+ "\n",
+ " def p0(self):\n",
+ " '''Momentum of synchronous particle.'''\n",
+ " return self.gamma_ref * beta(self.gamma_ref) * mass * c\n",
+ "\n",
+ " def update_gamma_ref(self):\n",
+ " '''Advance the energy of the synchronous particle\n",
+ " according to the synchronous phase by one turn.\n",
+ " '''\n",
+ " deltap_per_turn = charge * self.voltage / (\n",
+ " beta(self.gamma_ref) * c) * np.sin(self.phi_s)\n",
+ " new_p0 = self.p0() + deltap_per_turn\n",
+ " self.gamma_ref = gamma(new_p0)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "accepted-table",
+ "metadata": {},
+ "source": [
+ "To compute a synchronous phase, you may use these convenience functions to compute the $\\arcsin$ on the interval of [-π/2,π/2] .\n",
+ "\n",
+ "Remember you may likely want to find a synchronous phase on $\\varphi_s\\in[0,π]$ !"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "known-pitch",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def deltap_per_turn(Bdot, rho, circumference, p0):\n",
+ " Trev = circumference / (beta(gamma(p0)) * c)\n",
+ " return Bdot * rho * charge * Trev"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "micro-alpha",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def compute_phi_s(deltap_per_turn, p0, voltage):\n",
+ " '''Return *first* positive phase which matches the\n",
+ " given Δp/turn on the interval [0, π/2].\n",
+ " Do check whether you need to use π-φ_s for stability!\n",
+ " '''\n",
+ " return np.arcsin(\n",
+ " deltap_per_turn * beta(gamma(p0)) * c / (charge * voltage)\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "literary-absorption",
+ "metadata": {},
+ "source": [
+ "The tracking equations for the longitudinal plane read:\n",
+ "\n",
+ "$$\\begin{cases}\\,\n",
+ " z_{n+1} &= z_n - \\eta C \\left(\\cfrac{\\Delta p}{p_0}\\right)_n \\\\\n",
+ " (\\Delta p)_{n+1} &= (\\Delta p)_n + \\cfrac{q V}{(\\beta c)_n}\\cdot\\left(\\sin\\left(\\varphi_s - \\cfrac{2\\pi}{C}\\cdot hz_{n+1}\\right) - \\sin(\\varphi_s)\\right)\n",
+ "\\end{cases}$$\n",
+ "\n",
+ "We implement them with a leapfrog (half-drift + kick + half-drift) scheme."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "appointed-complex",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def track_one_turn(z_n, deltap_n, machine):\n",
+ " m = machine\n",
+ " # half drift\n",
+ " z_nhalf = z_n - m.eta(deltap_n) * deltap_n / m.p0() * m.circumference / 2\n",
+ " # rf kick\n",
+ " amplitude = charge * m.voltage / (beta(gamma(m.p0())) * c)\n",
+ " phi = m.phi_s - m.harmonic * 2 * np.pi * z_nhalf / m.circumference\n",
+ " \n",
+ " m.update_gamma_ref()\n",
+ " deltap_n1 = deltap_n + amplitude * (np.sin(phi) - np.sin(m.phi_s))\n",
+ " # half drift\n",
+ " z_n1 = z_nhalf - m.eta(deltap_n1) * deltap_n1 / m.p0() * m.circumference / 2\n",
+ " return z_n1, deltap_n1"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "included-meaning",
+ "metadata": {},
+ "source": [
+ "<b>1. Track particles given a stationary synchronous particle</b><br />\n",
+ "The `Machine` instance will keep track of the reference energy during the tracking by calling `update_gamma_ref()` once per turn:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "familiar-dividend",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m = Machine()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "welcome-latex",
+ "metadata": {},
+ "source": [
+ "Particles are tracked by their two longitudinal coordinates $(z, \\Delta p)$. The initial values are stored in `z_ini` and `deltap_ini` as `numpy.array`s. These should have `N` entries for $N$ particles.\n",
+ "\n",
+ "(You may use numpy helper functions such as `np.linspace` or `np.arange` for convenient initialisation!)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "utility-classic",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_turns = 3000\n",
+ "deltap_ini = np.linspace(0, 0.01 * m.p0(), 40)#np.array([0.]) #np.linspace(start, end, 20)\n",
+ "z_ini = np.zeros_like(deltap_ini)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "purple-implementation",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = len(z_ini)\n",
+ "assert (N == len(deltap_ini))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "formed-james",
+ "metadata": {},
+ "source": [
+ "To store the coordinate values during tracking, prepare some `n_turns` long 2D arrays with `N` entries per turn:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "included-negotiation",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "z = np.zeros((n_turns, N), dtype=np.float64)\n",
+ "deltap = np.zeros_like(z)\n",
+ "\n",
+ "z[0] = z_ini\n",
+ "deltap[0] = deltap_ini"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "under-blanket",
+ "metadata": {},
+ "source": [
+ "We would also like to store the reference gamma for each turn:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "played-marker",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "gammas = np.zeros(n_turns, dtype=np.float64)\n",
+ "gammas[0] = m.gamma_ref"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "atomic-uzbekistan",
+ "metadata": {},
+ "source": [
+ "Let's go, here's the tracking loop over the number of turns `n_turns`!"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "broken-inside",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i_turn in range(1, n_turns):\n",
+ " z[i_turn], deltap[i_turn] = track_one_turn(z[i_turn - 1], deltap[i_turn - 1], m)\n",
+ " gammas[i_turn] = m.gamma_ref"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "opening-stick",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(gammas)\n",
+ "plt.xlabel('Turns')\n",
+ "plt.ylabel('$\\gamma_{ref}$')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "authorized-robinson",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "plt.scatter(z, deltap / m.p0(), marker='.', s=0.5)\n",
+ "plt.xlim(-50,100)\n",
+ "plt.xlabel('$z$ [m]')\n",
+ "plt.ylabel('$\\Delta p/p_0$')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "altered-morocco",
+ "metadata": {},
+ "source": [
+ "<b>2. Determine the frequency of the synchrotron motion (via NAFF).</b>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "labeled-clear",
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "plt.plot(((z.T)[:19]).T)\n",
+ "plt.xlabel('turns')\n",
+ "plt.ylabel('$z$ [m]')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "sonic-anthropology",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "freqs_naff = []\n",
+ "\n",
+ "for signal in ((z.T)[:19]):\n",
+ " freq_naff = PyNAFF.naff(signal, turns=n_turns - 1, nterms=1)\n",
+ " try:\n",
+ " freq_naff = freq_naff[0, 1]\n",
+ " except IndexError:\n",
+ " freq_naff = 0\n",
+ " freqs_naff += [freq_naff]\n",
+ "\n",
+ "freqs_naff = np.array(freqs_naff)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "surgical-newport",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "plt.plot(np.max(((z.T)[:19]).T,axis=0), freqs_naff, c='r', marker='.')\n",
+ "\n",
+ "#plt.xticks([-np.pi/2, 0, np.pi/2, np.pi], [r\"$-\\pi/2$\", \"0\", r\"$\\pi/2$\", r\"$\\pi$\"])\n",
+ "plt.xlabel(r'Maximum amplitude z [m]')\n",
+ "plt.ylabel('Phase advance / \\n' + r'integration step $\\Delta t$ [$2\\pi$]');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "rotary-accordance",
+ "metadata": {},
+ "source": [
+ "<b>3. Track particles given an accelerating synchronous particle ($\\gamma=3.13$, typical acceleration rate of $\\dot{B}=2$ T/s, bending radius $\\rho=70.08$ m, rf voltage 200kV)</b>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "novel-collective",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "deltap_per_turn(2,70.08,2*np.pi*100,p_from_gamma(3.13))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "initial-cinema",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "phi_s=compute_phi_s(deltap_per_turn(2,70.08,2*np.pi*100,p_from_gamma(3.13)),p_from_gamma(3.13),200e3)\n",
+ "phi_s"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "streaming-pleasure",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m = Machine(phi_s=phi_s)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "working-dictionary",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_turns = 1000\n",
+ "deltap_ini = np.linspace(0, 0.01 * m.p0(), 20)#np.array([0.]) #np.linspace(start, end, 20)\n",
+ "z_ini = np.zeros_like(deltap_ini)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "anticipated-china",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = len(z_ini)\n",
+ "assert (N == len(deltap_ini))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "tested-jurisdiction",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "z = np.zeros((n_turns, N), dtype=np.float64)\n",
+ "deltap = np.zeros_like(z)\n",
+ "\n",
+ "z[0] = z_ini\n",
+ "deltap[0] = deltap_ini"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "accepting-northwest",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "gammas = np.zeros(n_turns, dtype=np.float64)\n",
+ "gammas[0] = m.gamma_ref"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "coordinated-evans",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i_turn in range(1, n_turns):\n",
+ " z[i_turn], deltap[i_turn] = track_one_turn(z[i_turn - 1], deltap[i_turn - 1], m)\n",
+ " gammas[i_turn] = m.gamma_ref"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "genetic-night",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(gammas)\n",
+ "plt.xlabel('Turns')\n",
+ "plt.ylabel('$\\gamma_{ref}$')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "tribal-minimum",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "plt.scatter(z, deltap / m.p0(), marker='.', s=0.5)\n",
+ "plt.xlabel('$z$ [m]')\n",
+ "plt.ylabel('$\\Delta p/p_0$')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "changed-check",
+ "metadata": {},
+ "source": [
+ "<b>4. Compute the transition energy for CERN PS</b>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "meaningful-quantum",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "gamma_t=1/np.sqrt(m.alpha_c)\n",
+ "gamma_t"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "based-piano",
+ "metadata": {},
+ "source": [
+ "<b>5. Set machine energy above transition / relativistic regime, track again.</b>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "raised-breach",
+ "metadata": {},
+ "source": [
+ "<b>Same synchronous phase => longitudinal defocusing</b>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "described-clarity",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m = Machine(gamma_ref=7,phi_s=phi_s)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "neutral-corruption",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_turns = 5000\n",
+ "deltap_ini = np.linspace(0, 0.01 * m.p0(), 20)#np.array([0.]) #np.linspace(start, end, 20)\n",
+ "z_ini = np.zeros_like(deltap_ini)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "conditional-moldova",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = len(z_ini)\n",
+ "assert (N == len(deltap_ini))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "skilled-concord",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "z = np.zeros((n_turns, N), dtype=np.float64)\n",
+ "deltap = np.zeros_like(z)\n",
+ "\n",
+ "z[0] = z_ini\n",
+ "deltap[0] = deltap_ini"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "injured-worker",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "gammas = np.zeros(n_turns, dtype=np.float64)\n",
+ "gammas[0] = m.gamma_ref"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "european-transportation",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i_turn in range(1, n_turns):\n",
+ " z[i_turn], deltap[i_turn] = track_one_turn(z[i_turn - 1], deltap[i_turn - 1], m)\n",
+ " gammas[i_turn] = m.gamma_ref"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "preceding-provincial",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(gammas)\n",
+ "plt.xlabel('Turns')\n",
+ "plt.ylabel('$\\gamma_{ref}$')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "sharp-denmark",
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "plt.scatter(z, deltap / m.p0(), marker='.', s=0.5)\n",
+ "plt.xlabel('$z$ [m]')\n",
+ "plt.ylabel('$\\Delta p/p_0$')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "confirmed-experiment",
+ "metadata": {},
+ "source": [
+ "<b>Same synchronous phase => longitudinal defocusing</b>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "regulated-delta",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m = Machine(gamma_ref=7,phi_s=np.pi-phi_s)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "exciting-picking",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_turns = 10000\n",
+ "deltap_ini = np.linspace(0, 0.01 * m.p0(), 20)#np.array([0.]) #np.linspace(start, end, 20)\n",
+ "z_ini = np.zeros_like(deltap_ini)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "under-engineering",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = len(z_ini)\n",
+ "assert (N == len(deltap_ini))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "bibliographic-tiffany",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "z = np.zeros((n_turns, N), dtype=np.float64)\n",
+ "deltap = np.zeros_like(z)\n",
+ "\n",
+ "z[0] = z_ini\n",
+ "deltap[0] = deltap_ini"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "tracked-hierarchy",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "gammas = np.zeros(n_turns, dtype=np.float64)\n",
+ "gammas[0] = m.gamma_ref"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "becoming-paraguay",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i_turn in range(1, n_turns):\n",
+ " z[i_turn], deltap[i_turn] = track_one_turn(z[i_turn - 1], deltap[i_turn - 1], m)\n",
+ " gammas[i_turn] = m.gamma_ref"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "distant-edgar",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(gammas)\n",
+ "plt.xlabel('Turns')\n",
+ "plt.ylabel('$\\gamma_{ref}$')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "final-hebrew",
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "plt.scatter(z, deltap / m.p0(), marker='.', s=0.5)\n",
+ "plt.xlim(right=100)\n",
+ "plt.ylim(bottom=-0.02)\n",
+ "plt.xlabel('$z$ [m]')\n",
+ "plt.ylabel('$\\Delta p/p_0$')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "coated-supplement",
+ "metadata": {},
+ "source": [
+ "<h3>Week 5: Visualization of synchrotron tune by particle tracking (slide 16)</h3>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "equipped-recording",
+ "metadata": {},
+ "source": [
+ "You can override the default PS parameters by supplying arguments to the `Machine(...)` instantiation, e.g. change the `gamma_ref` by `m = Machine(gamma_ref=20)` or the `phi_s` etc."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "corrected-repair",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m = Machine(voltage=200e3,phi_s=0.456)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "foreign-louisiana",
+ "metadata": {},
+ "source": [
+ "Now we define the kinetic and potential energy terms and the Hamiltonian function(al) itself:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "steady-arbor",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def T(deltap, machine):\n",
+ " '''Kinetic energy term in Hamiltonian.'''\n",
+ " return -0.5 * machine.eta(deltap) / machine.p0() * deltap**2\n",
+ "\n",
+ "def U(z, machine, beta_=None):\n",
+ " '''Potential energy term in Hamiltonian.\n",
+ " If beta is not given, compute it from synchronous particle.\n",
+ " '''\n",
+ " m = machine\n",
+ " if beta_ is None:\n",
+ " beta_ = beta(gamma(m.p0()))\n",
+ " ampl = charge * m.voltage / (beta_ * c * 2 * np.pi * m.harmonic)\n",
+ " phi = m.phi_s - 2 * np.pi * m.harmonic / m.circumference * z\n",
+ " # convenience: define z at unstable fixed point\n",
+ " z_ufp = -m.circumference * (np.pi - 2 * m.phi_s) / (2 * np.pi * m.harmonic)\n",
+ " # convenience: offset by potential value at unstable fixed point\n",
+ " # such that unstable fixed point (and separatrix) have 0 potential energy\n",
+ " return ampl * (-np.cos(phi) + \n",
+ " 2 * np.pi * m.harmonic / m.circumference * (z - z_ufp) * np.sin(m.phi_s) +\n",
+ " -np.cos(m.phi_s))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "flush-defense",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def hamiltonian(z, deltap, machine):\n",
+ " return T(deltap, machine) + U(z, machine, beta_=beta(gamma(machine.p0() + deltap)))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "foreign-swedish",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def plot_hamiltonian(machine, zleft=-50, zright=50, dpmax=0.01, cbar=True):\n",
+ " '''Plot Hamiltonian contours across (zleft, zright) and (-dpmax, dpmax).'''\n",
+ " Z, DP = np.meshgrid(np.linspace(zleft, zright, num=1000), \n",
+ " np.linspace(-dpmax, dpmax, num=1000))\n",
+ " H = hamiltonian(Z, DP * machine.p0(), machine) / machine.p0()\n",
+ " \n",
+ " plt.contourf(Z, DP, H, cmap=plt.get_cmap('hot_r'), levels=12,\n",
+ " zorder=0, alpha=0.5)\n",
+ " plt.xlabel('$z$ [m]')\n",
+ " plt.ylabel(r'$\\delta$')\n",
+ " if cbar:\n",
+ " colorbar = plt.colorbar(label=r'$\\mathcal{H}(z,\\Delta p)\\,/\\,p_0$')\n",
+ " colorbar.ax.axhline(0, lw=2, c='b')\n",
+ " plt.contour(Z, DP, H, colors='b', linewidths=2, levels=[0])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "planned-values",
+ "metadata": {},
+ "source": [
+ "Let's plot the Hamiltonian landscape on phase space:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "stone-sixth",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plot_hamiltonian(m)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "after-legislature",
+ "metadata": {},
+ "source": [
+ "Tracking just like last lecture:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "confused-sunrise",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_turns = 500\n",
+ "deltap_ini = np.linspace(0, 0.01 * m.p0(), 20)\n",
+ "z_ini = np.zeros_like(deltap_ini)\n",
+ "\n",
+ "N = len(z_ini)\n",
+ "assert (N == len(deltap_ini))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "global-marsh",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "z = np.zeros((n_turns, N), dtype=np.float64)\n",
+ "deltap = np.zeros_like(z)\n",
+ "\n",
+ "z[0] = z_ini\n",
+ "deltap[0] = deltap_ini"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "continent-headquarters",
+ "metadata": {},
+ "source": [
+ "Record the evolution of energy $\\gamma_\\mathrm{ref}$ and also Hamiltonian values of particles during tracking:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "purple-token",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "gammas = np.zeros(n_turns, dtype=np.float64)\n",
+ "gammas[0] = m.gamma_ref\n",
+ "\n",
+ "H_values = np.zeros_like(z)\n",
+ "H_values[0] = hamiltonian(z_ini, deltap_ini, m) / m.p0()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "unable-wages",
+ "metadata": {},
+ "source": [
+ "The tracking loop:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "finnish-south",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i_turn in range(1, n_turns):\n",
+ " z[i_turn], deltap[i_turn] = track_one_turn(z[i_turn - 1], deltap[i_turn - 1], m)\n",
+ " gammas[i_turn] = m.gamma_ref\n",
+ " H_values[i_turn] = hamiltonian(z[i_turn], deltap[i_turn], m) / m.p0()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "taken-configuration",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(gammas)\n",
+ "plt.xlabel('Turns')\n",
+ "plt.ylabel('$\\gamma_{ref}$');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "computational-specification",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.scatter(z, deltap / m.p0(), marker='.', s=0.5)\n",
+ "plt.xlabel('$z$ [m]')\n",
+ "plt.ylabel('$\\Delta p/p_0$')\n",
+ "plot_hamiltonian(m, zleft=-150, zright=50, dpmax=0.017)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "political-bench",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(H_values, c='C0');\n",
+ "plt.xlabel('Turns')\n",
+ "plt.ylabel(r'$\\mathcal{H}(z,\\Delta p)\\,/\\,p_0$');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "backed-sentence",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def plot_rf_overview():\n",
+ " z_range = np.linspace(-150, 40, num=1000)\n",
+ " # z location of unstable fixed point:\n",
+ " z_ufp = -m.circumference * (np.pi - 2 * m.phi_s) / (2 * np.pi * m.harmonic)\n",
+ "\n",
+ " fig, ax = plt.subplots(3, 1, figsize=(6, 10), sharex=True)\n",
+ "\n",
+ " plt.sca(ax[0])\n",
+ " plt.plot(z_range, 1e-3 * m.voltage * np.sin(m.phi_s - 2 * np.pi * m.harmonic / m.circumference * z_range))\n",
+ " plt.axhline(0, c='gray', lw=2)\n",
+ " plt.axhline(1e-3 * m.voltage * np.sin(m.phi_s), c='purple', lw=2, ls='--')\n",
+ " plt.axvline(0, c='purple', lw=2)\n",
+ " plt.axvline(z_ufp, c='red', lw=2)\n",
+ " plt.ylabel('rf wave $V(z)$ [kV]')\n",
+ "\n",
+ " plt.sca(ax[1])\n",
+ " plt.plot(z_range, 1e6 * U(z_range, m) / m.p0())\n",
+ " plt.axhline(0, c='gray', lw=2)\n",
+ " plt.ylabel(r'$U(z)\\,/\\,p_0\\cdot 10^6$')\n",
+ "\n",
+ " plt.scatter([z_ufp], [0], marker='*', c='white', edgecolor='red', zorder=10)\n",
+ " plt.scatter([0], [U(0, m) / m.p0()], marker='d', c='white', edgecolor='purple', zorder=10)\n",
+ "\n",
+ " plt.sca(ax[2])\n",
+ " plot_hamiltonian(m, zleft=z_range[0], zright=z_range[-1], cbar=False)\n",
+ " plt.scatter([z_ufp], [0], marker='*', c='white', edgecolor='red', zorder=10)\n",
+ " plt.scatter([0], [0], marker='d', c='white', edgecolor='purple')\n",
+ " plt.xlabel('$z$ [m]')\n",
+ " plt.ylabel('$\\delta$')\n",
+ " plt.subplots_adjust(hspace=0)\n",
+ " \n",
+ " return fig, ax"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "exposed-retrieval",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plot_rf_overview();"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "classical-proxy",
+ "metadata": {},
+ "source": [
+ "<h2 style=\"color: #e6541a;\">Exercises based on Tracking (Solutions not included)</h2>\n",
+ "\n",
+ "<div style=\"color: #e6541a; margin-top: 2em;\">\n",
+ "Modify the machine parameters and particle initial conditions to answer the following questions: <br />\n",
+ " \n",
+ "1. Do the Hamiltonian contours predict the tracked particle trajectories well?\n",
+ "---\n",
+ "2. Compare the stationary synchronous particle situation to the nonlinear pendulum: what is the meaning of the $\\pi$ phase offset of $\\varphi_s$ in terms of the pendulum? What state of the pendulum does the stable and what the unstable fixed point correspond to?\n",
+ "---\n",
+ "3. CERN PS accelerates protons from $\\gamma=3.1$ up to $\\gamma=27.7$. Is the transition energy crossed during acceleration? What does this mean for the setting of the synchronous phase $\\varphi_s$ in the control room? \n",
+ "---\n",
+ "4. What happens to phase focusing at the transition energy $\\gamma=\\gamma_\\mathrm{t}$? Can you explain why? (Think of the phase slippage mechanism.)\n",
+ "</div>\n",
+ "\n",
+ "$\\implies$ please make sure you understand and note down the answers to these questions: the concepts are central to accelerator physics and relevant for the exam."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "covered-toner",
+ "metadata": {},
+ "source": [
+ "<p style=\"color: #e6541a;\">$\\implies$ <i>Bonus: determine the synchrotron tune from tracking simulations (using NAFF) and compare to the derived formula for $Q_s$!</i></p>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "spatial-sodium",
+ "metadata": {},
+ "source": [
+ "<h3>Linear Congruential Generator (slide 19)</h3>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "armed-strain",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "class RandomNumberGenerator(object):\n",
+ " def __init__(self, M, a, c, seed):\n",
+ " self.M = M\n",
+ " self.a = a\n",
+ " self.c = c\n",
+ " self.xk = seed\n",
+ "\n",
+ " def generate(self):\n",
+ " xk1 = (self.a * self.xk + self.c) % self.M\n",
+ " self.xk = xk1\n",
+ " return xk1 / self.M"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "insured-sellers",
+ "metadata": {},
+ "source": [
+ "Instantiate the linear congruential generator by Lewis et al with a certain `seed`:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "sought-cathedral",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "prng_standard = RandomNumberGenerator(\n",
+ " M=2**31 - 1, \n",
+ " a=7**5, \n",
+ " c=0, \n",
+ " seed=12345)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "dense-comfort",
+ "metadata": {},
+ "source": [
+ "Generate a set of numbers from the sequence and analyse:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "french-modem",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "results = [prng_standard.generate() for i in range(10000)]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "continued-accordance",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.hist(results)\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel('#draws');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "heated-armenia",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "plt.scatter(results[:-1], results[1:], s=1, marker='.')\n",
+ "plt.xlabel('$x_k$')\n",
+ "plt.ylabel('$x_{k+1}$');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "packed-respondent",
+ "metadata": {},
+ "source": [
+ "<p style=\"color: #e6541a;\">$\\implies$ What happens when you change the parameters $M,a,c$?<br /><br />\n",
+ "Try e.g. $a=5$ or $M=2^{31}-2$...</p>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "decimal-photograph",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "prng_standard2 = RandomNumberGenerator(\n",
+ " M=2**31 - 2, \n",
+ " a=7**5, \n",
+ " c=0, \n",
+ " seed=12345)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "generic-array",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "results2 = [prng_standard2.generate() for i in range(10000)]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "adult-joseph",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.hist(results2)\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel('#draws');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ideal-dallas",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "plt.scatter(results2[:-1], results2[1:], s=1, marker='.')\n",
+ "plt.xlabel('$x_k$')\n",
+ "plt.ylabel('$x_{k+1}$');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "paperback-certification",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "prng_standard3 = RandomNumberGenerator(\n",
+ " M=2**31 - 1, \n",
+ " a=5, \n",
+ " c=0, \n",
+ " seed=12345)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "equipped-jewelry",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "results3 = [prng_standard3.generate() for i in range(10000)]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "german-retention",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.hist(results3)\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel('#draws');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "designed-tribune",
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "plt.scatter(results3[:-1], results3[1:], s=1, marker='.')\n",
+ "plt.xlabel('$x_k$')\n",
+ "plt.ylabel('$x_{k+1}$');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "authorized-carroll",
+ "metadata": {},
+ "source": [
+ "<h3>Box-Muller method (slide 20)</h3>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "seeing-invalid",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "prng_1 = RandomNumberGenerator(\n",
+ " M=2**31 - 1, \n",
+ " a=7**5, \n",
+ " c=0, \n",
+ " seed=12345)\n",
+ "\n",
+ "prng_2 = RandomNumberGenerator(\n",
+ " M=2**31 - 1, \n",
+ " a=7**5, \n",
+ " c=0, \n",
+ " seed=42)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "informal-capture",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def generate_normal():\n",
+ " xi1 = prng_1.generate()\n",
+ " xi2 = prng_1.generate()\n",
+ " r = np.sqrt(-2 * np.log(xi2))\n",
+ " x = r * np.cos(2 * np.pi * xi1)\n",
+ " y = r * np.sin(2 * np.pi * xi1)\n",
+ " return x, y"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "fewer-series",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "results = np.array(\n",
+ " [generate_normal() for i in range(10000)]\n",
+ ").flatten()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "thirty-wallace",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.hist(results, bins=20);"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "mighty-practice",
+ "metadata": {},
+ "source": [
+ "<h3>NumPy has it all...</h3>\n",
+ "\n",
+ "The `numpy` library implements all of these (based on a better behaved variant of the linear congruential generator):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "employed-newark",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.hist(np.random.random(size=10000));\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel('#draws');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "frequent-ottawa",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.hist(np.random.normal(size=10000), bins=20);\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel('#draws');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d1bb3dda",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "<h3>Initalization of longitudinal particle distribution: Interactive Tracking (slide 24)</h3>\n",
+ "\n",
+ "Initialise a bi-Gaussian distribution in the longitudinal phase-space plane for tracking!\n",
+ "\n",
+ "Refer once more to the CERN PS scenario, below transition and in a stationary rf bucket:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "68fb23fd",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m = Machine(gamma_ref=3.13, phi_s=0)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "df3f1231",
+ "metadata": {},
+ "source": [
+ "The length of the rf bucket corresponds to $C/h$:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7ec8568c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "m.circumference / m.harmonic"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "12e5014c",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "We choose an rms bunch length of $\\sigma_z=10\\,$m:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "83008665",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sigma_z = 10"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "36894af5",
+ "metadata": {},
+ "source": [
+ "The corresponding rms momentum difference $\\sigma_{\\Delta p}$ is given through the equilibrium condition (equal Hamiltonian values):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b96b8267",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sigma_deltap = np.sqrt(\n",
+ " 2 * m.p0() / -m.eta(0) * \n",
+ " charge * m.voltage * np.pi * m.harmonic / (beta(gamma(m.p0())) * c * m.circumference**2)\n",
+ ") * sigma_z"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0383dda6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sigma_dp = sigma_deltap / m.p0()\n",
+ "sigma_dp"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "26a48a18",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "plot_hamiltonian(m)\n",
+ "plt.scatter([sigma_z, 0], [0, sigma_dp], marker='*', c='k');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3b0ea4d8",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Back to tracking, out of the box:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9ba73ca6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = 1000\n",
+ "n_turns = 5000"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3eda594b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.random.seed(12345)\n",
+ "\n",
+ "z_ini = np.random.normal(loc=0, scale=sigma_z, size=N)\n",
+ "deltap_ini = np.random.normal(loc=0, scale=sigma_deltap, size=N)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a107fe88",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "z = np.zeros((n_turns, N), dtype=np.float64)\n",
+ "deltap = np.zeros_like(z)\n",
+ "\n",
+ "z[0] = z_ini\n",
+ "deltap[0] = deltap_ini"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "709fcf0d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i_turn in trange(1, n_turns):\n",
+ " z[i_turn], deltap[i_turn] = track_one_turn(z[i_turn - 1], deltap[i_turn - 1], m)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "06a939d9",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Now let's return to the outlined Monte-Carlo approach, analysing the results in terms of statistical moments.\n",
+ "\n",
+ "First the centroid:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4c9d702a",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "-"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "plt.plot(np.mean(z, axis=1))\n",
+ "\n",
+ "plt.xlabel('Turns')\n",
+ "plt.ylabel(r'$\\langle z \\rangle$');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3aa00a3e",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Then the rms beam size (bunch length):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1fdb9767",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(np.std(z, axis=1))\n",
+ "\n",
+ "plt.xlabel('Turns')\n",
+ "plt.ylabel(r'$\\sigma_z$');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ecbcca76",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "<h2>What happened?</h2>\n",
+ "\n",
+ "Let's look at the generated initial distribution of macro-particles:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "71433316",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plot_hamiltonian(m);\n",
+ "plt.scatter(z[0], deltap[0] / m.p0(), marker='.', s=1);"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0142e1f8",
+ "metadata": {},
+ "source": [
+ "$\\implies$ particles have been generated outside the rf bucket! "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "256c1d02",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "<h2>Rejection Sampling Method: Cut particles outside of the separatrix</h2>\n",
+ "\n",
+ "A solution to the problem of generating particles outside the separatrix: <i>reject</i> them at generation!\n",
+ "\n",
+ "$\\leadsto$ a word of caution: this approach modifies the effective distribution function (and the correspondingly generated effective rms values $\\sigma_z, \\sigma_\\delta$ become smaller)! "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "040a0303",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.random.seed(12345)\n",
+ "\n",
+ "z_ini = np.random.normal(loc=0, scale=sigma_z, size=N)\n",
+ "deltap_ini = np.random.normal(loc=0, scale=sigma_deltap, size=N)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "22cf6f11",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "<h2>Rejection Sampling Method: Resample these particles</h2>\n",
+ "\n",
+ "Hamiltonian values of particles outside separatrix are positive (below transition), $\\mathcal{H}>0$ (using the full nonlinear Hamiltonian)!\n",
+ "\n",
+ "NB: Due to the discrete kicks in the finite difference maps, the Hamiltonian is only an approximation for the separatrix: better remain at a few percent distance inside of it!"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f434f623",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "H_safetymargin = 0.05 * hamiltonian(0, 0, m)\n",
+ "\n",
+ "H_values = hamiltonian(z_ini, deltap_ini, m) - H_safetymargin\n",
+ "\n",
+ "while any(H_values >= 0):\n",
+ " mask_bad = H_values >= 0\n",
+ " N_bad = np.sum(mask_bad)\n",
+ " print (N_bad)\n",
+ " # re-initialise bad particles:\n",
+ " z_ini[mask_bad] = np.random.normal(loc=0, scale=sigma_z, size=N_bad)\n",
+ " deltap_ini[mask_bad] = np.random.normal(loc=0, scale=sigma_deltap, size=N_bad)\n",
+ " # re-evaluate rejection condition\n",
+ " H_values = hamiltonian(z_ini, deltap_ini, m) - H_safetymargin"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ee349b25",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Now we should be good to go!"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "be119eb5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plot_hamiltonian(m);\n",
+ "plt.scatter(z_ini, deltap_ini / m.p0(), marker='.', s=1);"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "af1c6162",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Tracking again..."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b812ae0e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "z = np.zeros((n_turns, N), dtype=np.float64)\n",
+ "deltap = np.zeros_like(z)\n",
+ "\n",
+ "z[0] = z_ini\n",
+ "deltap[0] = deltap_ini"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "644cc2bc",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i_turn in trange(1, n_turns):\n",
+ " z[i_turn], deltap[i_turn] = track_one_turn(z[i_turn - 1], deltap[i_turn - 1], m)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ddd7b0f6",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Looking again at the centroid:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c544cbb3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(np.mean(z, axis=1))\n",
+ "\n",
+ "plt.xlabel('Turns')\n",
+ "plt.ylabel(r'$\\langle z \\rangle$');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "580313a8",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Then the rms bunch length:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9e750dee",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "plt.plot(np.std(z, axis=1))\n",
+ "\n",
+ "plt.xlabel('Turns')\n",
+ "plt.ylabel(r'$\\sigma_z$');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1604ad58",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "<p style=\"color: #e6541a;\">$\\implies$ try larger choices of initial $\\sigma_z$ and see where the bunch length evolution saturates.</p>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6ebf9daa",
+ "metadata": {},
+ "source": [
+ "$\\leadsto$ another problem with initialising Gaussian distributions: at larger amplitudes $z$, the small-amplitude approximation with $\\mathcal{H}_\\mathrm{stat,small}$ necessarily breaks down! A Gaussian particle distribution is <b>not</b> in equilibrium for sufficiently large rms values in a nonlinear potential, the particles will <b>filament</b>! (...and the rms emittance will grow, as one can observe in the final equilibrium rms bunch length which is larger than the initial $\\sigma_z$!)\n",
+ "\n",
+ "$\\implies$ generally require full nonlinear Hamiltonian $\\mathcal{H}$ to construct PDF\n",
+ "\n",
+ "$$\\psi(\\mathcal{H})\\propto\\exp\\left(\\cfrac{\\mathcal{H}}{\\mathcal{H}_0}\\right)$$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "eb412faa",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "<h2>RMS Emittance</h2>\n",
+ "\n",
+ "Define statistical <b>rms emittance</b> as in Lecture 2:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a8ca51ed",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def emittance(z, deltap):\n",
+ " N = len(z)\n",
+ " \n",
+ " # subtract centroids\n",
+ " z = z - 1/N * np.sum(z)\n",
+ " deltap = deltap - 1/N * np.sum(deltap)\n",
+ " \n",
+ " # compute Σ matrix entries\n",
+ " z_sq = 1/N * np.sum(z * z)\n",
+ " deltap_sq = 1/N * np.sum(deltap * deltap)\n",
+ " crossterm = 1/N * np.sum(z * deltap)\n",
+ " \n",
+ " # determinant of Σ matrix\n",
+ " epsilon = np.sqrt(z_sq * deltap_sq - crossterm * crossterm)\n",
+ " return epsilon"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1743b226",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "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": "5febd72c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(epsn_z / e)\n",
+ "\n",
+ "plt.xlabel('Turns')\n",
+ "plt.ylabel('$\\epsilon_z$ [eV.s]');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d9417871",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "<p style=\"color: #e6541a;\">$\\implies$ try larger choices of initial $\\sigma_z$ and observe how the emittance evolves.</p>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "canadian-russian",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.9.18"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}