Quadrupole Example¶

Simple quadrupole example

In [1]:

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]:

# 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]:

# 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]:

np.float64(11.699046356605859)

Use Impact's built-in Gaussian particle generator¶

In [4]:

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]:

# Track
I2 = I.copy()
I2.configure()

# Track I2 = I.copy() I2.configure()

In [6]:

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]:

# Estimate for angle change for a 6 MeV/c momentum particle, offset by 1 mm.
ele["b1_gradient"] * ele["L_effective"] * 299792458 / 6e6 * 0.001

# Estimate for angle change for a 6 MeV/c momentum particle, offset by 1 mm. ele["b1_gradient"] * ele["L_effective"] * 299792458 / 6e6 * 0.001

Out[7]:

7.491813525419999e-05

In [8]:

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]:

I2.plot("mean_x")

I2.plot("mean_x")

Track beam¶

In [10]:

# Regular and Skew quads
I.run()

# Regular and Skew quads I.run()

In [11]:

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]:

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]:

np.float64(0.001000777236945671)

In [13]:

# 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 0x7fac58e2ffe0>

In [14]:

# 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 0x7fac584fe7b0>

In [15]:

PF.plot("x", "y")
PF.plot("delta_z", "delta_pz")

PF.plot("x", "y") PF.plot("delta_z", "delta_pz")