Algorithm & Features

The basic principles calculation in Python of KENV operation are shown here.

Including libraries

First connect the necessary libraries and configure the graphics parameters.

import numpy as np
import holoviews as hv
import warnings

hv.extension('matplotlib')
%output size=200 backend='matplotlib' fig='svg'
%opts Curve Scatter Area [fontsize={'title':14, 'xlabel':14, 'ylabel':14, 'ticks':14, 'legend': 14}]
%opts Area Curve [aspect=3 show_grid=True]
%opts Area  (linewidth=1 alpha=0.25)
%opts Curve (linewidth=2 alpha=0.5)

warnings.filterwarnings('ignore')

Calculation parameters

Let’s set up the area and the counting step.

dz = 0.005 # m
z_min = 0.7
z_max = 15
z = np.arange(z_min,z_max,dz) # m

Define an accelerator

Define an acceleration element (e.g. solenoid, or acceleration gap). The accelerator will be represented as a sequence of such elements.

class Element:
    def __init__(self, z0, MaxField, filename, name):
        self.z0 = z0
        self.MaxField = MaxField
        self.filename = filename
        self.name = name

Focusing solenoids:

Bz_beamline = {}
for   z0,         B0,        filename,       name in [
    # m           T                   Unique name
    [ 0.95,      0.001,    'Bz.dat', 'Sol. 1' ],
    [ 2.10,      0.030,    'Bz.dat', 'Sol. 2' ],
    [ 2.90,      0.001,    'Bz.dat', 'Sol. 3' ],
    [ 4.00,      0.030,    'Bz.dat', 'Sol. 4' ],
    [ 5.64,      0.030,    'Bz.dat', 'Sol. 5' ],
    [ 6.76,      0.001,    'Bz.dat', 'Sol. 6' ],
]:
    Bz_beamline[name] = Element(z0, B0, filename, name)

Bz_beamline['Sol. 3'].z0

2.9

Accelerating elements:

Ez_beamline = {}
for   z0,          E0,       filename,      name in [
    # m            MV/m                     Unique name
    [ 7.456,       -0.9,     'Ez.dat',  'Cavity 3'],
    [ 8.838,       -0.9,     'Ez.dat',  'Cavity 4'],
    [ 10.220,      -0.9,     'Ez.dat',  'Cavity 5'],
    [ 11.602,      -0.9,     'Ez.dat',  'Cavity 6'],
    [ 12.984,      -0.9,     'Ez.dat',  'Cavity 7'],
    [ 14.366,      -0.9,     'Ez.dat',  'Cavity 8'],
]:
    Ez_beamline[name] = Element(z0, E0, filename, name)

Reading field profiles

First, we just read the fields from the files and merge them into one array:

from scipy import interpolate

FieldFiles = {} # buffer for field files
def read_field(beamline, z):
    global FieldFiles
    F = 0
    for element in beamline.values():
        if not (element.filename in FieldFiles):
            print('Reading ' + element.filename)
            FieldFiles[element.filename] = np.loadtxt(
                element.filename
            )
        M = FieldFiles[element.filename]
        z_data = M[:,0]
        F_data = M[:,1]
        f = interpolate.interp1d(
            element.z0+z_data, element.MaxField*F_data,
            fill_value=(0, 0), bounds_error=False
        )
        F = F + f(z)
    return F

Bz = read_field(Bz_beamline, z)

Reading Bz.dat

To integrate the equation into the envelope, convert the array Bz into the function Bz(z)

Bz = interpolate.interp1d(z, Bz, fill_value=(0, 0),
                          bounds_error=False)

Now Bz is already the function that can quickly deduce the field at any points along the gas pedal.

Bz(2.11111)

array(0.02984424)

Let us construct the field \(Bz(z):\)

dim_z  = hv.Dimension('z',  unit='m', range=(0,z_max))
dim_Bz = hv.Dimension('Bz', unit='T', label='$B_z$')

z_Bz = hv.Area((z,Bz(z)), kdims=dim_z, vdims=dim_Bz)
z_Bz

Similarly for the accelerating electric field:

Ez = 1*read_field(Ez_beamline, z)

Reading Ez.dat

Ez = interpolate.interp1d(z, Ez, fill_value=(0, 0),
                          bounds_error=False)

dim_Ez = hv.Dimension('Ez', unit='MV/m', label=r'$E_z$')

z_Ez = hv.Area((z,Ez(z)), kdims=dim_z, vdims=dim_Ez)

z_Ez

Longitudinal beam dynamics

The equation for the longitudinal dynamics of the beam can be solved independently of the envelope equation, so that in the envelope equation the ready function of the beam energy on \(z\) is already used. Assuming that the speed of the electron is sufficiently close to the speed of light and hence its longitudinal coordinate \(z \approx ct\) and momentum \(p_z \approx \gamma mc\)

\[\frac{d\gamma}{dz} \approx \frac{eE_z}{mc^2},\]

Then it is sufficient to integrate the function \(E_z(z)\) once:

mc = 0.511 # MeV/c
E_0 = 2.000 # MeV
gamma_0 = E_0/mc + 1

from scipy import integrate

gamma = gamma_0 + integrate.cumtrapz(-Ez(z)/mc, z)

dim_gamma = hv.Dimension('gamma', label=r'$\gamma$',
                         range=(0,None))

z_gamma = hv.Curve((z,gamma), kdims=dim_z,
                   vdims=dim_gamma)

For further use in solving the envelope equation, convert the gamma array to a function:

gamma = interpolate.interp1d(
    z[1:],
    gamma,
    fill_value=(gamma[0], gamma[-1]),
    bounds_error=False
)

hv.Curve((z,gamma(z)), kdims=dim_z,
         vdims=dim_gamma)

Solving the envelope equation

The envelope equation for a circular beam of radius \(r\) with the Kapchinsky-Vladimirsky distribution with external focusing by linear fields:

\[\displaystyle r'' + \frac{1}{\beta^2\gamma} \gamma' r' + \frac{1}{2\beta^2\gamma}\gamma''r + kr - \frac{2I}{I_a (\beta\gamma)^3}\frac{1}{r} - \frac{\epsilon^2}{r^3} = 0.\]

Let’s go to the rms radius for the Kapchinsky-Vladimirsky distribution and the reweighted \(\displaystyle r_{rms} = \frac{r}{2}\) and \(\displaystyle P = \frac{2I}{I_a\beta^3\gamma^3}\). We get:

\[\displaystyle r_{rms}'' + \frac{1}{\beta^2\gamma} \gamma' r_{rms}' + \frac{1}{2\beta^2\gamma}\gamma''r_{rms} + kr_{rms} - \frac{P}{4r_{rms}} - \frac{\epsilon^2}{16{r_{rms}}^3} = 0.\]

Let \(\displaystyle y=\frac{dr_{rms}}{dz}, \displaystyle \frac{d\gamma }{dz}\approx \frac{e E_z}{m c^2}\), then:

\[\begin{split}\displaystyle \left\{\begin{matrix} \displaystyle \frac{dy}{dz}= - \frac{1}{\beta^2\gamma} \gamma' r_{rms}' - \frac{1}{2\beta^2\gamma}\gamma''r_{rms} - kr_{rms} + \frac{P}{4r_{rms}} + \frac{\epsilon^2}{16r_{rms}^3}\\ \displaystyle\frac{dr_{rms}}{dz} = y \end{matrix}\right.\end{split}\]

Let \(\vec X = \begin{bmatrix}y \\r_{rms} \\\end{bmatrix}\), now let’s make a differential equation \(X' = F(X).\)

\(k\) - solenoid stiffness:

\[k = \left ( \frac{eB_z}{2m_ec\beta\gamma} \right )^2 = \left ( \frac{e B_z}{2\beta\gamma\cdot 0.511\cdot 10^6 e \cdot \mathrm{volt}/c} \right )^2 = \left ( \frac{cB_z[\mathrm{T}]}{2\beta\gamma\cdot 0.511\cdot 10^6 \cdot \mathrm{volt}} \right )^2.\]

c = 299792458 # (m/s)

Let’s set the initial conditions:

I = 2000#2000.000 # A
Ia = 17000 # A
emitt_n = 1.80e-4#1.80e-5 # m
R_rms = 27.5e-3 #m r_rms
dR_rmsdz = 35.4e-3 #50e-3 dr_rms/dz ## 50/2*sqrt(2) = 35.4

from scipy.integrate import odeint
from scipy.misc import derivative

def dXdz(X, z):

    y     = X[0]
    sigma = X[1]

    g = gamma(z)

    dgdz = Ez(z)/mc
    d2gdz2 = derivative(lambda z:Ez(z),z, dz)/mc
    beta = np.sqrt(1-1/(g*g))
    K = ( c * Bz(z) / (2*0.511e6))**2
    emitt = emitt_n/(g*beta) # m
    P=2*I/(Ia*beta*beta*beta*g*g*g)
    dydz = P/(4*sigma) + \
           emitt*emitt/(sigma*sigma*sigma*16) - \
           K*sigma/(g*beta)**2 - \
           dgdz*y/(beta*beta*g) - \
           d2gdz2*sigma/(2*beta*beta*g)
    dsigmadz = y

    return [dydz, dsigmadz]

Solution function of the equation on the beam envelope:

def find_sigma(z):
    X0 = [dR_rmsdz, R_rms]

    sol = odeint(dXdz, X0, z, rtol = 0.001)
    sigma = sol[:, 1] # m

    return sigma # m

%opts Area.Beam [aspect=3 show_grid=True] (color='red' linewidth=1 alpha=0.3)

sigma = find_sigma(z)

dim_sig = hv.Dimension('r_{rms}', label="$r_{rms}$", unit='mm',
                       range=(0,150))

sigma_img = hv.Area((z,sigma*1e3), kdims=[dim_z],
                    vdims=[dim_sig], group='Beam')
sigma_img