# Exercise 14¶

Using the 5 FODO cells of Exercise 9, trasport the beam of Exercise 13 and plot its rms size and divergence along the line.

### SOLUTION¶

In [1]:
# Import the modules
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
import sympy as sy
from matplotlib import pyplot as plt
%matplotlib inline

D = lambda L: [(np.array([[1, L],[0, 1]]), L)]
Q = lambda f: [(np.array([[1, 0],[-1/f, 1]]), 0)]

In [2]:
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)

#prepare the optics
f=2.5
L_2=2
fodo_lattice= 5*D(L_2/2/5)+Q(-f)+10*D(L_2/(10))+Q(f)+5*D(L_2/2/5)
beamline_5FODO=5*fodo_lattice

#prepare the beam
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

beamlist=[(beam,0)]
for element in beamline_5FODO:
beamlist.append((element[0] @ beamlist[-1][0],beamlist[-1][1]+ element[1]))

plt.plot([x[1] for x in beamlist],[np.mean(x[0][0,:])for x 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.tick_params(axis='y', labelcolor='r')
ax2.plot([x[1] for x in beamlist],[np.std(x[0][1,:])for x in beamlist],'o-r',lw=3)
plt.title('Exercise 14');

In [3]:
plt.plot([x[1] for x in beamlist],[np.std(x[0][0,:])for x in beamlist],'-b',lw=3)
plt.grid(True)
plt.xlabel('s [m]')
plt.gca().set_ylabel('beam size [mm]', color='b')
plt.gca().tick_params(axis='y', labelcolor='b')


# Exercise 15¶

Starting from Exercise 14, what happens if you (a) increase or (b) reduce by a factor 2 the initial beam divercence (sigxp)?

### SOLUTIONS¶

In [4]:
# CASE a
f=2.5
L_2=2
fodo_lattice= 5*D(L_2/2/5)+Q(f)+10*D(L_2/(10))+Q(-f)+5*D(L_2/2/5)
beamline_5FODO=5*fodo_lattice

#prepare the beam
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

beamlist=[(beam,0)]
for element in beamline_5FODO:
beamlist.append((element[0] @ beamlist[-1][0],beamlist[-1][1]+ element[1]))

plt.plot([x[1] for x in beamlist],[np.mean(x[0][0,:])for x 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.tick_params(axis='y', labelcolor='r')
ax2.plot([x[1] for x in beamlist],[np.std(x[0][1,:])for x in beamlist],'o-r',lw=3)
plt.ylim([0,2])

plt.title('Exercise 15, case A')

# CASE b, just a copy of the first one with a small change...
f=2.5
L_2=2
fodo_lattice= 5*D(L_2/2/5)+Q(f)+10*D(L_2/(10))+Q(-f)+5*D(L_2/2/5)
beamline_5FODO=5*fodo_lattice

#prepare the beam
Npart=10000
beam=np.random.randn(2, Npart)
x0=0;
xp0=1
sigx=1;
sigxp=0.5/2;
beam[0,:]=sigx*beam[0,:]+x0
beam[1,:]=sigxp*beam[1,:]+xp0

beamlist=[(beam,0)]
for element in beamline_5FODO:
beamlist.append((element[0] @ beamlist[-1][0],beamlist[-1][1]+ element[1]))

plt.figure()
plt.plot([x[1] for x in beamlist],[np.mean(x[0][0,:])for x 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.tick_params(axis='y', labelcolor='r')
ax2.plot([x[1] for x in beamlist],[np.std(x[0][1,:])for x in beamlist],'o-r',lw=3)
plt.ylim([0,2])

plt.title('Exercise 15, case B')

Out[4]:
Text(0.5, 1.0, 'Exercise 15, case B')

It very important to note that the divergence of CASE A and CASE B are not proportional.

# Exercise 16¶

Using Equation from the Primer, display (a) the average position of the particles along the beam line. Likewise, (b) display the angular divergence. Compare with the result you found in Exercise 14.

### SOLUTIONS¶

In [5]:
#lattice
f=2.5
L_2=2
fodo_lattice= 5*D(L_2/2/5)+Q(f)+10*D(L_2/(10))+Q(-f)+5*D(L_2/2/5)
beamline_5FODO=5*fodo_lattice

#prepare the beam
Npart=10000
beam0=np.random.randn(2, Npart)
x0=0;
xp0=1
sigx=1;
sigxp=0.5;
beam0[0,:]=sigx*beam0[0,:]+x0
beam0[1,:]=sigxp*beam0[1,:]+xp0
beamlist=[(beam0,0)]

#prepare the sigma matrix
sigma0=np.array([[sigx**2, 0],[0, sigxp**2]])
sigmalist=[(sigma0,0)]

for element in beamline_5FODO:
beamlist.append((element[0] @ beamlist[-1][0],beamlist[-1][1]+ element[1]))

for element in beamline_5FODO:
sigmalist.append((element[0] @ sigmalist[-1][0] @ element[0].transpose() ,sigmalist[-1][1]+ element[1]))

plt.plot([x[1] for x in beamlist],[np.std(x[0][0,:])for x in beamlist],'-b',lw=3, label='particle tracking')
plt.plot([x[1] for x in sigmalist],[np.sqrt(x[0][0,0]) for x in sigmalist],'or', lw=3,label='sigma tracking')
plt.grid(True)
plt.legend(loc='best')
plt.xlabel('s [m]')
plt.ylabel('beam size [mm]');


# Exercise 17¶

Can you find an initial beam matrix sigma0 that reproduces itself at the end of the beam line?

### SOLUTION¶

This problem is not so simple. One could proceed with try and error. But there are three parameters to fix so it can be quite cumbersome.

In [43]:
from ipywidgets import interactive

def plotIt(sigma11, sigma22, sigma12):
#lattice
f=2.5
L_2=2
fodo_lattice= 5*D(L_2/2/5)+Q(f)+10*D(L_2/(10))+Q(-f)+5*D(L_2/2/5)
beamline_5FODO=5*fodo_lattice

#prepare the sigma matrix
sigma0=np.array([[sigma11,sigma12 ],[sigma12, sigma22]])
sigmalist=[(sigma0,0)]

for element in beamline_5FODO:
sigmalist.append((element[0] @ sigmalist[-1][0] @ element[0].transpose() ,sigmalist[-1][1]+ element[1]))

plt.plot([x[1] for x in sigmalist],[np.sqrt(x[0][0,0]) for x in sigmalist],'-b',lw=3)
plt.grid(True)

interactive_plot = interactive(plotIt,sigma11=(0,10,.1),sigma22=(0,10,.01),sigma12=(-10,10,.01),continuous_update=False)
output = interactive_plot.children[-1]
output.layout.height = '350px'
interactive_plot

In [44]:
def plotIt(sigma11):
#lattice
f=2.5
L_2=2
sigma22=1/sigma11
sigma12=0
# Very convenient to chose a point with no x-xp correlation
fodo_lattice= Q(2*f)+10*D(L_2/10)+Q(-f)+10*D(L_2/(10))+Q(2*f)
beamline_5FODO=5*fodo_lattice

#prepare the sigma matrix
sigma0=np.array([[sigma11,sigma12 ],[sigma12, sigma22]])
sigmalist=[(sigma0,0)]

for element in beamline_5FODO:
sigmalist.append((element[0] @ sigmalist[-1][0] @ element[0].transpose() ,sigmalist[-1][1]+ element[1]))

#plt.plot([x[1] for x in beamlist],[np.std(x[0][0,:])for x in beamlist],'-b',lw=3)
plt.plot([x[1] for x in sigmalist],[np.sqrt(x[0][0,0]) for x in sigmalist],'b-',lw=3)
plt.grid(True)

interactive_plot = interactive(plotIt,sigma11=(0,10,.1),continuous_update=False)
output = interactive_plot.children[-1]
output.layout.height = '350px'
interactive_plot


If the lattice is stable (and this is the case of our simple FODO, $f>L_{FODO}/4$), then only a specific sigma matrix (after having fixed the beam emittance, in our example we choose unitary emittance) shows a periodic behaviour. The beam associated to that sigma matrix is the matched beam.

# Periodical systems¶

Now we are going to explore the behaviour of the single particle (and later of the beam), not along the generic position s but turn-after-turn, i.e., $s_0,~s_0+L_{period},~s_0+2L_{period},...$ where $L_{period}$ is the period of our periodic lattice.

In [8]:
# Refer to
from IPython.display import Image
fig = Image(filename=('/Users/sterbini/CERNBox/2019/CAS/Vysoke_Tatry/Python/Friday/pag52.png'))
fig

Out[8]: