Autophase and Autophase and Scale examples¶

This includes general (slow) autophasing and scaling, as well as fast autophasing.

In [ ]:

%load_ext autoreload
%autoreload 2

%load_ext autoreload %autoreload 2

In [ ]:

from impact import Impact

import numpy as np

# Nicer plotting
import matplotlib.pyplot as plt
import matplotlib

%config InlineBackend.figure_format = 'retina'
matplotlib.rcParams["figure.figsize"] = (8, 4)

NUMPROCS = 8

from impact import Impact import numpy as np # Nicer plotting import matplotlib.pyplot as plt import matplotlib %config InlineBackend.figure_format = 'retina' matplotlib.rcParams["figure.figsize"] = (8, 4) NUMPROCS = 8

Make Impact object from the LCLS injector model:

In [ ]:

ifile = "templates/lcls_injector/ImpactT.in"

ifile = "templates/lcls_injector/ImpactT.in"

In [ ]:

I = Impact(ifile, verbose=True)
I.numprocs = 1

I = Impact(ifile, verbose=True) I.numprocs = 1

Phase and Scale the LCLS gun¶

In [ ]:

from impact.autophase import autophase_and_scale

from pmd_beamphysics import single_particle

P0 = single_particle(pz=1e-15, z=1e-15)

from impact.autophase import autophase_and_scale from pmd_beamphysics import single_particle P0 = single_particle(pz=1e-15, z=1e-15)

In [ ]:

autophase_and_scale(
    I,
    phase_ele_name="GUN",
    target=6e6,
    phase_range=(270, 360),
    scale_range=(10e6, 100e6),
    initial_particles=P0,
    verbose=True,
)

autophase_and_scale( I, phase_ele_name="GUN", target=6e6, phase_range=(270, 360), scale_range=(10e6, 100e6), initial_particles=P0, verbose=True, )

Check the energy:

In [ ]:

I.verbose = False
PF = I.track(P0, s=0.15)
PF["mean_energy"]

I.verbose = False PF = I.track(P0, s=0.15) PF["mean_energy"]

Examine this process using the debug flag. This will return the function used for phasing and scaling.

In [ ]:

ps_f, Itest = autophase_and_scale(
    I, phase_ele_name="GUN", target=6e6, initial_particles=P0, verbose=False, debug=True
)

ps_f, Itest = autophase_and_scale( I, phase_ele_name="GUN", target=6e6, initial_particles=P0, verbose=False, debug=True )

Plot various phases and scales:

In [ ]:

ptry = np.linspace(-100, 50, 30)
for sc in np.linspace(10e6, 100e6, 5):
    res = np.array([ps_f(p, sc) / 1e6 for p in ptry])
    plt.plot(ptry, res, label=f"{sc/1e6:0.2f} MV")
plt.title("Final energy for various phases and scales")
plt.ylabel("Final energy (MeV)")
plt.xlabel("phase (deg)")
plt.legend()

ptry = np.linspace(-100, 50, 30) for sc in np.linspace(10e6, 100e6, 5): res = np.array([ps_f(p, sc) / 1e6 for p in ptry]) plt.plot(ptry, res, label=f"{sc/1e6:0.2f} MV") plt.title("Final energy for various phases and scales") plt.ylabel("Final energy (MeV)") plt.xlabel("phase (deg)") plt.legend()

Make a 3D data and plot the surface

In [ ]:

X = np.linspace(-100, 50, 10)
Y = np.linspace(10e6, 100e6, 10)
X, Y = np.meshgrid(X, Y)


@np.vectorize
def f(phase, scale):
    return ps_f(phase, scale)


Z = f(X, Y)

X = np.linspace(-100, 50, 10) Y = np.linspace(10e6, 100e6, 10) X, Y = np.meshgrid(X, Y) @np.vectorize def f(phase, scale): return ps_f(phase, scale) Z = f(X, Y)

In [ ]:

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")

surf = ax.plot_surface(
    X, Y / 1e6, Z / 1e6, cmap=matplotlib.cm.coolwarm, linewidth=0, antialiased=True
)

# Add a color bar which maps values to colors.
# fig.colorbar(surf, shrink=0.5, aspect=5)

ax.set_xlabel("phase (deg)")
ax.set_ylabel("scale (MV)")
ax.set_zlabel("Final energy (MeV)")
plt.show()

fig = plt.figure() ax = fig.add_subplot(111, projection="3d") surf = ax.plot_surface( X, Y / 1e6, Z / 1e6, cmap=matplotlib.cm.coolwarm, linewidth=0, antialiased=True ) # Add a color bar which maps values to colors. # fig.colorbar(surf, shrink=0.5, aspect=5) ax.set_xlabel("phase (deg)") ax.set_ylabel("scale (MV)") ax.set_zlabel("Final energy (MeV)") plt.show()

Phase and scale LCLS linac sections¶

Linacs L0A and L0B are special, because they require 4 fieldmaps each to model the travelling wave structure. To tune these together, we need to add control groups.

These will control overall phases:

In [ ]:

I.add_group(
    "L0A",
    ele_names=["L0A_entrance", "L0A_body_1", "L0A_body_2", "L0A_exit"],
    var_name="theta0_deg",
    attributes="theta0_deg",
)
I.add_group(
    "L0B",
    ele_names=["L0B_entrance", "L0B_body_1", "L0B_body_2", "L0B_exit"],
    var_name="theta0_deg",
    attributes="theta0_deg",
)

I.add_group( "L0A", ele_names=["L0A_entrance", "L0A_body_1", "L0A_body_2", "L0A_exit"], var_name="theta0_deg", attributes="theta0_deg", ) I.add_group( "L0B", ele_names=["L0B_entrance", "L0B_body_1", "L0B_body_2", "L0B_exit"], var_name="theta0_deg", attributes="theta0_deg", )

These will control overall scaling, respecting the special factors:

In [ ]:

I.add_group(
    "L0A_scale",
    ele_names=["L0A_entrance", "L0A_body_1", "L0A_body_2", "L0A_exit"],
    var_name="rf_field_scale",
    factors=[0.86571945106805, 1, 1, 0.86571945106805],  # sin(k*d) with d = 3.5e-2 m
    absolute=True,
)

I.add_group(
    "L0B_scale",
    ele_names=["L0B_entrance", "L0B_body_1", "L0B_body_2", "L0B_exit"],
    var_name="rf_field_scale",
    factors=[0.86571945106805, 1, 1, 0.86571945106805],  # sin(k*d) with d = 3.5e-2 m
    absolute=True,
)

I["L0A_scale"]["rf_field_scale"] = 30e6
# I['L0A_scale'].__dict__

I.add_group( "L0A_scale", ele_names=["L0A_entrance", "L0A_body_1", "L0A_body_2", "L0A_exit"], var_name="rf_field_scale", factors=[0.86571945106805, 1, 1, 0.86571945106805], # sin(k*d) with d = 3.5e-2 m absolute=True, ) I.add_group( "L0B_scale", ele_names=["L0B_entrance", "L0B_body_1", "L0B_body_2", "L0B_exit"], var_name="rf_field_scale", factors=[0.86571945106805, 1, 1, 0.86571945106805], # sin(k*d) with d = 3.5e-2 m absolute=True, ) I["L0A_scale"]["rf_field_scale"] = 30e6 # I['L0A_scale'].__dict__

Now phase and scale L0A to 64 MeV:

In [ ]:

res_L0A = autophase_and_scale(
    I,
    phase_ele_name="L0A",
    scale_ele_name="L0A_scale",
    target=64e6,
    scale_range=(10e6, 100e6),
    initial_particles=P0,
    verbose=True,
)

res_L0A = autophase_and_scale( I, phase_ele_name="L0A", scale_ele_name="L0A_scale", target=64e6, scale_range=(10e6, 100e6), initial_particles=P0, verbose=True, )

Do the same for L0B:

In [ ]:

autophase_and_scale(
    I,
    phase_ele_name="L0B",
    scale_ele_name="L0B_scale",
    target=135e6,
    scale_range=(10e6, 100e6),
    initial_particles=P0,
    verbose=True,
)

autophase_and_scale( I, phase_ele_name="L0B", scale_ele_name="L0B_scale", target=135e6, scale_range=(10e6, 100e6), initial_particles=P0, verbose=True, )

Check the final energy and plot:

In [ ]:

I.track(P0, s=8.371612)["mean_energy"]

I.track(P0, s=8.371612)["mean_energy"]

In [ ]:

plt.plot(I.stat("mean_z"), I.stat("mean_kinetic_energy") / 1e6 + 0.511)

plt.plot(I.stat("mean_z"), I.stat("mean_kinetic_energy") / 1e6 + 0.511)

Fast autophase¶

This is a faster method that can find and set all relative phases by tracking the fields externally.

In [ ]:

%%time
I.autophase()

%%time I.autophase()

Sending in a dict will set these phases as it goes:

In [ ]:

I.verbose = True

I.autophase({"GUN": 1, "L0A": 2, "L0B": 3})

I.verbose = True I.autophase({"GUN": 1, "L0A": 2, "L0B": 3})

Autophase without scaling¶

Just phasing is simpler.

In [ ]:

from impact.autophase import autophase

ifile2 = "templates/apex_gun/ImpactT.in"

I2 = Impact(ifile2, verbose=False)

from impact.autophase import autophase ifile2 = "templates/apex_gun/ImpactT.in" I2 = Impact(ifile2, verbose=False)

In [ ]:

autophase(
    I2,
    ele_name="APEX_GUN",
    initial_particles=P0,
    metric="mean_kinetic_energy",
    verbose=True,
)

autophase( I2, ele_name="APEX_GUN", initial_particles=P0, metric="mean_kinetic_energy", verbose=True, )

In [ ]:

phase_f, Itest = autophase(
    I2,
    ele_name="APEX_GUN",
    metric="mean_kinetic_energy",
    initial_particles=P0,
    debug=True,
)

phase_f, Itest = autophase( I2, ele_name="APEX_GUN", metric="mean_kinetic_energy", initial_particles=P0, debug=True, )

In [ ]:

# Phases to try
ptry = np.linspace(0, 360, 60)

energies = np.array([phase_f(p) / 1e3 for p in ptry])

plt.plot(ptry, energies)
plt.ylim(0, 800)
plt.title("Final energy for various phases in the APEX gun")
plt.ylabel("Final kinetic energy (keV)")
plt.xlabel("phase (deg)")

# Phases to try ptry = np.linspace(0, 360, 60) energies = np.array([phase_f(p) / 1e3 for p in ptry]) plt.plot(ptry, energies) plt.ylim(0, 800) plt.title("Final energy for various phases in the APEX gun") plt.ylabel("Final kinetic energy (keV)") plt.xlabel("phase (deg)")

Autophase with alternative metric, and bunch tracking with space charge.¶

The above uses mean_energy as the metric to maximize. Alternatively, one might want to minimize energy spread. This is accomplished by passing maximize=False and metric='sigma_pz' or similar.

In [ ]:

from distgen import Generator


ifile = "templates/lcls_injector/ImpactT.in"
gfile = "templates/lcls_injector/distgen.yaml"

G = Generator(gfile)
G["n_particle"] = 2000
G.run()
P0 = G.particles

from distgen import Generator ifile = "templates/lcls_injector/ImpactT.in" gfile = "templates/lcls_injector/distgen.yaml" G = Generator(gfile) G["n_particle"] = 2000 G.run() P0 = G.particles

In [ ]:

%%time
I = Impact(ifile, initial_particles=P0, verbose=False)
I.stop = 0.16
I.numprocs = NUMPROCS
I.run()

%%time I = Impact(ifile, initial_particles=P0, verbose=False) I.stop = 0.16 I.numprocs = NUMPROCS I.run()

In [ ]:

phase_f, Itest = autophase(
    I,
    ele_name="GUN",
    metric="sigma_pz",
    maximize=False,
    initial_particles=P0,
    debug=True,
    verbose=True,
)

phase_f, Itest = autophase( I, ele_name="GUN", metric="sigma_pz", maximize=False, initial_particles=P0, debug=True, verbose=True, )

In [ ]:

I.particles["final_particles"].plot("z", "pz")

I.particles["final_particles"].plot("z", "pz")

In [ ]:

# Phases to try
ptry = np.linspace(290, 310, 20)

sigma_pzs = np.array([phase_f(p) for p in ptry])

plt.plot(ptry, sigma_pzs)
# plt.ylim(0, 800)
# plt.title('Final energy for various phases in the APEX gun')
# plt.ylabel('Final kinetic energy (keV)')
plt.xlabel("phase (deg)")

# Phases to try ptry = np.linspace(290, 310, 20) sigma_pzs = np.array([phase_f(p) for p in ptry]) plt.plot(ptry, sigma_pzs) # plt.ylim(0, 800) # plt.title('Final energy for various phases in the APEX gun') # plt.ylabel('Final kinetic energy (keV)') plt.xlabel("phase (deg)")

In [ ]:

phase_f(293.5)

phase_f(293.5)

In [ ]:

Itest.particles["final_particles"].plot("z", "pz")

Itest.particles["final_particles"].plot("z", "pz")

In [ ]:

phase_f, Itest = autophase(
    I,
    ele_name="GUN",
    metric="sigma_pz",
    maximize=False,
    initial_particles=P0,
    debug=True,
    s_stop=1.45,
    verbose=True,
)
# Phases to try
ptry = np.linspace(270, 290, 30)

sigma_pzs = np.array([phase_f(p) for p in ptry])

plt.plot(ptry, sigma_pzs)
# plt.ylim(0, 800)
# plt.title('Final energy for various phases in the APEX gun')
# plt.ylabel('Final kinetic energy (keV)')
plt.xlabel("phase (deg)")

phase_f, Itest = autophase( I, ele_name="GUN", metric="sigma_pz", maximize=False, initial_particles=P0, debug=True, s_stop=1.45, verbose=True, ) # Phases to try ptry = np.linspace(270, 290, 30) sigma_pzs = np.array([phase_f(p) for p in ptry]) plt.plot(ptry, sigma_pzs) # plt.ylim(0, 800) # plt.title('Final energy for various phases in the APEX gun') # plt.ylabel('Final kinetic energy (keV)') plt.xlabel("phase (deg)")

In [ ]:

phase_f(280.0)

phase_f(280.0)

In [ ]:

Itest.particles["final_particles"].plot("z", "pz")

Itest.particles["final_particles"].plot("z", "pz")