Solenoid Example¶

Simple solenoid example

In [1]:

from impact import Impact

import numpy as np
import os

import matplotlib.pyplot as plt

%config InlineBackend.figure_format='retina'

from impact import Impact import numpy as np import os import matplotlib.pyplot as plt %config InlineBackend.figure_format='retina'

In [2]:

ifile1d = "../templates/solenoid/ImpactT_solenoid_1d.in"
ifile2d = "../templates/solenoid/ImpactT_solenoid_2d.in"

os.path.exists(ifile1d), os.path.exists(ifile2d)

ifile1d = "../templates/solenoid/ImpactT_solenoid_1d.in" ifile2d = "../templates/solenoid/ImpactT_solenoid_2d.in" os.path.exists(ifile1d), os.path.exists(ifile2d)

Out[2]:

(True, True)

Use Impact's built-in Gaussian particle generator¶

In [3]:

I1 = Impact(ifile1d)
I2 = Impact(ifile2d)

# Turn off SC
I1["total_charge"] = 0
I2["total_charge"] = 0
print(I1)

I1 = Impact(ifile1d) I2 = Impact(ifile2d) # Turn off SC I1["total_charge"] = 0 I2["total_charge"] = 0 print(I1)

================ Impact-T Summary ================
10000 particles
1 bunch of electrons
total charge: 0.0 pC
Distribution type: gauss3
Free space start
Processor domain: 1 x 1 = 1 CPUs
Space charge grid: 32 x 32 x 32
Maximum time steps: 1000000
Reference Frequency: 1000000000.0 Hz
Initial reference time: 0.0 s
Simulation starting from the beginning
=================================================
Impact-T configured in /tmp/tmpocpimvqp

In [4]:

%%time
I1["total_charge"] = 0
I1.run()
I1.plot()

%%time I1["total_charge"] = 0 I1.run() I1.plot()

CPU times: user 107 ms, sys: 19.3 ms, total: 127 ms
Wall time: 10.5 s

In [5]:

%%time
I2["total_charge"] = 0
I2.run()
I2.plot()

%%time I2["total_charge"] = 0 I2.run() I2.plot()

CPU times: user 324 ms, sys: 21.3 ms, total: 346 ms
Wall time: 1.69 s

In [6]:

# The 2D version keeps the field internally as a FieldMesh
I2.fieldmaps["1T912.T7"]["field"].plot_onaxis()

# The 2D version keeps the field internally as a FieldMesh I2.fieldmaps["1T912.T7"]["field"].plot_onaxis()

Single particle tracking¶

In [7]:

%%time
P1 = I1.track1(s=0.5, z0=0, x0=0.019, pz0=3e6)

%%time P1 = I1.track1(s=0.5, z0=0, x0=0.019, pz0=3e6)

CPU times: user 34.6 ms, sys: 8.45 ms, total: 43 ms
Wall time: 75.7 ms

Compare 1D and 2D maps¶

In [8]:

X0 = 0.003

I1.track1(s=0.4, x0=X0, pz0=3e6)
I2.track1(s=0.4, x0=X0, pz0=3e6)

k1 = "mean_z"
k2 = "mean_x"

x1 = I1.stat(k1)
y1 = I1.stat(k2)

x2 = I2.stat(k1)
y2 = I2.stat(k2)

fig, ax = plt.subplots(figsize=(16, 9))

ax.plot(x1, y1, color="black", label="1D fieldmap")
ax.plot(x2, y2, color="red", linestyle="--", label="2D fieldmap")
ax.legend()

X0 = 0.003 I1.track1(s=0.4, x0=X0, pz0=3e6) I2.track1(s=0.4, x0=X0, pz0=3e6) k1 = "mean_z" k2 = "mean_x" x1 = I1.stat(k1) y1 = I1.stat(k2) x2 = I2.stat(k1) y2 = I2.stat(k2) fig, ax = plt.subplots(figsize=(16, 9)) ax.plot(x1, y1, color="black", label="1D fieldmap") ax.plot(x2, y2, color="red", linestyle="--", label="2D fieldmap") ax.legend()

Out[8]:

<matplotlib.legend.Legend at 0x7febd2655400>

In [9]:

I2.ele["SOL1"]

I2.ele["SOL1"]

Out[9]:

{'description': 'name:SOL1',
 'original': ' 0.4 0 0 3 -0.02678 .2 912 0 0 0 0 0 0 /!name:SOL1',
 'L': 0.4,
 'type': 'solenoid',
 'zedge': -0.02678,
 'b_field': 0.2,
 'filename': '1T912.T7',
 'radius': 0.0,
 's': 0.37322,
 'name': 'SOL1'}

In [10]:

fig, ax = plt.subplots(figsize=(12, 8))

k1 = "mean_z"
k2 = "mean_x"

f1 = 1e3
f2 = 1e3
u1 = "mm"
u2 = "mm"

for X0 in np.linspace(0, 0.018, 10):
    I1.track1(s=0.4, x0=X0, pz0=3e6)
    I2.track1(s=0.4, x0=X0, pz0=3e6)

    x1 = I1.stat(k1)
    y1 = I1.stat(k2)

    x2 = I2.stat(k1)
    y2 = I2.stat(k2)

    if X0 == 0:
        label1 = "1D fieldmap"
        label2 = "2D fieldmap"
    else:
        label1 = None
        label2 = None

    ax.plot(x1 * f1, y1 * f2, color="black", label=label1)
    ax.plot(x2 * f1, y2 * f2, color="red", linestyle="--", label=label2)

ax.set_ylim(0, 18)
ax.set_xlabel(f"{k1} ({u1})")
ax.set_ylabel(f"{k2} ({u2})")
ax.legend()

fig, ax = plt.subplots(figsize=(12, 8)) k1 = "mean_z" k2 = "mean_x" f1 = 1e3 f2 = 1e3 u1 = "mm" u2 = "mm" for X0 in np.linspace(0, 0.018, 10): I1.track1(s=0.4, x0=X0, pz0=3e6) I2.track1(s=0.4, x0=X0, pz0=3e6) x1 = I1.stat(k1) y1 = I1.stat(k2) x2 = I2.stat(k1) y2 = I2.stat(k2) if X0 == 0: label1 = "1D fieldmap" label2 = "2D fieldmap" else: label1 = None label2 = None ax.plot(x1 * f1, y1 * f2, color="black", label=label1) ax.plot(x2 * f1, y2 * f2, color="red", linestyle="--", label=label2) ax.set_ylim(0, 18) ax.set_xlabel(f"{k1} ({u1})") ax.set_ylabel(f"{k2} ({u2})") ax.legend()

Out[10]:

<matplotlib.legend.Legend at 0x7febd2614920>

Track beam¶

In [11]:

I1 = Impact(ifile1d)
I2 = Impact(ifile2d)

# Turn off SC
I1["total_charge"] = 0
I2["total_charge"] = 0
I1.run()
I2.run()

I1 = Impact(ifile1d) I2 = Impact(ifile2d) # Turn off SC I1["total_charge"] = 0 I2["total_charge"] = 0 I1.run() I2.run()

In [12]:

I1.output["stats"].keys()

I1.output["stats"].keys()

Out[12]:

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

PI = I1.particles["initial_particles"]
PF = I1.particles["final_particles"]
PI["sigma_y"]

PI = I1.particles["initial_particles"] PF = I1.particles["final_particles"] PI["sigma_y"]

Out[13]:

np.float64(0.0010106410879665382)

In [14]:

# Compare these.
key1 = "mean_z"
key2 = "sigma_x"
units1 = str(I1.units(key1))
units2 = str(I1.units(key2))
plt.xlabel(key1 + f" ({units1})")
plt.ylabel(key2 + f" ({units2})")
plt.plot(I1.stat(key1), I1.stat(key2))
plt.scatter(
    [I1.particles[name][key1] for name in I1.particles],
    [I2.particles[name][key2] for name in I2.particles],
    color="red",
)

# Compare these. key1 = "mean_z" key2 = "sigma_x" units1 = str(I1.units(key1)) units2 = str(I1.units(key2)) plt.xlabel(key1 + f" ({units1})") plt.ylabel(key2 + f" ({units2})") plt.plot(I1.stat(key1), I1.stat(key2)) plt.scatter( [I1.particles[name][key1] for name in I1.particles], [I2.particles[name][key2] for name in I2.particles], color="red", )

Out[14]:

<matplotlib.collections.PathCollection at 0x7febd33bfe60>

In [15]:

# Compare these.
key1 = "mean_z"
key2 = "sigma_x"
units1 = str(I1.units(key1))
units2 = str(I1.units(key2))
plt.xlabel(key1 + f" ({units1})")
plt.ylabel(key2 + f" ({units2})")
plt.plot(I1.stat(key1), I1.stat(key2), label="1D solenoid")
plt.scatter(
    [I1.particles[name][key1] for name in I1.particles],
    [I1.particles[name][key2] for name in I1.particles],
    color="red",
)
key2 = "sigma_y"
plt.plot(I2.stat(key1), I2.stat(key2), label="2D solenoid")
plt.scatter(
    [I2.particles[name][key1] for name in I2.particles],
    [I2.particles[name][key2] for name in I2.particles],
    color="green",
)
plt.legend()

# Compare these. key1 = "mean_z" key2 = "sigma_x" units1 = str(I1.units(key1)) units2 = str(I1.units(key2)) plt.xlabel(key1 + f" ({units1})") plt.ylabel(key2 + f" ({units2})") plt.plot(I1.stat(key1), I1.stat(key2), label="1D solenoid") plt.scatter( [I1.particles[name][key1] for name in I1.particles], [I1.particles[name][key2] for name in I1.particles], color="red", ) key2 = "sigma_y" plt.plot(I2.stat(key1), I2.stat(key2), label="2D solenoid") plt.scatter( [I2.particles[name][key1] for name in I2.particles], [I2.particles[name][key2] for name in I2.particles], color="green", ) plt.legend()

Out[15]:

<matplotlib.legend.Legend at 0x7febd1e9d400>

In [16]:

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

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