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.set_ylabel("divergence [mrad]", color='r')
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.set_ylabel("divergence [mrad]", color='r')
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.set_ylabel("divergence [mrad]", color='r')
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]: