Skip to content
Snippets Groups Projects
Unverified Commit 1ec966c7 authored by tegenolf's avatar tegenolf Committed by GitHub
Browse files

Add notebook to week 7

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