diff --git a/week7/nmap7.ipynb b/week7/nmap7.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..91ae034faca6415c2b398320a0d47d872e215635
--- /dev/null
+++ b/week7/nmap7.ipynb
@@ -0,0 +1,1204 @@
+{
+ "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 7</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 config7 import (np, plt)\n",
+ "from scipy.constants import m_p, e, c\n",
+ "%matplotlib inline"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "58404467-c2fd-491c-87f5-42848bd0e7d9",
+ "metadata": {},
+ "source": [
+ "<h2>Betatron Matrices</h2>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "734808ec-8325-4524-bb69-f80ffdcdeb88",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def M_drift(L):\n",
+ " return np.array([\n",
+ " [1, L],\n",
+ " [0, 1]\n",
+ " ])\n",
+ "\n",
+ "def M_dip_x(L, rho0):\n",
+ " return np.array([\n",
+ " [np.cos(L / rho0), rho0 * np.sin(L / rho0)],\n",
+ " [-1 / rho0 * np.sin(L / rho0), np.cos(L / rho0)]\n",
+ " ])\n",
+ "\n",
+ "def M_dip_y(L, rho0):\n",
+ " return M_drift(L)\n",
+ "\n",
+ "def M_quad_x(L, k):\n",
+ " ksq = np.sqrt(k + 0j)\n",
+ " return np.array([\n",
+ " [np.cos(ksq * L), 1 / ksq * np.sin(ksq * L)],\n",
+ " [-ksq * np.sin(ksq * L), np.cos(ksq * L)]\n",
+ " ]).real\n",
+ "\n",
+ "def M_quad_y(L, k):\n",
+ " ksq = np.sqrt(k + 0j)\n",
+ " return np.array([\n",
+ " [np.cosh(ksq * L), 1 / ksq * np.sinh(ksq * L)],\n",
+ " [ksq * np.sinh(ksq * L), np.cosh(ksq * L)]\n",
+ " ]).real"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0f9532eb-7616-4c5f-8da8-b6b567c08fa3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def track(M, u, up):\n",
+ " '''Apply M to each individual [u;up] vectors value.'''\n",
+ " return np.einsum('ij,...j->i...', M, np.vstack((u, up)).T)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0381ff44-6119-48c2-ac61-6090afc89239",
+ "metadata": {},
+ "source": [
+ "<h2>Thin Sextupole Kick</h2>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "cdd4c0b3-4476-4570-a8fe-403f8899f457",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def track_sext_4D(x, xp, y, yp, mL):\n",
+ " xp += 0.5 * mL * (y * y - x * x)\n",
+ " yp += mL * x * y\n",
+ " return x, xp, y, yp"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c6fa6d11-def6-43ca-a779-26c4ee6ea6a3",
+ "metadata": {},
+ "source": [
+ "<h2>Simulation Examples</h2>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b14536da-d463-40a7-a860-a2848c40cd0b",
+ "metadata": {},
+ "source": [
+ "<h3>1. Simulating a drift:</h3>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "397f2880-55a6-4982-87cd-a33658c59b02",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.random.seed(12345)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d2bf96b4-d127-4b20-bfc7-304cdd34296d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = 100\n",
+ "sig_x = 5e-3\n",
+ "sig_xp = 2e-3\n",
+ "\n",
+ "x = np.random.normal(0, sig_x, N)\n",
+ "xp = np.random.normal(0, sig_xp, N)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ff3193f4-0fdf-4424-ab75-06a30aa9ca4b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ds = 0.01\n",
+ "\n",
+ "D = M_drift(ds)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "778ff776-2c57-4b9a-bd67-e1c288d04a0a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.scatter(x, xp, c='C0', s=10, marker='.')\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel(\"$x'$\");"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "db556ca7-1ed8-4738-a239-7cf81d452397",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for s in np.arange(-1, 1, ds):\n",
+ " x, xp = track(D, x, xp)\n",
+ " plt.scatter(np.ones(N) * s, x, c='C0', s=1, marker='.')\n",
+ "plt.xlabel('$s$')\n",
+ "plt.ylabel('$x$');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "19f4b1ad-48a8-484d-84ec-19a067ebd2fd",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.scatter(x, xp, c='C0', s=10, marker='.')\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel(\"$x'$\");"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "889afdf2-0d85-45a5-bf07-a3d9b9adef01",
+ "metadata": {},
+ "source": [
+ "Same simulation again with correlated $x$, $x'$:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5bc67b5c-1b1f-4237-a618-f1f87e142a78",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.random.seed(12345)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "aac570be-19ca-4b54-91f1-c5c9be78bdef",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = 100\n",
+ "sig_x = 5e-3\n",
+ "sig_xp = 2e-3\n",
+ "\n",
+ "x = np.random.normal(0, sig_x, N)\n",
+ "xp = np.random.normal(0, sig_xp / 2, N) - x * sig_x / sig_xp * 0.4"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1171d3cb-2d00-43ae-9216-782f5763a2cd",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ds = 0.01\n",
+ "\n",
+ "D = M_drift(ds)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "58e4a7d2-eeac-40a4-9d75-a9beb83fdb59",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.scatter(x, xp, c='C0', s=10, marker='.')\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel(\"$x'$\");"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "18bf22ea-3ad4-4cac-b49f-049f8b1b1c38",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for s in np.arange(-1, 1, ds):\n",
+ " x, xp = track(D, x, xp)\n",
+ " plt.scatter(np.ones(N) * s, x, c='C0', s=1, marker='.')\n",
+ "plt.xlabel('$s$')\n",
+ "plt.ylabel('$x$');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "273234c6-a1e0-4234-ab9c-6064953b0fbf",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.scatter(x, xp, c='C0', s=10, marker='.')\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel(\"$x'$\");"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "be0e88c6-2d85-4691-b99b-edb38c16aea5",
+ "metadata": {},
+ "source": [
+ "$\\implies$ particles move (drift) for- and backward along $x$ depending on their momentum (angle) $x'$!"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e5524471-e30c-4b50-ba47-1d5b370cecb4",
+ "metadata": {},
+ "source": [
+ "<h3>Simulating a quadrupole in focusing plane:</h3>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3049203e-f542-47bf-9621-c59e18643502",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.random.seed(12345)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "059135e6-d6e7-497e-b7d0-2fd4d9297586",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = 100\n",
+ "sig_x = 5e-3\n",
+ "sig_xp = 2e-3\n",
+ "\n",
+ "x = np.random.normal(0, sig_x, N)\n",
+ "xp = np.random.normal(0, sig_xp, N)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9be96dec-e05a-4465-bd94-e55b4b1eae86",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ds = 0.01\n",
+ "\n",
+ "Qx = M_quad_x(ds, 10)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "922eebd2-8388-4384-8d41-4fc829ba485a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.scatter(x, xp, c='C0', s=10, marker='.')\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel(\"$x'$\");"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5edb87a1-5be7-444a-8a95-a143f1f10b84",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for s in np.arange(-1, 1, ds):\n",
+ " x, xp = track(Qx, x, xp)\n",
+ " plt.scatter(np.ones(N) * s, x, c='C0', s=1, marker='.')\n",
+ "plt.xlabel('$s$')\n",
+ "plt.ylabel('$x$');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "962c27f1-dab1-41ca-a8d8-0544cf443003",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.scatter(x, xp, c='C0', s=10, marker='.')\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel(\"$x'$\");"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8a81b0ab-45be-42fd-8415-f7667dc1f09b",
+ "metadata": {},
+ "source": [
+ "<h3>Simulating a quadrupole in defocusing plane:</h3>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "00dee1a9-b23d-4cee-bd89-8f0b78b7dae2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.random.seed(12345)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2f03e443-fa18-40e8-9cbe-bf46d29292d8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = 100\n",
+ "sig_x = 5e-3\n",
+ "sig_xp = 2e-3\n",
+ "\n",
+ "x = np.random.normal(0, sig_x, N)\n",
+ "xp = np.random.normal(0, sig_xp, N)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ad35d573-53fd-40d1-bc50-48596f77d333",
+ "metadata": {},
+ "source": [
+ "Note the negative sign for the $k$ to obtain defocusing in the horizontal plane:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ad441245-5e29-4b96-b8f0-96ef70776fd3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ds = 0.01\n",
+ "\n",
+ "Qx = M_quad_x(ds, -10)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9d9d1f03-2985-4690-b670-ae2d36e0c1cb",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.scatter(x, xp, c='C0', s=10, marker='.')\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel(\"$x'$\");"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d2a6b561-d8c4-4447-b16c-4a71157c73e9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for s in np.arange(-1, 1, ds):\n",
+ " x, xp = track(Qx, x, xp)\n",
+ " plt.scatter(np.ones(N) * s, x, c='C0', s=1, marker='.')\n",
+ "plt.xlabel('$s$')\n",
+ "plt.ylabel('$x$');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "540af12a-7a57-4a4e-a9d5-a113e2ab2d68",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.scatter(x, xp, c='C0', s=10, marker='.')\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel(\"$x'$\");"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "439555a8-0b0d-49a8-b561-1cb0045735b2",
+ "metadata": {},
+ "source": [
+ "$\\implies$ the problem of constant continous quadrupole fields lies in that one plane focuses while the other plane defocuses!"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a8d3958b-f6ab-4775-bd76-4cc33aedeeaa",
+ "metadata": {},
+ "source": [
+ "<h3>4. Simulating a FODO cell:</h3>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "bbf3246b-6564-4b79-963f-1abe960b8fe1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.random.seed(12345)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5857bc0d-cea9-44d3-9c39-b74b3ae86a63",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = 100\n",
+ "sig_x = 5e-3\n",
+ "sig_xp = 3e-4\n",
+ "\n",
+ "x = np.random.normal(0, sig_x, N)\n",
+ "xp = np.random.normal(0, sig_xp, N)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3f1c7998-4938-4682-8bde-e820f77cb1d3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ds = 0.1\n",
+ "k = 0.2\n",
+ "\n",
+ "D = M_drift(ds)\n",
+ "Qfx = M_quad_x(ds, k)\n",
+ "Qdx = M_quad_x(ds, -k)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1a6760e4-f0e0-4d63-abd6-ee18220e99ea",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.scatter(x, xp, c='C0', s=10, marker='.')\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel(\"$x'$\");"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b9edc3dc-be1a-441a-b6a2-e806f7fbe73d",
+ "metadata": {},
+ "source": [
+ "We assume a total FODO cell length of 10m and a length of each quadrupole magnet of 1m.\n",
+ "Tracking the FODO cell in the horizontal plane, starting from the center of the focusing quadrupole (=horizontally focusing!):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2612d244-8c07-42df-9d04-94db7a391a52",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# 1/2 focusing quad\n",
+ "for s in np.arange(0, 0.5001, ds):\n",
+ " x, xp = track(Qfx, x, xp)\n",
+ " plt.scatter(np.ones(N) * s, x, c='C0', s=1, marker='.')\n",
+ "# drift\n",
+ "for s in np.arange(0.5, 4.5001, ds)[1:]:\n",
+ " x, xp = track(D, x, xp)\n",
+ " plt.scatter(np.ones(N) * s, x, c='C0', s=1, marker='.')\n",
+ "# defocusing quad\n",
+ "for s in np.arange(4.5, 5.5001, ds)[1:]:\n",
+ " x, xp = track(Qdx, x, xp)\n",
+ " plt.scatter(np.ones(N) * s, x, c='C0', s=1, marker='.')\n",
+ "# drift\n",
+ "for s in np.arange(5.5, 9.5001, ds)[1:]:\n",
+ " x, xp = track(D, x, xp)\n",
+ " plt.scatter(np.ones(N) * s, x, c='C0', s=1, marker='.')\n",
+ "# 1/2 focusing quad\n",
+ "for s in np.arange(9.5, 10.0001, ds)[1:]:\n",
+ " x, xp = track(Qfx, x, xp)\n",
+ " plt.scatter(np.ones(N) * s, x, c='C0', s=1, marker='.')\n",
+ "plt.xlabel('$s$')\n",
+ "plt.ylabel('$x$');"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "edf21022-fd62-4faa-8f43-3a874af87a1b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.scatter(x, xp, c='C0', s=10, marker='.')\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel(\"$x'$\");"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9d5e3a0b-aade-491b-b76a-c566dfc926ea",
+ "metadata": {},
+ "source": [
+ "What about the vertical plane now? The quadrupoles have their function inverted, a horizontally focusing quadrupole defocuses in the vertical plane, so the same lattice looks like \"D-O-F-O\" with respect to the vertical plane:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8c812ee8-e125-4511-9d77-18d4dd9eaeff",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sig_y = 5e-3\n",
+ "sig_yp = 3e-4\n",
+ "\n",
+ "y = np.random.normal(0, sig_y, N)\n",
+ "yp = np.random.normal(0, sig_yp, N)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4fbf3925-8597-4bc1-a914-77b1c4dcd824",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Qfy = M_quad_y(ds, k)\n",
+ "Qdy = M_quad_y(ds, -k)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "df82b0c1-608e-48c1-b6a2-66925175ca07",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# 1/2 vertically defocusing quad (horizontally focusing, so \"Qf\")\n",
+ "for s in np.arange(0, 0.5001, ds):\n",
+ " y, yp = track(Qfy, y, yp)\n",
+ " plt.scatter(np.ones(N) * s, y, c='C0', s=1, marker='.')\n",
+ "# drift\n",
+ "for s in np.arange(0.5, 4.5001, ds)[1:]:\n",
+ " y, yp = track(D, y, yp)\n",
+ " plt.scatter(np.ones(N) * s, y, c='C0', s=1, marker='.')\n",
+ "# defocusing quad (horizontally defocusing, so \"Qd\")\n",
+ "for s in np.arange(4.5, 5.5001, ds)[1:]:\n",
+ " y, yp = track(Qdy, y, yp)\n",
+ " plt.scatter(np.ones(N) * s, y, c='C0', s=1, marker='.')\n",
+ "# drift\n",
+ "for s in np.arange(5.5, 9.5001, ds)[1:]:\n",
+ " y, yp = track(D, y, yp)\n",
+ " plt.scatter(np.ones(N) * s, y, c='C0', s=1, marker='.')\n",
+ "# 1/2 focusing quad\n",
+ "for s in np.arange(9.5, 10.0001, ds)[1:]:\n",
+ " y, yp = track(Qfy, y, yp)\n",
+ " plt.scatter(np.ones(N) * s, y, c='C0', s=1, marker='.')\n",
+ "plt.xlabel('$s$')\n",
+ "plt.ylabel('$y$');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "631ece08-d991-4155-826d-ab4af9a81c8e",
+ "metadata": {},
+ "source": [
+ "$\\implies$ can ensure quasi-harmonic motion in both (!) transverse planes! Transverse confinement of beam by alternating-gradient (AG) focusing! This is the principle behind synchrotrons!\n",
+ "\n",
+ "How does this look like over long time scales? Let us build the one-cell matrix and track for many cells:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "be2815b2-c594-46ed-bb2b-398ce24f595f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "M_cell = M_quad_x(0.5, k) # 1/2 focusing quad\n",
+ "M_cell = M_cell.dot(M_drift(4)) # drift\n",
+ "M_cell = M_cell.dot(M_quad_x(1, -k)) # defocusing quad\n",
+ "M_cell = M_cell.dot(M_drift(4)) # drift\n",
+ "M_cell = M_cell.dot(M_quad_x(0.5, k)) # 1/2 focusing quad"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3bece1ba-2df3-4faa-b4c2-7d0be633947e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_cells = 100"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "fc50b41e-c032-4b0d-8efe-8f6f66b8f949",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.random.seed(12345)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "82993ffc-0ab9-4c9e-be54-7d90c6458fef",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = 100\n",
+ "sig_x = 5e-3\n",
+ "sig_xp = 3e-4\n",
+ "\n",
+ "x = np.random.normal(0, sig_x, N)\n",
+ "y = np.random.normal(0, sig_xp, N)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "abf4158c-668b-4475-89ea-cd60ac80187a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i in range(n_cells):\n",
+ " x, xp = track(M_cell, x, xp)\n",
+ "\n",
+ " plt.scatter(np.ones(N) * i, x, c='C0', s=1, marker='.')\n",
+ " plt.scatter([i], [x[-1]], c='r', s=10, marker='.')\n",
+ "plt.xlabel('Cells')\n",
+ "plt.ylabel('$x$');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5922d0ba-a635-44e3-be5a-1f85eb9924dd",
+ "metadata": {},
+ "source": [
+ "$\\implies$ we observe regular motion, amplitudes remain bound! It looks like the magnet configuration is stable and the beam is well confined!\n",
+ "\n",
+ "What about the phase-space trajectories at this position in the lattice (a so-called Poincaré section)?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "952fc008-ceb0-47b8-916f-edc22b95be33",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i in range(n_cells):\n",
+ " x, xp = track(M_cell, x, xp)\n",
+ " plt.scatter(x[::10], xp[::10], c='C0', s=10, marker='.')\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel(\"$x'$\")\n",
+ "plt.gca().set_aspect(np.diff(plt.xlim()) / np.diff(plt.ylim()));"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "cdc8f2e0-289a-4abf-8901-6735a6033cb9",
+ "metadata": {},
+ "source": [
+ "$\\implies$ the circles indicate linear bound motion!"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "93edce31-4b3b-4294-a3b2-e1c63ca915f3",
+ "metadata": {},
+ "source": [
+ "<h3>5. Simulating a FODO cell with increasingly strong $k$:</h3>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "658c32c7-a96e-479d-9288-5845f5e4cd99",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "k = 0.431"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6907cbf8-734c-4180-82e3-2b693e8f669e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "M_cell = M_quad_x(0.5, k) # 1/2 focusing quad\n",
+ "M_cell = M_cell.dot(M_drift(4)) # drift\n",
+ "M_cell = M_cell.dot(M_quad_x(1, -k)) # defocusing quad\n",
+ "M_cell = M_cell.dot(M_drift(4)) # drift\n",
+ "M_cell = M_cell.dot(M_quad_x(0.5, k)) # 1/2 focusing quad"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5ae035af-5616-4627-8bc6-21418613eae9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_cells = 100"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c52dd957-66de-4aa0-8ae8-d4d7daa36555",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = 100\n",
+ "sig_x = 5e-3\n",
+ "sig_xp = 3e-4\n",
+ "\n",
+ "x = np.random.normal(0, sig_x, N)\n",
+ "xp = np.random.normal(0, sig_xp, N)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b3cdb0cb-6826-4fd6-a8a5-1f3524c0bb4b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i in range(n_cells):\n",
+ " x, xp = track(M_cell, x, xp)\n",
+ " plt.scatter(np.ones(N) * i, x, c='C0', s=1, marker='.')\n",
+ " plt.scatter([i], [x[-1]], c='r', s=10, marker='.')\n",
+ "plt.xlabel('Cells')\n",
+ "plt.ylabel('$x$');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e05df3a7-3bcd-4fb8-b6ce-3a6335d8cb1c",
+ "metadata": {},
+ "source": [
+ "$\\implies$ motion becomes unstable! Is the one-cell matrix a \"valid\" (symplectic) betatron matrix?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b0caacc4-b24e-4543-a7d9-9ea62a6b76bc",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.linalg.det(M_cell)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4bf0efe8-0181-4993-a08e-0700ee2711db",
+ "metadata": {},
+ "source": [
+ "$\\implies$ the matrix obeys $\\det(\\mathcal{M})=1$ and is thus symplectic. But what about the eigenvalues?\n",
+ "\n",
+ "Solve the characteristic polynomial of the one-cell matrix, $\\det(\\mathcal{M}-\\lambda\\mathbb{1})=0$ for $\\lambda$:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "51533d43-b457-4368-b070-bff2973a74dd",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.linalg.eigvals(M_cell)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f8944c88-5aa1-45e2-96d4-41b1a0ea27fd",
+ "metadata": {},
+ "source": [
+ "$\\implies$ we find one $|\\lambda|>1$! If one absolute eigenvalue becomes larger than unity, the magnet configuration becomes unstable! That explains the instability (exponential divergence) here! Equivalently one finds $|\\mathrm{Tr}(\\mathcal{M})|>2$:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8e92e564-514c-4587-a7bc-0c8dce9146c9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.trace(M_cell)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "4b8c2711-50ef-4df9-9f93-073232be83ff",
+ "metadata": {},
+ "source": [
+ "What happens to a single particle in phase space in the Poincaré section?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e12d7300-a398-4913-b398-f457dd138d44",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i in range(10):\n",
+ " x, xp = track(M_cell, x, xp)\n",
+ " plt.scatter(x[0], xp[0], c='C0', s=10, marker='.')\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel(\"$x'$\")\n",
+ "plt.gca().set_aspect(np.diff(plt.xlim()) / np.diff(plt.ylim()));"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b3e86ef0-fb42-485b-b8f9-96bf4d934917",
+ "metadata": {},
+ "source": [
+ "<h3>6. Simulating a FODO cell with a sextupole:</h3>\n",
+ "\n",
+ "We go back to the stable FODO cell configuration and add a thin sextupole magnet after 1/4 of the lattice, between the first focusing and the second defocusing quadrupole!\n",
+ "\n",
+ "The sextupole kick provides a non-linearity in the potential that confines the particles. At large enough amplitude, the non-linear term dominates the particles are no longer bound / confined!\n",
+ "\n",
+ "We need to track in 4D phase-space (full transverse plane with both x and y), as the sextupole provides coupling terms:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9c298f2d-a26e-4d43-bbc3-7b77613ec6d5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.random.seed(12345)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e44d0e59-8e22-4810-a759-e75d6cd97346",
+ "metadata": {},
+ "source": [
+ "We need a first matrix 1/4 of the cell until the sextupole, one for the $x$ (`M_cell_x_1`) and another one for the $y$ plane (`M_cell_y_1`). Then a second matrix each to track $x$ and $y$ for the remaining 3/4 of the cell (`M_cell_x_2`, `M_cell_y_2`):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0ea158e9-44f9-4a83-b13f-3ecb5ce1cfe8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "k = 0.2\n",
+ "mL = 1.8\n",
+ "\n",
+ "# horizontal plane:\n",
+ "M_cell_x_1 = M_quad_x(0.5, k) # 1/2 focusing quad\n",
+ "M_cell_x_1 = M_cell_x_1.dot(M_drift(2)) # drift\n",
+ "## here sits the sextupole\n",
+ "M_cell_x_2 = M_drift(2) # drift\n",
+ "M_cell_x_2 = M_cell_x_2.dot(M_quad_x(1, -k)) # defocusing quad\n",
+ "M_cell_x_2 = M_cell_x_2.dot(M_drift(4)) # drift\n",
+ "M_cell_x_2 = M_cell_x_2.dot(M_quad_x(0.5, k)) # 1/2 focusing quad\n",
+ "\n",
+ "# vertical plane:\n",
+ "M_cell_y_1 = M_quad_y(0.5, k) # 1/2 focusing quad\n",
+ "M_cell_y_1 = M_cell_y_1.dot(M_drift(2)) # drift\n",
+ "## here sits the sextupole\n",
+ "M_cell_y_2 = M_drift(2) # drift\n",
+ "M_cell_y_2 = M_cell_y_2.dot(M_quad_y(1, -k)) # defocusing quad\n",
+ "M_cell_y_2 = M_cell_y_2.dot(M_drift(4)) # drift\n",
+ "M_cell_y_2 = M_cell_y_2.dot(M_quad_y(0.5, k)) # 1/2 focusing quad"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f06799a9-eb63-4f33-bef5-73c16e870654",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "n_cells = 1000"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "fca54824-a91b-4dc9-bd6a-a6df87b8268b",
+ "metadata": {},
+ "source": [
+ "Initialize the tranverse particle distribution:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e67863be-d15d-4c4d-b853-d63e95a8230f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "N = 100\n",
+ "sig_x = 5e-3\n",
+ "sig_xp = 3e-4\n",
+ "sig_y = 5e-3\n",
+ "sig_yp = 3e-4\n",
+ "\n",
+ "x = np.random.normal(0, sig_x, N)\n",
+ "xp = np.random.normal(0, sig_xp, N)\n",
+ "y = np.random.normal(0, sig_y, N)\n",
+ "yp = np.random.normal(0, sig_yp, N)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ee6e84bb-1358-40f0-8b64-107f124e0f27",
+ "metadata": {},
+ "source": [
+ "Let us record the horizontal phase-space coordinates during the tracking:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1dffadf6-8c46-4578-8582-2b293f80b41d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "rec_x = np.zeros((n_cells, N), dtype=x.dtype)\n",
+ "rec_xp = np.zeros_like(rec_x)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d74ad2b3-829e-4e3d-a9a7-b5d1d4754da5",
+ "metadata": {},
+ "source": [
+ "Let us set the last particle to the same phase-space coordinates as the first particle up to a very small epsilon, for later:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "33f453cd-6b17-4985-bb7c-cad68305870e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "x[-1] = x[0]\n",
+ "xp[-1] = xp[0]\n",
+ "y[-1] = y[0]\n",
+ "yp[-1] = yp[0] * 1.001"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "675484c2-9409-47e8-a7e5-fd7acdee114f",
+ "metadata": {},
+ "source": [
+ "Let's go, the full 4D tracking loop comes here:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3e0c4055-eb74-4881-9b18-3ce99f54f2c7",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i in range(n_cells):\n",
+ " # initial 1/4 of the cell\n",
+ " x, xp = track(M_cell_x_1, x, xp)\n",
+ " y, yp = track(M_cell_y_1, y, yp)\n",
+ " # sextupole\n",
+ " x, xp, y, yp = track_sext_4D(x, xp, y, yp, mL)\n",
+ " # remaining 3/4 of the cell\n",
+ " x, xp = track(M_cell_x_2, x, xp)\n",
+ " y, yp = track(M_cell_y_2, y, yp)\n",
+ " \n",
+ " plt.scatter(np.ones(N) * i, x, c='C0', s=1, marker='.')\n",
+ " plt.scatter([i], [x[-1]], c='r', s=10, marker='.')\n",
+ " \n",
+ " rec_x[i] = x\n",
+ " rec_xp[i] = xp\n",
+ "plt.xlabel('Cells')\n",
+ "plt.ylabel('$x$');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "2230cbcd-e0cf-4009-babe-67ed9cd105d7",
+ "metadata": {},
+ "source": [
+ "How does the horizontal phase-space (Poincaré map) look like?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "da308026-1c30-403d-9282-068952010d93",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i in range(n_cells):\n",
+ " plt.scatter(rec_x[:100,::10], rec_xp[:100,::10], c='C0', s=1, marker='.')\n",
+ "plt.xlabel('$x$')\n",
+ "plt.ylabel(\"$x'$\")\n",
+ "plt.gca().set_aspect(np.diff(plt.xlim()) / np.diff(plt.ylim()));"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "31983881-511a-4281-acf1-e348266ab2e1",
+ "metadata": {},
+ "source": [
+ "$\\implies$ single particles do not maintain the same (linear) ampitude (radius in $x-x'$) during the tracking!\n",
+ "\n",
+ "A single particle looks like this:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "efd05d93-505b-4d57-adbe-4729b014dbce",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(rec_x[:, 0])\n",
+ "plt.xlabel('Turns')\n",
+ "plt.ylabel('$x$');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "04445659-6999-406f-a2b6-e0e859a6ee49",
+ "metadata": {},
+ "source": [
+ "$\\implies$ distorted regular motion (see the asymmetry between positive and negative $x$ values in the oscillation), the particle is still bound but the sextupole deforms the phase-space topology from the regular circles we observed for purely linear tracking.\n",
+ "\n",
+ "Remember, the last particle was just a copy from the first particle with a slightly increased $y'$ value. Let us investigate the difference in their horizontal position during the tracking:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8c9ab0ac-a81d-488a-bd56-387d949db3cc",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(np.abs(rec_x[:,0] - rec_x[:,-1]))\n",
+ "plt.yscale('log')\n",
+ "plt.xlabel('Turns')\n",
+ "plt.ylabel(r'$|\\Delta x|$')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "622c588f-aa4b-4333-a9b3-e2f93aaba518",
+ "metadata": {},
+ "source": [
+ "$\\implies$ for finite sextupole strength, we observe on average an exponential increase. This points to a finite positive maximum Lyapunov exponent, which is an early indicator of deterministic chaos.\n",
+ "\n",
+ "All in all, the thin sextupole magnet in the lattice\n",
+ "- provides a non-linearity in the potential which the particles see\n",
+ "- distorts the regular motion in phase-space\n",
+ "- leads to a change of the (linear) amplitude in phase space\n",
+ "- provides deterministic chaos, in particular at larger amplitudes (positive Lyapunov exponent!)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e10b9110-3c8c-40e2-8087-9676ca171b4f",
+ "metadata": {},
+ "source": [
+ "$\\implies$ repeat this exercise with a zero sextupole strength $m=0$ to confirm these insights for yourself!\n",
+ "\n",
+ "Hint: in order to observe a meaningful result in the last plot, add a factor `* 1.001` to `xp` for the last particle to see an effect. Due to the absent coupling between $x$ and $y$, there will be no difference in the $x$ motion of both particles without the sextupole!"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "86f36f7c-b066-42c2-acc7-7e62db426ec3",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "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
+}