Wakefield example¶

Simple 1 m drift with a wakefield.

This verifies that the analytic formula uses is SLAC-PUB-9663 Eq. 8

In [1]:

import numpy as np
import os

import matplotlib.pyplot as plt
import matplotlib

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

import numpy as np import os import matplotlib.pyplot as plt import matplotlib matplotlib.rcParams["figure.figsize"] = (8, 4) %config InlineBackend.figure_format='retina'

In [2]:

# locate the drift template
from impact import Impact

ifile = "../templates/wakefield/ImpactT.in"
os.path.exists(ifile)

# locate the drift template from impact import Impact ifile = "../templates/wakefield/ImpactT.in" os.path.exists(ifile)

Out[2]:

True

In [3]:

# gamma for 1 GeV
1e9 / 0.511e6

# gamma for 1 GeV 1e9 / 0.511e6

Out[3]:

1956.9471624266146

Use Impact's built-in Gaussian particle generator¶

In [4]:

I = Impact(ifile)
I.header["Np"] = 10000
I.header["Nx"] = 32
I.header["Ny"] = 32
I.header["Nz"] = 32
I.header["Dt"] = 10e-12

I = Impact(ifile) I.header["Np"] = 10000 I.header["Nx"] = 32 I.header["Ny"] = 32 I.header["Nz"] = 32 I.header["Dt"] = 10e-12

In [5]:

I.lattice

I.lattice

Out[5]:

[{'description': 'name:sc_off',
  'original': '0 0 0 -8 0 -1 0 /!name:sc_off',
  'type': 'spacecharge',
  's': 0.0,
  'is_on': False,
  'name': 'sc_off'},
 {'type': 'comment',
  'description': '!0 0 0 -5 0 0 -1000.0 /!name:2d_to_3d_spacecharge',
  'name': '2d_to_3d_spacecharge'},
 {'description': 'name:wakefield_1',
  'original': '0 -1 0 -6 1 1 0.0 1.0 0.0116 0.0292 0.035 /!name:wakefield_1',
  'type': 'wakefield',
  's_begin': 0.0,
  's': 1.0,
  'method': 'analytical',
  'iris_radius': 0.0116,
  'gap': 0.0292,
  'period': 0.035,
  'name': 'wakefield_1'},
 {'description': 'name:drift_1',
  'original': '1.0 0 0 0 1.0 0.15 /!name:drift_1',
  'L': 1.0,
  'type': 'drift',
  'zedge': 1.0,
  'radius': 0.15,
  's': 2.0,
  'name': 'drift_1'},
 {'description': 'name:stop_1',
  'original': '0 0 0 -99 0 0.0 1 /!name:stop_1',
  'type': 'stop',
  's': 1.0,
  'name': 'stop_1'},
 {'type': 'comment', 'description': '', 'name': 'comment_2'}]

In [6]:

I.run()

I.run()

In [7]:

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

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

Out[7]:

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

PI = I.particles["initial_particles"]
PF = I.particles["final_particles"]
PI, PF

PI = I.particles["initial_particles"] PF = I.particles["final_particles"] PI, PF

Out[8]:

(<ParticleGroup with 10000 particles at 0x7f511e68a9c0>,
 <ParticleGroup with 10000 particles at 0x7f511e6e3b30>)

In [9]:

PI.plot("delta_z", "delta_pz")
PF.plot("delta_z", "delta_pz")

PI.plot("delta_z", "delta_pz") PF.plot("delta_z", "delta_pz")

In [10]:

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

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

In [11]:

# np.savetxt('/Users/chrisonian/Scratch/impactwake.dat', np.array([PF['z'], PF['pz']]).T)

# np.savetxt('/Users/chrisonian/Scratch/impactwake.dat', np.array([PF['z'], PF['pz']]).T)

Make particles in distgen¶

In [12]:

from distgen import Generator

YAML = """
n_particle: 20000
random_type: hammersley
species: electron
start:
  tstart:
    units: sec
    value: 0
  type: time
total_charge:
  units: nC
  value: 1
r_dist:
  sigma_xy:
    units: mm
    value: .001
  type: radial_gaussian
z_dist:
  avg_z:
    units: mm
    value: 0
  sigma_z:
    units: mm
    value: 0.1
  type: gaussian
  

transforms:
  setPz:
    type: set_avg pz
    avg_pz: 
      value: 1
      units: GeV/c
  
"""
G = Generator(YAML)
G.run()
P = G.particles

from distgen import Generator YAML = """ n_particle: 20000 random_type: hammersley species: electron start: tstart: units: sec value: 0 type: time total_charge: units: nC value: 1 r_dist: sigma_xy: units: mm value: .001 type: radial_gaussian z_dist: avg_z: units: mm value: 0 sigma_z: units: mm value: 0.1 type: gaussian transforms: setPz: type: set_avg pz avg_pz: value: 1 units: GeV/c """ G = Generator(YAML) G.run() P = G.particles

In [13]:

I = Impact(ifile, initial_particles=P, verbose=False)
I.run()
PF2 = I.particles["final_particles"]

I = Impact(ifile, initial_particles=P, verbose=False) I.run() PF2 = I.particles["final_particles"]

In [14]:

PF2.plot("x", "px")
PF2.plot("delta_z", "delta_pz")

PF2.plot("x", "px") PF2.plot("delta_z", "delta_pz")

Compare¶

In [15]:

for k in ["x", "z", "pz"]:
    plt.hist(PF[k], density=True, bins=100, label="Impact-T generator", alpha=0.5)
    plt.hist(PF2[k], density=True, bins=100, label="Distgen generator", alpha=0.5)
    plt.xlabel(k)
    plt.legend()
    plt.show()

for k in ["x", "z", "pz"]: plt.hist(PF[k], density=True, bins=100, label="Impact-T generator", alpha=0.5) plt.hist(PF2[k], density=True, bins=100, label="Distgen generator", alpha=0.5) plt.xlabel(k) plt.legend() plt.show()

Checking the wakefield with SLAC-PUB-9663¶

Impact-T seems to use Eq. * from SLAC-PUB-9663, Karl Bane (2003).

https://www.slac.stanford.edu/pubs/slacpubs/9500/slac-pub-9663.pdf

In [16]:

# Define alpha function for the s00 calc.


def alpha(g):
    """
    SLAC-PUB-9663 equation (5)

    """
    a1 = 0.4648
    return 1 - a1 * np.sqrt(g) - (1 - 2 * a1) * g


def bane_wake(z, a=0.0116, g=0.0292, L=0.035):
    s00 = g / 8 * (a / (alpha(g / L) * L)) ** 2

    #   'iris_radius': 0.0116,
    #   'gap': 0.0292,
    #   'period': 0.035,

    Z0c_over_pi = 120 * 299792458.0  # Ohm m/s

    WL = Z0c_over_pi / a**2 * np.exp(-np.sqrt(z / s00))

    return WL


def bane_wake2(z, a=0.0116, g=0.0292, L=0.035):
    """
    From SLAC-PUB-11829
    """

    s1 = 0.41 * a**1.8 * g**1.6 / L**2.4

    Z0c_over_pi = 120 * 299792458.0  # Ohm m/s

    WL = Z0c_over_pi / a**2 * np.exp(-np.sqrt(z / s1))

    return WL


plt.xlabel("z (m)")
plt.ylabel("Wake (V/C)")
plt.yscale("log")

dzz = 0.00001
zz = np.arange(0, 0.01, dzz)

plt.yscale("log")
plt.plot(zz, bane_wake(zz), label="SLAC-PUB-9663", color="red")
plt.plot(zz, bane_wake2(zz), label="SLAC-PUB-11829", color="green")
plt.legend()

# Define alpha function for the s00 calc. def alpha(g): """ SLAC-PUB-9663 equation (5) """ a1 = 0.4648 return 1 - a1 * np.sqrt(g) - (1 - 2 * a1) * g def bane_wake(z, a=0.0116, g=0.0292, L=0.035): s00 = g / 8 * (a / (alpha(g / L) * L)) ** 2 # 'iris_radius': 0.0116, # 'gap': 0.0292, # 'period': 0.035, Z0c_over_pi = 120 * 299792458.0 # Ohm m/s WL = Z0c_over_pi / a**2 * np.exp(-np.sqrt(z / s00)) return WL def bane_wake2(z, a=0.0116, g=0.0292, L=0.035): """ From SLAC-PUB-11829 """ s1 = 0.41 * a**1.8 * g**1.6 / L**2.4 Z0c_over_pi = 120 * 299792458.0 # Ohm m/s WL = Z0c_over_pi / a**2 * np.exp(-np.sqrt(z / s1)) return WL plt.xlabel("z (m)") plt.ylabel("Wake (V/C)") plt.yscale("log") dzz = 0.00001 zz = np.arange(0, 0.01, dzz) plt.yscale("log") plt.plot(zz, bane_wake(zz), label="SLAC-PUB-9663", color="red") plt.plot(zz, bane_wake2(zz), label="SLAC-PUB-11829", color="green") plt.legend()

Out[16]:

<matplotlib.legend.Legend at 0x7f5111bcb410>

In [17]:

# Compare with particles
sigma = 0.0001
Qtot = -1e-9  # C


def density(z, sigma=0.0001):
    return 1 / (np.sqrt(2 * np.pi) * sigma) * np.exp(-0.5 * (z / sigma) ** 2)


dz = sigma / 10
zlist = np.arange(-6 * sigma, 6 * sigma, dz)

# Check normalization
np.sum(density(zlist)) * dz

# Compare with particles sigma = 0.0001 Qtot = -1e-9 # C def density(z, sigma=0.0001): return 1 / (np.sqrt(2 * np.pi) * sigma) * np.exp(-0.5 * (z / sigma) ** 2) dz = sigma / 10 zlist = np.arange(-6 * sigma, 6 * sigma, dz) # Check normalization np.sum(density(zlist)) * dz

Out[17]:

np.float64(0.9999999979663955)

In [18]:

def total_bane_wake(z):
    W = bane_wake(zz)
    return np.sum(W * density(zz + z) * dzz)

def total_bane_wake(z): W = bane_wake(zz) return np.sum(W * density(zz + z) * dzz)

In [19]:

plt.xlabel("z")
plt.ylabel(r"$\Delta p_z$ (eV/c)")
plt.scatter(
    PF["delta_z"], PF["pz"] - PF["max_pz"], marker="x", label="Impact-T tracking"
)
plt.plot(
    zlist,
    Qtot * np.array([total_bane_wake(z) for z in zlist]),
    color="blue",
    label="SLAC-PUB-9663 equation (8)",
)
plt.title("Integrated total wake comparison")
plt.legend()

plt.xlabel("z") plt.ylabel(r"$\Delta p_z$ (eV/c)") plt.scatter( PF["delta_z"], PF["pz"] - PF["max_pz"], marker="x", label="Impact-T tracking" ) plt.plot( zlist, Qtot * np.array([total_bane_wake(z) for z in zlist]), color="blue", label="SLAC-PUB-9663 equation (8)", ) plt.title("Integrated total wake comparison") plt.legend()

Out[19]:

<matplotlib.legend.Legend at 0x7f5115d25550>

Comparison with Wakefield file¶

Many codes will use a wakefield file, with a list of z and single particle wake in V/C

In [20]:

wfile = "Sz_p5um_10mm_per35mm_cell.sdds"
reffile = os.path.join("../templates/wakefield", wfile)
reffile

wfile = "Sz_p5um_10mm_per35mm_cell.sdds" reffile = os.path.join("../templates/wakefield", wfile) reffile

Out[20]:

'../templates/wakefield/Sz_p5um_10mm_per35mm_cell.sdds'

In [21]:

!head -n 8 ../templates/wakefield/Sz_p5um_10mm_per35mm_cell.sdds

!head -n 8 ../templates/wakefield/Sz_p5um_10mm_per35mm_cell.sdds

SDDS1
&column name=z, units=m, type=double,  &end
&column name=W, units=V/C, type=double,  &end
&column name=t, units=s, type=double,  &end
&data mode=ascii, &end
! page number 1
               20001
0.00000000e+000  9.22430234e+012  0.00000000e+000

In [22]:

# Load the file
edat = np.loadtxt(reffile, skiprows=7).T
zw = edat[0]
dzw = np.mean(np.diff(zw))
W_from_file = edat[1] / 35e-3  # Convert to per m

plt.ylabel("W (V/C)")
plt.xlabel("z (m)")
plt.yscale("log")
plt.plot(zw, W_from_file, label=wfile)
plt.plot(zw, np.array([bane_wake(z) for z in zw]), label="SLAC-PUB-9663 equation (8)")
plt.scatter(
    zw,
    np.array([bane_wake(z) for z in zw]),
    label="SLAC-PUB-11829 equation (12)",
    color="red",
)
plt.legend()

# Load the file edat = np.loadtxt(reffile, skiprows=7).T zw = edat[0] dzw = np.mean(np.diff(zw)) W_from_file = edat[1] / 35e-3 # Convert to per m plt.ylabel("W (V/C)") plt.xlabel("z (m)") plt.yscale("log") plt.plot(zw, W_from_file, label=wfile) plt.plot(zw, np.array([bane_wake(z) for z in zw]), label="SLAC-PUB-9663 equation (8)") plt.scatter( zw, np.array([bane_wake(z) for z in zw]), label="SLAC-PUB-11829 equation (12)", color="red", ) plt.legend()

Out[22]:

<matplotlib.legend.Legend at 0x7f511e4dcbc0>

In [23]:

def total_wake_from_file(z):
    return np.sum(W_from_file * density(zw + z) * dzw)


total_wake_from_file(0) * 1e-9

def total_wake_from_file(z): return np.sum(W_from_file * density(zw + z) * dzw) total_wake_from_file(0) * 1e-9

Out[23]:

np.float64(106087.97587646772)

In [24]:

plt.xlabel("z (m)")
plt.ylabel(r"$\Delta p_z$ (eV/c)")
plt.scatter(
    PF["delta_z"], PF["pz"] - PF["max_pz"], marker="x", label="Impact-T tracking"
)
plt.plot(
    zlist,
    Qtot * np.array([total_wake_from_file(z) for z in zlist]),
    color="red",
    label=wfile,
)
plt.plot(
    zlist,
    Qtot * np.array([total_bane_wake(z) for z in zlist]),
    color="blue",
    label="SLAC-PUB-9663 equation (8)",
)
plt.legend()

plt.xlabel("z (m)") plt.ylabel(r"$\Delta p_z$ (eV/c)") plt.scatter( PF["delta_z"], PF["pz"] - PF["max_pz"], marker="x", label="Impact-T tracking" ) plt.plot( zlist, Qtot * np.array([total_wake_from_file(z) for z in zlist]), color="red", label=wfile, ) plt.plot( zlist, Qtot * np.array([total_bane_wake(z) for z in zlist]), color="blue", label="SLAC-PUB-9663 equation (8)", ) plt.legend()

Out[24]:

<matplotlib.legend.Legend at 0x7f5115d27500>