Quadrupole Example¶
Simple quadrupole example
In [1]:
Copied!
from impact import Impact
from pmd_beamphysics.units import mec2
import numpy as np
import os
import matplotlib.pyplot as plt
%config InlineBackend.figure_format='retina'
from impact import Impact
from pmd_beamphysics.units import mec2
import numpy as np
import os
import matplotlib.pyplot as plt
%config InlineBackend.figure_format='retina'
In [2]:
Copied!
# locate the drift template
ifile = '../templates/quadrupole/ImpactT.in'
os.path.exists(ifile)
# locate the drift template
ifile = '../templates/quadrupole/ImpactT.in'
os.path.exists(ifile)
Out[2]:
True
In [3]:
Copied!
# calculate gamma*beta
Etot = 6e6 #eV
gamma = Etot/mec2
GB = np.sqrt(gamma**2 -1)
GB
# calculate gamma*beta
Etot = 6e6 #eV
gamma = Etot/mec2
GB = np.sqrt(gamma**2 -1)
GB
Out[3]:
11.699046356605859
Use Impact's built-in Gaussian particle generator¶
In [4]:
Copied!
I = Impact(ifile)
I.header['Np'] = 100000
I.header['Nx'] = 32
I.header['Ny'] = 32
I.header['Nz'] = 32
I.header['Dt'] = 10e-12
I.header['Bcurr'] = 0
I.header['zmu2'] = GB
# set normal and skew quads
I.ele['CQ01']['b1_gradient'] = 0.00714 # T/m
I.ele['SQ01']['b1_gradient'] = 0
I = Impact(ifile)
I.header['Np'] = 100000
I.header['Nx'] = 32
I.header['Ny'] = 32
I.header['Nz'] = 32
I.header['Dt'] = 10e-12
I.header['Bcurr'] = 0
I.header['zmu2'] = GB
# set normal and skew quads
I.ele['CQ01']['b1_gradient'] = 0.00714 # T/m
I.ele['SQ01']['b1_gradient'] = 0
Single particle tracking¶
In [5]:
Copied!
# Track
I2 = I.copy()
I2.configure()
# Track
I2 = I.copy()
I2.configure()
In [6]:
Copied!
ele = I2.ele['CQ01']
ele
ele = I2.ele['CQ01']
ele
Out[6]:
{'description': 'name:CQ01', 'original': '0.36 0 0 1 0.01601 0.00714 0.210 0.0254 0.0 0.0 0.0 0.0 0 /!name:CQ01', 'L': 0.36, 'type': 'quadrupole', 'zedge': 0.01601, 'b1_gradient': 0.00714, 'L_effective': 0.21, 'radius': 0.0254, 'x_offset': 0.0, 'y_offset': 0.0, 'x_rotation': 0.0, 'y_rotation': 0.0, 'z_rotation': 0.0, 's': 0.37601, 'name': 'CQ01'}
In [7]:
Copied!
# Estimate for angle change for a 6 MeV/c momentum particle, offset by 1 mm.
ele['b1_gradient']*ele['L_effective']*299792458 / 6e6 * .001
# Estimate for angle change for a 6 MeV/c momentum particle, offset by 1 mm.
ele['b1_gradient']*ele['L_effective']*299792458 / 6e6 * .001
Out[7]:
7.491813525419999e-05
In [8]:
Copied!
P2 = I2.track1(s=2.2, z0 = 0, x0=0.001, pz0=6e6)
P2.xp
P2 = I2.track1(s=2.2, z0 = 0, x0=0.001, pz0=6e6)
P2.xp
Out[8]:
array([7.51244699e-05])
In [9]:
Copied!
I2.plot('mean_x')
I2.plot('mean_x')
Track beam¶
In [10]:
Copied!
# Regular and Skew quads
I.run()
# Regular and Skew quads
I.run()
In [11]:
Copied!
I.output['stats'].keys()
I.output['stats'].keys()
Out[11]:
dict_keys(['t', 'mean_z', 'mean_gamma', 'mean_beta', 'max_r', 'sigma_gamma', 'mean_x', 'sigma_x', 'norm_emit_x', 'mean_y', 'sigma_y', 'norm_emit_y', 'sigma_z', 'norm_emit_z', 'max_amplitude_x', 'max_amplitude_y', 'max_amplitude_z', 'loadbalance_min_n_particle', 'loadbalance_max_n_particle', 'n_particle', 'moment3_x', 'moment3_y', 'moment3_z', 'moment4_x', 'moment4_y', 'moment4_z', 'cov_x__x', 'cov_x__y', 'cov_x__z', 'cov_y__y', 'cov_y__z', 'cov_z__z', 'mean_kinetic_energy', 'cov_x__px', 'cov_y__py', 'cov_z__pz', 'cov_x__py', 'cov_x__pz', 'cov_px__px', 'cov_y__px', 'cov_px__py', 'cov_z__px', 'cov_px__pz', 'cov_y__pz', 'cov_py__py', 'cov_z__py', 'cov_py__pz', 'cov_pz__pz'])
In [12]:
Copied!
PI = I.particles['initial_particles']
PF = I.particles['final_particles']
PI['sigma_y']
PI = I.particles['initial_particles']
PF = I.particles['final_particles']
PI['sigma_y']
Out[12]:
0.001000777236945671
In [13]:
Copied!
# Compare these.
key1 = 'mean_z'
key2 = 'sigma_x'
units1 = str(I.units(key1))
units2 = str(I.units(key2))
plt.xlabel(key1+f' ({units1})')
plt.ylabel(key2+f' ({units2})')
plt.plot(I.stat(key1), I.stat(key2))
plt.scatter(
[I.particles[name][key1] for name in I.particles],
[I.particles[name][key2] for name in I.particles], color='red')
# Compare these.
key1 = 'mean_z'
key2 = 'sigma_x'
units1 = str(I.units(key1))
units2 = str(I.units(key2))
plt.xlabel(key1+f' ({units1})')
plt.ylabel(key2+f' ({units2})')
plt.plot(I.stat(key1), I.stat(key2))
plt.scatter(
[I.particles[name][key1] for name in I.particles],
[I.particles[name][key2] for name in I.particles], color='red')
Out[13]:
<matplotlib.collections.PathCollection at 0x12e300e30>
In [14]:
Copied!
# Compare these.
key1 = 'mean_z'
key2 = 'sigma_x'
units1 = str(I.units(key1))
units2 = str(I.units(key2))
plt.xlabel(key1+f' ({units1})')
plt.ylabel(key2+f' ({units2})')
plt.plot(I.stat(key1), I.stat(key2))
plt.scatter(
[I.particles[name][key1] for name in I.particles],
[I.particles[name][key2] for name in I.particles], color='red')
key2 = 'sigma_y'
plt.plot(I.stat(key1), I.stat(key2))
plt.scatter(
[I.particles[name][key1] for name in I.particles],
[I.particles[name][key2] for name in I.particles], color='green')
# Compare these.
key1 = 'mean_z'
key2 = 'sigma_x'
units1 = str(I.units(key1))
units2 = str(I.units(key2))
plt.xlabel(key1+f' ({units1})')
plt.ylabel(key2+f' ({units2})')
plt.plot(I.stat(key1), I.stat(key2))
plt.scatter(
[I.particles[name][key1] for name in I.particles],
[I.particles[name][key2] for name in I.particles], color='red')
key2 = 'sigma_y'
plt.plot(I.stat(key1), I.stat(key2))
plt.scatter(
[I.particles[name][key1] for name in I.particles],
[I.particles[name][key2] for name in I.particles], color='green')
Out[14]:
<matplotlib.collections.PathCollection at 0x12e4b9070>
In [15]:
Copied!
PF.plot('x', 'y')
PF.plot('delta_z', 'delta_pz')
PF.plot('x', 'y')
PF.plot('delta_z', 'delta_pz')