An overview of the Hands-On lattice calculations

We have 7 h scheduled this week (https://cas.web.cern.ch/sites/cas.web.cern.ch/files/programmes/vysoke-tatry-2019-programme_5.pdf).

In [15]:
from IPython.display import Image
fig = Image(filename=('/Users/sterbini/CERNBox/2019/CAS/Vysoke_Tatry/Python/Wednesday/programme_5.png'))
fig
Out[15]:

The goal is to complement with hand-on numerical exercises the lectures of Transverse Linear Beam Dynamics:

  • helping you to visualize the concepts,
  • introducing you to the basic numerical approach.

On Wednesday 11.09 we will cover the following topics

  • setup of the enviroment and familiarization with python
  • first tracking program (1D): definition of the lattice elements, lattice list and beam list
  • from a single particle to a beam

On Friday 13.09 we will cover

  • tracking of second order moments of the particle ensemble
  • periodic solution of the lattice
  • matching ("try and error")

On Saturday 14.09 we will generalize the code including

  • off momentum (dispersion function, 3x3 matrices)
  • 2D motion (horizontal AND vertical plane, 4x4 matrices)
  • an example of coupling

We will fine tune the pace depending on the needs.

The format is quite simple: we are going to solve the exercises presented in https://indico.cern.ch/event/808940/contributions/3553546/attachments/1904762/3145489/CAS_Optics_Primer.pdf

You will try first on your own for few minutes and after we will solve and comment together.

The installation procedure for Python3 (it will be used also next week for the longitudinal hands-on) can be found at https://codimd.web.cern.ch/s/HkfQy3YbB

Do not hesitate to ask Werner, Volker, Andrea or Guido if you have questions, doubts or problems!

An ice-breaker example

Just to get familiar with a bit of syntax

In [89]:
# Import the modules
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
import sympy as sy
In [67]:
# Matrix definition
Omega=np.array([[0, 1],[-1,0]])
M=np.array([[1, 0],[1,1]])
M
Out[67]:
array([[1, 0],
       [1, 1]])
In [68]:
# Sum and multiplication of matrices
Omega - M.T @ Omega @ M
# M.T means the "traspose of M".
Out[68]:
array([[0, 0],
       [0, 0]])
In [71]:
# Function definition
def Q(f=1):
    return np.array([[1, 0],[-1/f,1]])

np.linalg.det(Q(10))
Out[71]:
1.0
In [77]:
np.linalg.det?

Signature: np.linalg.det(a)
Docstring:
Compute the determinant of an array.

Parameters
----------
a : (..., M, M) array_like
    Input array to compute determinants for.

Returns
-------
det : (...) array_like
    Determinant of `a`.

See Also
--------
slogdet : Another way to represent the determinant, more suitable
  for large matrices where underflow/overflow may occur.

Notes
-----

.. versionadded:: 1.8.0

Broadcasting rules apply, see the `numpy.linalg` documentation for
details.

The determinant is computed via LU factorization using the LAPACK
routine z/dgetrf.

Examples
--------
The determinant of a 2-D array [[a, b], [c, d]] is ad - bc:

>>> a = np.array([[1, 2], [3, 4]])
>>> np.linalg.det(a)
-2.0

Computing determinants for a stack of matrices:

>>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ])
>>> a.shape
(3, 2, 2)
>>> np.linalg.det(a)
array([-2., -3., -8.])
File:      /anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py
Type:      function
In [75]:
# check a simple plot
plt.plot([0,10],[0,10],'ob-')
plt.xlabel('My x-label [arb. units]')
plt.ylabel('My y-label [arb. units]')
plt.title('My title')
Out[75]:
Text(0.5, 1.0, 'My title')
In [78]:
# an sns plot
sns.set(style="ticks")
rs = np.random.RandomState(10)
x = rs.normal(size=10000)
y = rs.normal(size=10000)
sns.jointplot(x, y, kind="hex");

Linear Optics and matrices

In [19]:
# from Wolfgang's pag.26
fig = Image(filename=('/Users/sterbini/CERNBox/2019/CAS/Vysoke_Tatry/Python/Wednesday/pag26.png'))
fig
Out[19]:

Exercise 1

Show that multiplying two drift matrices, one with $L_1$ and the other with $L_2$ in the upper right corner, produces a matrix with the sum of the distances in the upper right corner.

SOLUTION

Even if we can prove it by hand, we take the opportunity to show some symbolic computation in python.

In [21]:
# we introduce the symbols
L_1=sy.Symbol('L_1');
L_2=sy.Symbol('L_2');

# we define the matrices
DRIFT_1=sy.Matrix([[1,L_1],[0,1]])
DRIFT_2=sy.Matrix([[1,L_2],[0,1]])

# we multiply the matrices. 
# NOTA BENE: the @ operator is the "multiplication between matrix" in py3
DRIFT_2 @ DRIFT_1
Out[21]:
Matrix([
[1, L_1 + L_2],
[0,         1]])

A similar analysis can be done for two thin quadrupoles.

In [85]:
# we introduce the symbols
f_1=sy.Symbol('f_1');
f_2=sy.Symbol('f_2');

# we define the matrices
Q_1=sy.Matrix([[1,0],[1/f_1,1]])
Q_2=sy.Matrix([[1,0],[1/f_2,1]])

# we multiply the matrices. 
# NOTA BENE: the @ operator is the "multiplication between matrix" in py3
Q_2 @ Q_1
Out[85]:
Matrix([
[            1, 0],
[1/f_2 + 1/f_1, 1]])

Exercise 2

How do you describe a ray that is parallel to the optical axis?

SOLUTION

The phase space vector of a particle parallel ot the optical axis can be represented (1D) by \begin{equation} X=\left( \begin{array}{c} x \\ 0 \end{array} \right) \end{equation}

Exercise 3

How do you describe a ray that is on the optical axis?

SOLUTION

The phase space vector of a particle sitting on the optical axis can be represented by \begin{equation} X=\left( \begin{array}{c} 0 \\ x' \end{array} \right). \end{equation} We assumed that "on the optical axis" referred only to the position and not to the particle angle, x'.

Exercise 4

Show by multiplying the respective matrices that a parallel ray, which first passes through a lens with focal length $f$ and then moves on a straight line, actually crosses the optical axis at a distance $L=f$ downstream of the lens. Hint: think a little extra about ordering of the matrices!

Exercise 5

Set f=3 and numerically verify what you found in Exercise 4, namely that parallel rays cross the axis after a distance $L=f.$

Exercise 6

Recall that the imaging equation for a lens is $1/b+1/g=1/f,$ which corresponds to a system of one focusing lens with focal length $f,$ sandwiched between drift spaces with length $g$ and $b,$ respectively. Write a beamline description that corresponds to this system. We will later return to it and analyze it.

Exercise 7

Define a FODO beamline and prepare initial coordinates that describe a particle that is on the optical axis, but has an initial angle $x'$ and plot the position $x$ along the beam line.

Exercise 8

Plot the angle $x'$ along the beam line.

Exercise 9

Plot both the position $x$ and the angle $x'$ through five cells.

Exercise 10

Plot the position $x$ through 100 cells, play with different values of the focal length F and explore whether you can make the oscillations grow.

Exercise 11

Use the beam line for the imaging system you prepared in Exercise 6 and launch a particle with $x_0=0$ and an angle of $x'_0=1\,$mrad at one end. Verify that this particle crosses the center of the beam pipe at the exit of the beam line, provided that $b,g,$ and $f$ satisfy the imaging equation that is shown in Exercise 6.

Exercise 12

Define an ensemble of 1000 particles with an arbitrary first order ($x0$, $xpo$) and second order momenta( $x_{rms}$ and $x'_{rms}$) Verify the angular divergence of the beam is the one set.

Exercise 13

Trasport the beam distribution of Exercise 12 in a drift of length 1 m. Compare the initial and final distribution.

Test of linearity. Scale the input vector by 17 times the month of your birthday (85 if you are born in May) and verify that the output vector from the matrix multiplication has changed by the same factor.

Now launch 3 particles such that they define a trianle of surface A. Verify that this linear trasport preserve the area of the tiangle.