diff --git a/week6/nmap6.ipynb b/week6/nmap6.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..6bee13bc3e00937df4bc03d2ad089cc65509d0c6
--- /dev/null
+++ b/week6/nmap6.ipynb
@@ -0,0 +1,1720 @@
+{
+ "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 6</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 config6 import (np, plt, plot_rfwave, tqdm, trange, \n",
+    "                     beta, gamma, Machine, track_one_turn, \n",
+    "                     charge, mass, emittance, hamiltonian,\n",
+    "                     plot_hamiltonian, plot_rf_overview, \n",
+    "                     plot_dist, plot_mp)\n",
+    "from scipy.constants import m_p, e, c\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": [
+    "<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):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "transparent-webmaster",
+   "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\n",
+    "sigma_deltap"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "exempt-lender",
+   "metadata": {},
+   "source": [
+    "<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"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "polyphonic-wagner",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_hamiltonian(m)\n",
+    "plt.scatter(z, deltap / m.p0(), marker='.', s=1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "extended-surfing",
+   "metadata": {},
+   "source": [
+    "<h3> Compute duration of ramp</h3>\n",
+    "\n",
+    "$(\\Delta\\gamma)_\\mathrm{turn} = \\cfrac{\\Delta E_\\mathrm{tot}}{m_0c^2} = \\cfrac{qV\\sin(\\varphi_s)}{m_0c^2}$\n",
+    "\n",
+    "such that accelerating from $\\gamma_\\mathrm{ref}=3.1$ to $\\gamma_\\mathrm{ref}=27.7$ takes as many turns as\n",
+    "\n",
+    "$n_\\mathrm{turns}=\\cfrac{27.7-3.1}{(\\Delta\\gamma)_\\mathrm{turn}}$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "quiet-rover",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dgamma_per_turn = charge * m.voltage * np.sin(m.phi_s) / (mass * c**2)\n",
+    "n_turns = int(np.ceil((27.7-3.13) / dgamma_per_turn))\n",
+    "n_turns"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "tropical-culture",
+   "metadata": {},
+   "source": [
+    "Record longitudinal emittance during tracking:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "sufficient-negative",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "epsn_z = np.zeros(n_turns, dtype=np.float64)\n",
+    "epsn_z[0] = emittance(z, deltap)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "worth-currency",
+   "metadata": {},
+   "source": [
+    "Let's go tracking!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cooked-consciousness",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for i_turn in trange(1, n_turns):\n",
+    "    z, deltap = track_one_turn(z, deltap, m)\n",
+    "    epsn_z[i_turn] = emittance(z, deltap)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "informal-agent",
+   "metadata": {},
+   "source": [
+    "We have reached the extraction energy of $\\gamma=27.7$:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "front-brunei",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m.gamma_ref"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "given-absence",
+   "metadata": {},
+   "source": [
+    "<h3>How did we do?</h3>\n",
+    "\n",
+    "Check the rms emittance:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "sharp-rescue",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.plot(np.arange(n_turns) / 100000, epsn_z / e)\n",
+    "plt.xlabel('Turns [100k]')\n",
+    "plt.ylabel('$\\epsilon_z$ [eV.s]');"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "typical-collection",
+   "metadata": {},
+   "source": [
+    "$\\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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "furnished-links",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_hamiltonian(m);\n",
+    "plt.scatter(z % 600, deltap / m.p0(), marker='.', s=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "electrical-curve",
+   "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\n",
+    "\n",
+    "m = Machine()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "related-privacy",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for i_turn in trange(1,int(n_turns/100)):\n",
+    "    z, deltap = track_one_turn(z, deltap, m)\n",
+    "    epsn_z[i_turn] = emittance(z, deltap)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "boring-comedy",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(m.gamma_ref)\n",
+    "plot_hamiltonian(m);\n",
+    "plt.scatter(z, deltap / m.p0(), marker='.', s=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "faced-princess",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m.gamma_ref"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "constant-newark",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.sqrt(1/m.alpha_c)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "amber-retreat",
+   "metadata": {},
+   "source": [
+    "<h3>Crossing transition</h3>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "limited-reset",
+   "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\n",
+    "\n",
+    "m = Machine()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "visible-school",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "phi_s_1 = np.pi - m.phi_s # synchronous phase after transition calculated from synchronous phase before transition\n",
+    "\n",
+    "for i_turn in trange(1, n_turns):\n",
+    "    z, deltap = track_one_turn(z, deltap, m)\n",
+    "    \n",
+    "    epsn_z[i_turn] = emittance(z, deltap)\n",
+    "    \n",
+    "    condition = m.gamma_ref > np.sqrt(1/m.alpha_c) # condition to change synchronous phase: gamma above gamma_transition\n",
+    "    \n",
+    "    if condition:\n",
+    "        m.phi_s = phi_s_1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "capable-ladder",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.plot(np.arange(n_turns) / 100000, epsn_z / e)\n",
+    "plt.xlabel('Turns [100k]')\n",
+    "plt.ylabel('$\\epsilon_z$ [eV.s]');"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "velvet-technician",
+   "metadata": {},
+   "source": [
+    "$\\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)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "71dad750-f7c1-4381-9f7e-407a97c07e5b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from PyHEADTAIL.trackers.rf_bucket import RFBucket"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8d6ccb37-ce27-4195-bac1-00f89f3698fa",
+   "metadata": {},
+   "source": [
+    "We define a convenience function to provide `RFBucket` instance given a `Machine` instance:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b0499317-39f1-494c-a31e-0f78bae28411",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def get_pyht_rfbucket(machine):\n",
+    "    m = machine\n",
+    "    deltap_per_turn = charge * m.voltage / (beta(gamma(m.p0())) * c) * np.sin(m.phi_s)\n",
+    "    rfb = RFBucket(m.circumference, m.gamma_ref, mass, charge, [m.alpha_c], deltap_per_turn, [m.harmonic], [m.voltage], [np.pi + m.phi_s])\n",
+    "    # PyHEADTAIL has a different convention for the phase and is offset by pi compared to our lecture\n",
+    "    return rfb"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7a8af807-0422-4664-b509-5aa3d92ff35e",
+   "metadata": {},
+   "source": [
+    "<h3>Visualising the Distributions in the RF Bucket</h3>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bd8760ec-8bc6-4ef3-ab92-f4128d60a5de",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m = Machine(gamma_ref=3.13)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9c8820d4-2aaf-4cfb-976f-6765b6827395",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rfb = get_pyht_rfbucket(m)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5356375e-6127-4cff-8451-7e8c534621a0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from PyHEADTAIL.particles.rfbucket_matching import (ThermalDistribution, WaterbagDistribution, ParabolicDistribution)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "884fbafd-c180-4942-97e5-4abb2078e45b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sigma_z = 8"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "63a1dc42-5f32-44af-8529-87ac778762b9",
+   "metadata": {},
+   "source": [
+    "Computing the (initial) guess for $\\mathcal{H}_0$ based on the small-amplitude approximation:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f4cf9828-321f-44a4-994d-1a5977634959",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "H0 = rfb.guess_H0(sigma_z, from_variable='sigma')\n",
+    "H0"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "58189c7a-755e-4202-8963-450df183ef18",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_dist(ThermalDistribution, rfb, H0=0.05*H0);"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7b5a10a6-3dac-41a0-8cd4-2e5944803f16",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_dist(WaterbagDistribution, rfb, H0=2*H0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "27b62728-ba57-4ef9-9f01-e78dab20da03",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_dist(ParabolicDistribution, rfb, H0=H0*2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d1f9bbc3-5614-45de-9a86-4218daa1eef1",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d82d4bfc-2b95-46fd-abc8-ef7fee8ad21e",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "189c0360-5951-4a3f-b5fb-23aef3be6699",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cb45ec0d-9701-4b32-abf3-b1a50c89553e",
+   "metadata": {},
+   "source": [
+    "Matching algorithm implemented in the `RFBucketMatcher` class in `PyHEADTAIL`\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "82aad28a-6034-4a3c-97cd-586649bc787c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from PyHEADTAIL.particles.generators import RFBucketMatcher"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dbe611ba-b719-4882-b37e-9afe74c0b014",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rfb_matcher = RFBucketMatcher(rfb, ThermalDistribution, sigma_z=sigma_z)\n",
+    "rfb_matcher.integrationmethod = 'cumtrapz'"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d1955099-6310-4904-a4a0-7165b6ab72fb",
+   "metadata": {},
+   "source": [
+    "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",
+    "print(\"Small-amplitude approximation: \" + str(rfb.guess_H0(sigma_z, from_variable='sigma')))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f9d774c3-3dbb-4abb-9ed6-3b421c720f98",
+   "metadata": {},
+   "source": [
+    "$\\implies$ For smaller target $\\hat{\\sigma}_z$ the values are closer together."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "19929620-574d-41f1-bf8b-88b3d4fd10e8",
+   "metadata": {},
+   "source": [
+    "Let's have a look at the generated macro-particle distribution:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "98bb9062-c240-4fb9-8f14-e07d6d4d77a1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_mp(z, delta, rfb);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "23d74231-fd49-44d9-a127-501e0bf2b99a",
+   "metadata": {},
+   "source": [
+    "Does the rms bunch length of the macro-particle distribution match the chosen target $\\hat{\\sigma}_z$?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b92822f5-582e-47bd-a46b-144e040a8043",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(\"RMS bunch length: \" + str(np.std(z)))\n",
+    "print(\"Target sigma_z: \" + str(sigma_z))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "acc29308-7c1f-4afa-9575-ed0189255f9b",
+   "metadata": {},
+   "source": [
+    "<h2>Root-finding</h2>\n",
+    "\n",
+    "<h3>Let's use scipy.optimize ...</h3>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c8d090fd-eb15-47f6-9391-11e2de309d32",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipy.optimize import brentq, newton"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9807ca23-733f-4440-925f-56fe27b873bd",
+   "metadata": {},
+   "source": [
+    "At what $x$ does $\\cfrac{1}{\\sqrt{x}+1}$ take the value `0.4`?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e676a536-1fab-4762-a7a6-61fd173928b0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def f(x):\n",
+    "    return 1 / (np.sqrt(x) + 1) - 0.4"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e5184f7d-a7f1-4d7a-b304-e9a877df8b20",
+   "metadata": {},
+   "source": [
+    "Brent-Dekker algorithm with interval $x\\in[0, 4]$:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f5fb0a95-f8d4-496f-adbe-1e3eb1bbd87f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "brentq(f, 0, 4)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2812695d-dc9b-4548-8203-a7a6fb4e4e32",
+   "metadata": {},
+   "source": [
+    "<h3>... and Newton's Secant method</h3>"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "384e3c1f-916a-4938-bee4-eb5d30f3e806",
+   "metadata": {},
+   "source": [
+    "Secant method with initial values $x_0=0$, $x_1=10^{-4}$:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a81bdb20-73d6-4a38-9b72-b8f12e89a03a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "newton(f,0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b3a79de2-4c83-44a2-9af0-8a4de02f71c1",
+   "metadata": {},
+   "source": [
+    "<h3>How does the function look like?</h3>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a7a4b9ac-d8b6-4173-a6cc-71677a83e364",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "x = np.linspace(0, 4, 1000)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e98a563d-618c-4ced-9104-86491335242e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.plot(x, f(x))\n",
+    "plt.axhline(0, c='k', lw=2)\n",
+    "plt.xlabel('$x$')\n",
+    "plt.ylabel('$f(x)$');"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9f85826c-5f43-42cc-b8cc-3a10d1192801",
+   "metadata": {},
+   "source": [
+    "<h3>Implementing the Secant method</h3>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2d9a869d-2580-4674-9746-97e47deb177f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def secant_method(f, x0, x1, iterations, rtol=2e-12):\n",
+    "    \"\"\"Return the root calculated using the secant method.\"\"\"\n",
+    "    for i in range(iterations):\n",
+    "        x2 = x1 - f(x1) * (x1 - x0) / float(f(x1) - f(x0))\n",
+    "        if np.abs(x2 - x1) / x1 < rtol:\n",
+    "            break\n",
+    "        x0, x1 = x1, x2\n",
+    "    return x2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e0e38a3b-5b40-42cb-902b-14deba599c1e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "secant_method(f, 0, 1, 100)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d845eda7-42f8-4e85-8b8c-edefc28eb0db",
+   "metadata": {},
+   "source": [
+    "(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):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bbbcee55-3dd5-407a-8f59-d44b974578e0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m = Machine(phi_s=0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fb8c6faf-d25d-474b-adca-fa4cb50263a7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sigma_z = 13.5"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "45d704d6-a229-40bf-b83a-c59a319f1932",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def generate_gaussian_in_rfbucket(N, sigma_z, machine, seed=12345, margin=0.05):\n",
+    "    '''Generate a bi-Gaussian distribution with N macro-particles,\n",
+    "    rms bunch length sigma_z and a matched sigma_deltap via the\n",
+    "    machine settings.\n",
+    "    '''\n",
+    "    np.random.seed(seed)\n",
+    "    m = machine\n",
+    "    \n",
+    "    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\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)\n",
+    "    \n",
+    "    H_safetymargin = margin * hamiltonian(0, 0, m)\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\n",
+    "    \n",
+    "    return z_ini, deltap_ini"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1eb0b492-b2fc-4ecd-8de0-55fa81dd81d8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "N = 10000\n",
+    "n_turns = 5000"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f9021cbc-4e26-4847-88bd-9fda118c36ed",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "z_ini, deltap_ini = generate_gaussian_in_rfbucket(N, sigma_z, m)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "423f9dc9-c246-42da-861b-e0255dc34f4d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_mp(z_ini, deltap_ini / m.p0(), rfb=get_pyht_rfbucket(m), n_bins=20);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a3c290ef-e08e-43d3-b28a-c788ca72b5b1",
+   "metadata": {},
+   "source": [
+    "Track the bi-Gaussian distribution..."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4d49d207-236d-434b-9097-29941e8a5ba1",
+   "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": "df20e1d7-7748-4c96-be03-0aec3cdf1322",
+   "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": "607dd837-74b5-4258-9f70-5fb4098b2858",
+   "metadata": {},
+   "source": [
+    "The rms emittance evolution:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ae3aa012-7673-4e4a-9176-37a6e8185161",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "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)\n",
+    "ylim_d = 1.1 * np.max(np.abs(ylim_m - 4 * np.pi * epsn_z / e))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e7121496-7d6e-4e0c-bfe4-fe93d0e24da6",
+   "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": "b83c18bd-7a7a-4411-b6e5-14950c44fe8e",
+   "metadata": {},
+   "source": [
+    "$\\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)$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c30eff4c-3b47-4758-9b69-f015eb8ffa44",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rfb = get_pyht_rfbucket(m)\n",
+    "\n",
+    "rfb_matcher = RFBucketMatcher(rfb, ThermalDistribution, sigma_z=sigma_z)\n",
+    "rfb_matcher.integrationmethod = 'cumtrapz'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0a60a6f1-876e-4f9d-b0ae-dae66b218188",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "z_ini, delta_ini, _, _ = rfb_matcher.generate(N)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "382bc26d-321d-4dbf-bb65-dbbe49827bbd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "deltap_ini = delta_ini * m.p0()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8ef7c870-cd39-4fcd-9c23-160e3980ccb2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_mp(z_ini, delta_ini, rfb, n_bins=20);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5f818b3e-df01-4551-b519-14436b114c86",
+   "metadata": {},
+   "source": [
+    "Track the matched thermal distribution..."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b160b0a7-4ad2-46fd-8a81-7b6ab778917b",
+   "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": "e28afb33-0562-4128-ab00-5aa24608509f",
+   "metadata": {
+    "scrolled": true
+   },
+   "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": "a9a97d61-7998-465d-afbe-d3f6c00acddb",
+   "metadata": {},
+   "source": [
+    "The rms emittance evolution:"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fae06e09-37bf-4032-8f6f-43d44cb321b6",
+   "metadata": {},
+   "source": [
+    "The rms emittance evolution:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4e1e39a1-5f92-47af-840a-7ca55f48adc4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "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}$)!"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4a335edf-1cb0-4a9b-9499-3679ce396e5d",
+   "metadata": {},
+   "source": [
+    "<h2>Physical emittance growth: Dipole injection mismatch</h2>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dc77bb22-52f8-4273-bb53-3a09c3e7e324",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m = Machine(phi_s=0)\n",
+    "\n",
+    "sigma_z = 8\n",
+    "\n",
+    "z_ini, deltap_ini = generate_gaussian_in_rfbucket(N, sigma_z, m, margin=0.15)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "baef42d4-f6e7-4af1-a6d6-ffba009ffcbb",
+   "metadata": {},
+   "source": [
+    "Simulate a dipole injection mismatch (e.g. when the rf phase is not well synchronised between the injector and the synchrotron):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "471d66b4-0574-4d6f-a167-4bd7ba001156",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "z_ini -= 0.5 * sigma_z"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "53e01f17-642f-43bf-ac8c-30c0cb7d9c41",
+   "metadata": {},
+   "source": [
+    "4 meter mismatch in $z$ correspond to a phase mismatch of 16 degree:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b952ef70-5a26-4354-b2d9-21c652db06cc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "4 / (m.circumference / m.harmonic) * 360"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "39c3dce3-1df1-43df-828f-653860c2d656",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_mp(z_ini, deltap_ini / m.p0(), rfb=get_pyht_rfbucket(m), n_bins=20);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "95944d0d-3053-4639-a9ef-fb1bd9ed1126",
+   "metadata": {},
+   "source": [
+    "$\\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:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5beb74ef-1408-4738-9117-c82aa6fae6d4",
+   "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": "5d4168ce-213a-4695-ab0a-55e91f42c02d",
+   "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": "0d5c6857-6693-4dae-85b4-ac7dd8fb1071",
+   "metadata": {},
+   "source": [
+    "<h3>Centroid results</h3>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3677d73a-542a-4965-8f16-38e60b4e5e8c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.plot(np.mean(z, axis=1))\n",
+    "\n",
+    "plt.xlabel('Turns')\n",
+    "plt.ylabel(r'$\\langle z \\rangle$ [m]');"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e0334d3b-bd04-49e0-bd90-4d8fc8e8bea9",
+   "metadata": {},
+   "source": [
+    "$\\implies$ exponential decay of the initial offset (due to the non-linearity of the rf bucket)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d7107e24-9854-40d1-9523-95f5d1e24822",
+   "metadata": {},
+   "source": [
+    "<h3>RMS bunch length results</h3>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "14148630-7eff-49ca-b1e4-4dacc1cd97e0",
+   "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": "6e0b8011-7571-4280-8b6c-a989b885b094",
+   "metadata": {},
+   "source": [
+    "$\\implies$ saturation of the rms bunch length growth"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "645f0d2b-344a-4830-b015-ff3843c50e78",
+   "metadata": {},
+   "source": [
+    "<h3>RMS emittance results</h3>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d4e0f5d1-8833-4595-bacb-70023f184a7c",
+   "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": "e0f375a6-28ad-4c39-813f-c85cb6021f73",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.plot(epsn_z / e)\n",
+    "\n",
+    "plt.xlabel('Turns')\n",
+    "plt.ylabel('$\\epsilon_z$ [eV.s]');"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7788ba39-d406-46d2-a5f0-87f2762eb1db",
+   "metadata": {},
+   "source": [
+    "$\\implies$ in this example, 10% emittance growth as a result of the 4 meter injection offset."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "27c3287f-29da-46ff-bed8-e4808f3de18c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_mp(z[-1], deltap[-1] / m.p0(), rfb=get_pyht_rfbucket(m), n_bins=40);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d40cdc11-5ebb-4090-976e-f99aab567ab5",
+   "metadata": {},
+   "source": [
+    "$\\implies$ the filamentation of the macro-particle distribution is clearly visible!"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6523b6d7-2984-4eb8-b8d3-2968a93cba41",
+   "metadata": {},
+   "source": [
+    "<h2>Physical emittance growth: Quadrupole injection mismatch</h2>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ea3648e7-7b00-4d17-b0bb-6d0f9fa41b45",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "m = Machine(phi_s=0)\n",
+    "\n",
+    "sigma_z = 8\n",
+    "\n",
+    "z_ini, deltap_ini = generate_gaussian_in_rfbucket(N, sigma_z, m, margin=0.15)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "09a54392-f60b-4563-9c10-793328382fd9",
+   "metadata": {},
+   "source": [
+    "Simulate a quadrupole injection mismatch (e.g. when the rf voltage (rf bucket height) is not matched between the injector and the synchrotron):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a77e1f8c-89e6-44a7-b7ed-ef37ca41476f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "deltap_ini *= 0.5"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "edba9fe9-61e1-4343-906c-abcb62a76b2e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_mp(z_ini, deltap_ini / m.p0(), rfb=get_pyht_rfbucket(m), n_bins=20);"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "19ba693a-f643-47fb-ad0b-47632f7f4940",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0360ab47-0c87-44ea-a45f-fd9b1700b531",
+   "metadata": {},
+   "source": [
+    " $\\implies$ note the squeezed rms momentum spread, the contours of the macro-particle density are no longer matched to the Hamiltonian contours."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5d37bd16-d02d-4d62-a19b-ffeb5c948c9c",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fd445783-47b6-4315-a067-a5d516fe4424",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cc5856eb-9c81-44ba-9457-2a18cc706937",
+   "metadata": {},
+   "source": [
+    "Tracking the mismatched distribution of macro-particles:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2ad36092-7775-49e4-b1a4-3796e097cf00",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "42df2535-76cf-4647-b217-ff65e1ca0ab2",
+   "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": "8550522f-11b0-4b53-bf27-e44e1f3f1696",
+   "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": "faa79c31-bcab-4ab8-b198-b3dd8114e897",
+   "metadata": {},
+   "source": [
+    "<h3>Centroid results</h3>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2a1d11bc-c25c-4f61-84a4-0b1aad5be9e2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.plot(np.mean(z, axis=1))\n",
+    "\n",
+    "plt.xlabel('Turns')\n",
+    "plt.ylabel(r'$\\langle z \\rangle$ [m]');"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f5b846ea-c561-4d05-96b7-86a6a18ff723",
+   "metadata": {},
+   "source": [
+    "$\\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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "02cb8496-e4d3-4131-bda8-f88957c8ff75",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_mp(z[-1], deltap[-1] / m.p0(), rfb=get_pyht_rfbucket(m), n_bins=40);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "58404467-c2fd-491c-87f5-42848bd0e7d9",
+   "metadata": {},
+   "source": [
+    "$\\implies$ again, the filamentation of the macro-particle distribution is clearly visible!"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.20"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}