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 [1]:
from IPython.display import Image
fig = Image(filename=('/Users/sterbini/CERNBox/2019/CAS/Vysoke_Tatry/Python/Wednesday/programme_5.png'))
fig
Out[1]:

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 [2]:
# Import the modules
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
import sympy as sy
In [3]:
# Matrix definition
Omega=np.array([[0, 1],[-1,0]])
M=np.array([[1, 0],[1,1]])
M
Out[3]:
array([[1, 0],
       [1, 1]])
In [4]:
# Sum and multiplication of matrices
Omega - M.T @ Omega @ M
# M.T means the "traspose of M".
Out[4]:
array([[0, 0],
       [0, 0]])
In [5]:
aux=np.array([1,2.3,3,])
In [6]:
# Function definition
def Q(f=1):
    return np.array([[1, 0],[-1/f,1]])

np.linalg.det(Q(10))
Out[6]:
1.0
In [7]:
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 [8]:
# 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[8]:
Text(0.5, 1.0, 'My title')
In [9]:
# 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 [10]:
# from Wolfgang's pag.26
fig = Image(filename=('/Users/sterbini/CERNBox/2019/CAS/Vysoke_Tatry/Python/Wednesday/pag26.png'))
fig
Out[10]:

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 [11]:
# 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_1 @ DRIFT_2
Out[11]:
Matrix([
[1, L_1 + L_2],
[0,         1]])

A similar analysis can be done for two thin quadrupoles.

In [12]:
# 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[12]:
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!

SOLUTION

We can still use the symbolic package of python (or do it by hand).

In [13]:
# definition of the focal length symbol
f=sy.Symbol('f');
x=sy.Symbol('x');

# Focusing THIN quadrupole of focal length f
Q=sy.Matrix([[1,0], [-1/f,1]])

# Drift of length f
DRIFT=sy.Matrix([[1,f], [0,1]])

# from Exercise 2, we propagate a particle in parallel to the optical axis
DRIFT @ Q @ sy.Matrix([[x],[0]])
Out[13]:
Matrix([
[   0],
[-x/f]])

and this is on the beam axis (see Exercise 3), whatever is the initial offset, x. This is the meaning of the focal length.

In [14]:
# Remember
fig = Image(filename=('/Users/sterbini/CERNBox/2019/CAS/Vysoke_Tatry/Python/Wednesday/pag14.png'))
fig
Out[14]:

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.$

SOLUTION

In [15]:
# 1D case and thin linear elements
Q=lambda f : np.array([[1, 0],[-1/f,1]])
D=lambda L : np.array([[1, L],[0,1]])
myF=3 # [m]
X=np.array([[1],[0]])
D(myF) @ Q(myF) @ X
Out[15]:
array([[ 0.        ],
       [-0.33333333]])

Indeed we found back the expected results.

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.

SOLUTION

For the first time we introduce beamlines. For representing a beamline can be imaged as a sequece of elements (drift, quadrupoles,...) carrying the linear matrix associated to the element and its length. A possible solution to implement these requirements is to represent the beamline as a list of tuples. Each tuples has two values:

  1. the linear matrix representing the elements
  2. the length of the element.

This approach is quite convenient for two reasons:

  1. we can use the algebra of the list (for python list we can sum two (or more) lists and multiply a list by a scalar). This will ease significantly the beamline composition.
  2. tuples are immutable object, therefore they memory performance is very good.

Now that we have tha data structure to represent a beamline, we can cast a single element a an atomic beamline. So in the following you have the implementation.

In [ ]:
In [16]:
# The drift as a sequence of a single tuple
D = lambda L: [(np.array([[1, L],[0, 1]]), L)]
In [17]:
2*D(1.2)
Out[17]:
[(array([[1. , 1.2],
         [0. , 1. ]]), 1.2), (array([[1. , 1.2],
         [0. , 1. ]]), 1.2)]
In [18]:
# The quadrupole as a sequence of a single tuple
Q = lambda f: [(np.array([[1, 0],[-1/f, 1]]), 0)]

# our parameters
b=3
g=2
f=1/(1/b+1/g)
# VERY IMPORTANT
# This is NOT the SUM of the matrices! It is the concatenation of the list!
beamline=D(b)+Q(f)+D(g)
beamline
Out[18]:
[(array([[1, 3],
         [0, 1]]), 3), (array([[ 1.        ,  0.        ],
         [-0.83333333,  1.        ]]), 0), (array([[1, 2],
         [0, 1]]), 2)]
In [19]:
import pandas as pd
pd.DataFrame(beamline)
Out[19]:
0 1
0 [[1, 3], [0, 1]] 3
1 [[1.0, 0.0], [-0.8333333333333333, 1.0]] 0
2 [[1, 2], [0, 1]] 2

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.

In [20]:
# Refer to 
fig = Image(filename=('/Users/sterbini/CERNBox/2019/CAS/Vysoke_Tatry/Python/Wednesday/pag42.png'))
fig
Out[20]:

SOLUTION

The first part of the exercise is to define the FODO cell starting. A FODO cell is one of the most fundamental lattice topologies. In its simpler version, it consists on two thin quadrupoles of opposite focal length, f, interspaced by two drifts of length L/2 (L is the total length of the cell). We can start the FODO from an arbitrary point inside it. We decided to start describing it from the center of drift. For the moment we can define arbitrarly the f and L paramenters: we chowwse 2.5 m and 1 m respectively.

In [21]:
f=2.5
L_2=1
beamline=5*D(L_2/2/5)+Q(f)+10*D(L_2/(10))+Q(-f)+5*D(L_2/2/5)

pd.DataFrame(beamline)
Out[21]:
0 1
0 [[1.0, 0.1], [0.0, 1.0]] 0.1
1 [[1.0, 0.1], [0.0, 1.0]] 0.1
2 [[1.0, 0.1], [0.0, 1.0]] 0.1
3 [[1.0, 0.1], [0.0, 1.0]] 0.1
4 [[1.0, 0.1], [0.0, 1.0]] 0.1
5 [[1.0, 0.0], [-0.4, 1.0]] 0.0
6 [[1.0, 0.1], [0.0, 1.0]] 0.1
7 [[1.0, 0.1], [0.0, 1.0]] 0.1
8 [[1.0, 0.1], [0.0, 1.0]] 0.1
9 [[1.0, 0.1], [0.0, 1.0]] 0.1
10 [[1.0, 0.1], [0.0, 1.0]] 0.1
11 [[1.0, 0.1], [0.0, 1.0]] 0.1
12 [[1.0, 0.1], [0.0, 1.0]] 0.1
13 [[1.0, 0.1], [0.0, 1.0]] 0.1
14 [[1.0, 0.1], [0.0, 1.0]] 0.1
15 [[1.0, 0.1], [0.0, 1.0]] 0.1
16 [[1.0, 0.0], [0.4, 1.0]] 0.0
17 [[1.0, 0.1], [0.0, 1.0]] 0.1
18 [[1.0, 0.1], [0.0, 1.0]] 0.1
19 [[1.0, 0.1], [0.0, 1.0]] 0.1
20 [[1.0, 0.1], [0.0, 1.0]] 0.1
21 [[1.0, 0.1], [0.0, 1.0]] 0.1

Please pay attention that in order to sample more finely the particle trajectory in the drift we split them in pieces of 0.2 m. Now we need to solve the second part of the exercise, that is the actual tracking along the beamline we just defined. To record the position of the particle(s) around the sequence we will use a data structure identical to the one used for the beamline: a list of tuples. In this case, however, the meaning of the item in the tuples is different. Each tuples has two values:

  1. the position in the trace space of the particle(s).
  2. the corresponding s-position in the beamline.

We call it the beam-list variable.

We will see in the following cell how, starting from the initial position of the particle, we can trasport it along the beam line.

In [22]:
X=np.array([[0],[0.001]])
# the starting beamlist contains only a position (s=0)
beamlist=[(X,0)]
beamlist
Out[22]:
[(array([[0.   ],
         [0.001]]), 0)]

The proper tracking will append to the initial beamlist the transformed coordinate in the trace-space and the incremented s-coordinate.

In [23]:
for element in beamline:
    beamlist.append((element[0] @ beamlist[-1][0],beamlist[-1][1]+ element[1]))
In [24]:
import pandas as pd
pd.DataFrame(beamlist, columns=['coordinates', 's-position'])
Out[24]:
coordinates s-position
0 [[0.0], [0.001]] 0.0
1 [[0.0001], [0.001]] 0.1
2 [[0.0002], [0.001]] 0.2
3 [[0.00030000000000000003], [0.001]] 0.3
4 [[0.0004], [0.001]] 0.4
5 [[0.0005], [0.001]] 0.5
6 [[0.0005], [0.0008]] 0.5
7 [[0.00058], [0.0008]] 0.6
8 [[0.00066], [0.0008]] 0.7
9 [[0.00074], [0.0008]] 0.8
10 [[0.00082], [0.0008]] 0.9
11 [[0.0009], [0.0008]] 1.0
12 [[0.00098], [0.0008]] 1.1
13 [[0.00106], [0.0008]] 1.2
14 [[0.00114], [0.0008]] 1.3
15 [[0.00122], [0.0008]] 1.4
16 [[0.0013], [0.0008]] 1.5
17 [[0.0013], [0.00132]] 1.5
18 [[0.001432], [0.00132]] 1.6
19 [[0.0015639999999999999], [0.00132]] 1.7
20 [[0.0016959999999999998], [0.00132]] 1.8
21 [[0.0018279999999999998], [0.00132]] 1.9
22 [[0.00196], [0.00132]] 2.0
In [25]:
# just to show the first 5 rows
pd.DataFrame(beamlist, columns=['coordinates', 's-position']).head()
Out[25]:
coordinates s-position
0 [[0.0], [0.001]] 0.0
1 [[0.0001], [0.001]] 0.1
2 [[0.0002], [0.001]] 0.2
3 [[0.00030000000000000003], [0.001]] 0.3
4 [[0.0004], [0.001]] 0.4
In [26]:
from matplotlib import pyplot as plt
%matplotlib inline
params = {'legend.fontsize': 'x-large',
          'figure.figsize': (15, 5),
         'axes.labelsize': 'x-large',
         'axes.titlesize':'x-large',
         'xtick.labelsize':'x-large',
         'ytick.labelsize':'x-large'}
plt.rcParams.update(params)
plt.plot([i[1] for i in beamlist],[i[0][0]*1000 for i in beamlist],'o-b',lw=3)
plt.grid(True)
plt.xlabel('s [m]')
plt.ylabel('x [mm]')
plt.title('Exercise 7');

Exercise 8

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

SOLUTION

In [27]:
plt.plot([i[1] for i in beamlist],[i[0][1]*1000 for i in beamlist],'o-r',lw=3)
plt.grid(True)
plt.xlabel('s [m]')
plt.ylabel("x' [mrad]");
plt.title('Exercise 8');

Exercise 9

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

SOLUTION

In [28]:
beamline_5FODO=5*beamline
X=np.array([[0],[0.001]])
# the starting beamlist contains only a position (s=0)
beamlist=[(X,0)]
for element in beamline_5FODO:
    beamlist.append((element[0] @ beamlist[-1][0],beamlist[-1][1]+ element[1]))
    
plt.plot([i[1] for i in beamlist],[i[0][0]*1000 for i in beamlist],'o-b',lw=3)
plt.grid(True)
plt.xlabel('s [m]')
plt.gca().set_ylabel('x [mm]', color='b')
plt.gca().tick_params(axis='y', labelcolor='b')

ax2 = plt.gca().twinx()  # instantiate a second axes that shares the same x-axis
ax2.set_ylabel("x' [mrad]", color='r')
ax2.tick_params(axis='y', labelcolor='r')
ax2.plot([i[1] for i in beamlist],[i[0][1]*1000 for i in beamlist],'o-r',lw=3)
plt.title('Exercise 9');

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.

In [29]:
f=.4000
L_2=1
fodo_lattice= 5*D(L_2/2/5)+Q(f)+10*D(L_2/(10))+Q(-f)+5*D(L_2/2/5)
beamline=100*fodo_lattice
X=np.array([[0],[0.001]])
# the starting beamlist contains only a position (s=0)
beamlist=[(X,0)]
for element in beamline:
    beamlist.append((element[0] @ beamlist[-1][0],beamlist[-1][1]+ element[1]))
    
plt.plot([i[1] for i in beamlist],[i[0][0]*1000 for i in beamlist],'o-b',lw=3)
plt.grid(True)
plt.xlabel('s [m]')
plt.ylabel('x [mm]')

ax2 = plt.gca().twinx()  # instantiate a second axes that shares the same x-axis
ax2.set_ylabel("x' [mrad]", color='r')
ax2.tick_params(axis='y', labelcolor='r')
ax2.plot([i[1] for i in beamlist],[i[0][1]*1000 for i in beamlist],'o-r',lw=3)

plt.title('Exercise 10');

It it important to observe that when the focal length approaches the stability limit (only for $f>L/4$ the machine is stable), then the oscillation start to grow during the transport in the lattice.

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.

In [30]:
# a bit of parameters
b=3
L=5
g=L-b
f=1/(1/b+1/g)

# the usual approach

beamline=10*D(b/10)+Q(f)+10*D(g/10)

X=np.array([[0],[0.001]])
# the starting beamlist contains only a position (s=0)
beamlist=[(X,0)]
for element in beamline:
    beamlist.append((element[0] @ beamlist[-1][0],beamlist[-1][1]+ element[1]))
    
plt.plot([i[1] for i in beamlist],[i[0][0]*1000 for i in beamlist],'o-b',lw=3)
plt.grid(True)
plt.xlabel('s [m]')
plt.ylabel('x [mm]')

ax2 = plt.gca().twinx()  # instantiate a second axes that shares the same x-axis
ax2.set_ylabel("x' [mrad]", color='r') 
ax2.tick_params(axis='y', labelcolor='r')
ax2.plot([i[1] for i in beamlist],[i[0][1]*1000 for i in beamlist],'o-r',lw=3)

plt.title('Exercise 11');

We can do a bit animation showing that, moving the position of the quadruple and keeping fulfilled the imaging equation the particle entering with x=0 will exit with x=0.

In [31]:
# jupyter labextension install @jupyter-widgets/jupyterlab-manager
from ipywidgets import interactive

def plotIt(b):
    L=5
    g=L-b
    f=1/(1/b+1/g)

    # the usual approach

    beamline=10*D(b/10)+Q(f)+10*D(g/10)

    X=np.array([[0],[0.001]])
    # the starting beamlist contains only a position (s=0)
    beamlist=[(X,0)]
    for element in beamline:
        beamlist.append((element[0] @ beamlist[-1][0],beamlist[-1][1]+ element[1]))

    plt.plot([i[1] for i in beamlist],[i[0][0]*1000 for i in beamlist],'o-b',lw=3)
    plt.grid(True)
    plt.xlabel('s [m]')
    plt.ylabel('x [mm]')
    plt.ylim([0,5])


    ax2 = plt.gca().twinx()  # instantiate a second axes that shares the same x-axis
    ax2.set_ylabel("x' [mrad]", color='r') 
    ax2.tick_params(axis='y', labelcolor='r')
    ax2.plot([i[1] for i in beamlist],[i[0][1]*1000 for i in beamlist],'o-r',lw=3)

    plt.title('Exercise 11, animated')
    plt.ylim([-5,5])
interactive_plot = interactive(plotIt,b=(2,4,.1),continuous_update=False)
output = interactive_plot.children[-1]
output.layout.height = '350px'
interactive_plot

Exercise 12

Define an ensemble of 1000 particle 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.

In [32]:
# Refer to 
fig = Image(filename=('/Users/sterbini/CERNBox/2019/CAS/Vysoke_Tatry/Python/Wednesday/pag44.png'))
fig
Out[32]:

SOLUTION

In [33]:
# simple solution
Npart=1000
beam=np.random.randn(2, Npart)
x0=0;
xp0=1
sigx=1; 
sigxp=0.5; 
beam[0,:]=sigx*beam[0,:]+x0
beam[1,:]=sigxp*beam[1,:]+xp0
np.std(beam[1,:])
Out[33]:
0.5134883501995198

We did not find exactly back the 0.5 value due to the statistical error coming with the finite numbers of particles. We can explore the dependence of the statistical error as function of the number of particles by a numerical simulation.

In [34]:
myDivergence=[]
myDivergenceStatisticalError=[]
Npart_range=np.round(np.logspace(2,6,10))
Npart_range=Npart_range.astype(int)
for Npart in Npart_range:
    aux=[]
    for j in range(10):
        beam=np.random.randn(2, Npart)
        x0=0;
        xp0=1
        sigx=1; 
        sigxp=0.5; 
        beam[0,:]=sigx*beam[0,:]+x0
        beam[1,:]=sigxp*beam[1,:]+xp0
        aux.append(np.std(beam[1,:]))
    myDivergence.append(np.mean(aux))
    myDivergenceStatisticalError.append(np.std(aux))
plt.errorbar(np.log10(Npart_range),myDivergence,myDivergenceStatisticalError,lw=5)
plt.plot([2,6],[.5,.5],':g',lw=5)
plt.grid(True)
plt.xlabel('log10(N_part)')
plt.ylabel('Beam divergence of the ensemble');

Depending on the required precision we need to select the convenient number of particles.

Exercise 13

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

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 triangle of surface A. Verify that this linear trasport preserve the area of the triangle.

In [35]:
# simple solution
Npart=10000
beam=np.random.randn(2, Npart)
x0=0;
xp0=1
sigx=1; 
sigxp=0.5; 
beam[0,:]=sigx*beam[0,:]+x0
beam[1,:]=sigxp*beam[1,:]+xp0
In [36]:
beamAfterDrift=D(1)[0][0]@beam
In [37]:
# Before the drift
import seaborn as sns
g = sns.jointplot(x=beam[0,:], y=beam[1,:], kind="hex", color="b")
g.set_axis_labels("x [mm]", "x' [mrad]");
In [38]:
# After the drift
import seaborn as sns
g = sns.jointplot(x=beamAfterDrift[0,:], y=beamAfterDrift[1,:], kind="hex", color="b")
g.set_axis_labels("x [mm]", "x' [mrad]");
In [39]:
# just to have fun
from ipywidgets import interactive
import seaborn as sns

D = lambda L: [(np.array([[1, L],[0, 1]]), L)]

Npart=10000
beam=np.random.randn(2, Npart)
x0=0;
xp0=1
sigx=1; 
sigxp=0.5;  
beam[0,:]=sigx*beam[0,:]+x0
beam[1,:]=sigxp*beam[1,:]+xp0
beamAfterDrift=D(10)[0][0]@beam

def plotIt(drift_length):
    beamAfterDrift=D(drift_length)[0][0]@beam
    g = sns.jointplot(x=beamAfterDrift[0,:], y=beamAfterDrift[1,:], kind="hex", color="b")
    g.set_axis_labels("x [mm]", "x' [mrad]");
    g.ax_marg_x.set_xlim(-20, 30)

interactive_plot = interactive(plotIt,drift_length=(0,10,1),continuous_update=False)
output = interactive_plot.children[-1]
output.layout.height = '500px'
interactive_plot
In [40]:
# test of linearity
month=2
v1=D(1)[0][0]@(17*month*beam)
v2=17*month*D(1)[0][0]@(beam)
np.array_equal(v1,v2)
Out[40]:
True

Area preservation

Now let us consider three points of the distributions and let us verify that the area of the triangle before and after the transformation is preserved.

In [41]:
# Remember Suzie's lecture
fig = Image(filename=('/Users/sterbini/CERNBox/2019/CAS/Vysoke_Tatry/Python/Wednesday/susie.png'))
fig
Out[41]:
In [42]:
# this is the formula of the area of the triangle given the three coordiantes of the vertices
from numpy.linalg import norm
def area(a, b, c) :
    return 0.5 * norm( np.cross( b-a, c-a))
In [43]:
# we take the first three particles in the beam before the trasport (arbitrary choice)
area(beam[:,0],beam[:,1],beam[:,2])
Out[43]:
0.570209999854786
In [44]:
# after transport
area(beamAfterDrift[:,0],beamAfterDrift[:,1],beamAfterDrift[:,2])
Out[44]:
0.570209999854786
In [45]:
#first triangle
plt.plot(beam[:,:3][0],beam[:,:3][1],'o-b')
plt.plot(beam[:,[2,0]][0],beam[:,[2,0]][1],'o-b', label='before the drift')

#second triangle
plt.plot(beamAfterDrift[:,:3][0],beamAfterDrift[:,:3][1],'o-r')
plt.plot(beamAfterDrift[:,[2,0]][0],beamAfterDrift[:,[2,0]][1],'o-r',label='after the drift')

plt.legend()
plt.xlabel('x [mm]')
plt.ylabel("x' [mrad]")
plt.grid(True)
In [46]:
# just to have fun
from ipywidgets import interactive


def plotIt(i):
    #first triagle
    plt.plot(beam[0,:],beam[1,:],'.b',alpha=.01)
    plt.plot(beamAfterDrift[0,:],beamAfterDrift[1,:],'.r',alpha=.01)

    mySmallBeam=beam[:,(0+i):(3+i)];
    plt.plot(mySmallBeam[:,:3][0],mySmallBeam[:,:3][1],'o-b',lw=3)
    plt.plot(mySmallBeam[:,[2,0]][0],mySmallBeam[:,[2,0]][1],'o-b',lw=3, label='before the drift')

    #second triagle
    mySmallBeamAfterDrift=beamAfterDrift[:,(0+i):(3+i)];

    plt.plot(mySmallBeamAfterDrift[:,:3][0],mySmallBeamAfterDrift[:,:3][1],'o-r',lw=3)
    plt.plot(mySmallBeamAfterDrift[:,[2,0]][0],mySmallBeamAfterDrift[:,[2,0]][1],'o-r',lw=3,label='after the drift')
    plt.xlim(-5,25)
    plt.ylim(-.5,3)
    plt.legend(loc=1)
    plt.xlabel('x [mm]')
    plt.ylabel("x' [mrad]")
    plt.grid(True)

interactive_plot = interactive(plotIt,i=(0,100,1),continuous_update=False)
output = interactive_plot.children[-1]
output.layout.height = '500px'
interactive_plot

Now remember the first few instuctions we ran today:

In [47]:
# Matrix definition
Omega=np.array([[0, 1],[-1,0]])
M=np.array([[1, 1],[0,1]])

# Sum and multiplication of matrices
Omega - M.T @ Omega @ M
# M.T means the "traspose of M".
# This is to verify that the matrix is symplectic.
Out[47]:
array([[0, 0],
       [0, 0]])

Covariance matrix of the beam distribution

In [48]:
# very important to set the bias=True
np.cov(beam,bias=True)
Out[48]:
array([[ 0.98402864, -0.00156091],
       [-0.00156091,  0.24846869]])
In [49]:
# The meaning of the covariance matrix or sigma matrix is the following (see Primer)
x_mean=np.mean(beam[0,:])
xp_mean=np.mean(beam[1,:])
sigma11= np.mean([(i-x_mean)**2 for i in beam[0]])
sigma22= np.mean([(i-xp_mean)**2 for i in beam[1]])
sigma12=sigma21= np.mean([(i-x_mean)*(j-xp_mean) for i,j in zip(beam[0],beam[1])])
print(f'''sigma11 is {sigma11}
sigma22 is {sigma22}
sigma12 is {sigma12}
''')

sigma11 is 0.9840286443650799
sigma22 is 0.2484686892734233
sigma12 is -0.001560910923525561
In [50]:
# The rms emittance is
print(f'rms emittance before: {np.sqrt(np.linalg.det(np.cov(beam,bias=True)))} mm mrad')
print(f'rms emittance after : {np.sqrt(np.linalg.det(np.cov(beamAfterDrift,bias=True)))} mm mrad')

rms emittance before: 0.4944672598160407 mm mrad
rms emittance after : 0.4944672598160354 mm mrad
In [51]:
# What is the impact of the numerical error if we consider very long drifts...numerical stability!!!
myRange=np.logspace(0,7.8,500)
plt.semilogx(myRange,[np.sqrt(np.linalg.det(np.cov(D(i)[0][0]@beam,bias=True))) for i in myRange],'o-b')
plt.grid(True)
plt.xlabel('drift length [m]')
plt.ylabel('rms beam emittance [mm mrad]');

/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:3: RuntimeWarning: invalid value encountered in sqrt
  This is separate from the ipykernel package so we can avoid doing imports until