Quickstart Guide

Getting Started

Let's start by constructing a simple FODO cell Beamline, which consists of a quadrupole magnet that focuses the beam in the horizontal plane (defocuses in the vertical), followed by a drift (empty space), then a quadrupole that defocuses in the horizontal plane (focuses in the vertical), and finally another drift.

To do this, we will define four LineElements corresponding to each of these objects. The lengths of each object are specified by L (in meters), and the quadrupole strengths can be set using the property Kn1, where n means the "normal" multipole (s would be "skew") and 1 means 1st order multipole (quadrupole).

using Beamlines

qf = Quadrupole(Kn1=0.36, L=0.5)
d = Drift(L=1.2)
qd = Quadrupole(Kn1=-0.36, L=0.5)

fodo = Beamline([qf, d, qd, d])

Beamline:
 species_ref = Inferred
 p_over_q_ref = Inferred

  Index   Name   Kind         s [m]   L [m] 
  1              Quadrupole   0       0.5
  2              Drift        0.5     1.2
  3              Quadrupole   1.7     0.5
  4              Drift        2.2     1.2

As you can see, because we did not specify a reference particle species species_ref, nor a signed reference magnetic rigidity p_over_q_ref, both of those show Inferred. This means that it will infer those values from either a preceeding Beamline, or simply leave it up to a tracking code to decide. One can also specify E_ref or pc_ref instead of p_over_q_ref for convenience. Let's do that now, and specify electrons with a total reference energy of 18 GeV:

fodo = Beamline([qd, d, qd, d], species_ref=Species("electron"), E_ref=18e9)

Beamline:
 species_ref = electron
 E_ref = 1.8e10

  Index   Name   Kind         s [m]   L [m] 
  1              Quadrupole   0       0.5
  2              Drift        0.5     1.2
  3              Quadrupole   1.7     0.5
  4              Drift        2.2     1.2

Beamlines.jl uses the AtomicAndPhysicalConstants.jl package for specifying particle species, and so any species defined by that package may be provided.

The rest of the output looks ok, except for the fact that the name column is empty! This is because we didn't specify a name property for each element. It would often be convenient if we can make the variable symbols (e.g. qf, d, etc.) automatically fill in the name field for each element. We can do exactly this by wrapping the element definitions in a @elements block:

@elements begin
    qf = Quadrupole(Kn1=0.36, L=0.5)
    d = Drift(L=1.2)
    qd = Quadrupole(Kn1=-0.36, L=0.5)
end

fodo = Beamline([qf, d, qd, d], species_ref=Species("electron"), E_ref=18e9)

Beamline:
 species_ref = electron
 E_ref = 1.8e10

  Index   Name   Kind         s [m]   L [m] 
  1       qf     Quadrupole   0       0.5
  2       d      Drift        0.5     1.2
  3       qd     Quadrupole   1.7     0.5
  4       d      Drift        2.2     1.2

Much better!

Python users may use the dict-based naming function elements instead.

@elements ...;
@elements begin
  ...
end

julia> @elements qf = Quadrupole();

julia> qf.name
"qf"

julia> @elements begin
         my_drift = Drift(L=0.5);
         my_sextupole = Sextupole(Kn2L=3.2);
       end;

julia> my_drift.name
"my_drift"

julia> my_sextupole.name
"my_sextupole"

elements(eles_dict)

Multiple LineElements in (multiple) Beamlines

Continuing with the FODO cell example above, note that fodo contains two instances of the line element d. Therefore, if the length of d is changed, then both instances of d will see this new, changed length:

d.L = 2.0
println(fodo.line[2].L)
println(fodo.line[4].L)

2.0
2.0

However, both drifts in fodo are unique elements. We can check this using the === operator:

println(fodo.line[2] === fodo.line[4])
println(d === fodo.line[2])
println(d === fodo.line[4])

false
false
false

Under the hood, when an element is placed in a Beamline, a shallow copy of that element is created that points to the "parent" element, from which it inherits its parameters. So, in this above example when the "get" fodo.line[2].L is executed, the code goes to the parent element d and returns d.L. "Sets", such as fodo.line[2].L = 10, will also pass through from the child to the parent:

fodo.line[2].L = 3.0
println(d)
println(fodo.line[4].L)
println(d === fodo.line[4].L)

LineElement:
  UniversalParams
   kind            = "Drift"
   name            = "d"
   L               = 3.0
   tracking_method = SciBmadStandard(
     radiation_damping_on      = false,
     radiation_fluctuations_on = false,
     ibs_damping_on            = false,
     ibs_fluctuations_on       = false,
   )


3.0
false

The only case where a child element can have parameters different from its parent is when a given parameter group is contained within the child. For example, fodo.line[2] and fodo.line[4] both have their own instance of BeamlineParams, from which we can extract things like beamline_index, s, and s_downstream. On the other hand, the parent element d does not have a BeamlineParams.

println("beamline_index:")
println(fodo.line[2].beamline_index)
println(fodo.line[4].beamline_index)

println("s_downstream:")
println(fodo.line[2].s_downstream)
println(fodo.line[4].s_downstream)

# This will error:
try
d.beamline_index
catch err
println(err)
end

beamline_index:
2
4
s_downstream:
3.5
7.0
ErrorException("  Unable to get key beamline_index from LineElement: element is not in a Beamline. \n  If you placed this element in a Beamline, use `findchildren` to find \n  the child instances of this element in a given Beamline.\n")

The parent element can be retrieved using parent:

fodo.line[2].parent

LineElement:
  UniversalParams
   kind            = "Drift"
   name            = "d"
   L               = 3.0
   tracking_method = SciBmadStandard(
     radiation_damping_on      = false,
     radiation_fluctuations_on = false,
     ibs_damping_on            = false,
     ibs_fluctuations_on       = false,
   )

Sometimes it can be a pain to find exactly where in a beamline are the corresponding child elements. As such, one can index the Beamline directly with a LineElement to obtain a vector of all child elements:

children = fodo[d];

2-element Vector{LineElement}:
 LineElement:
  BeamlineParams          InheritParams
   beamline_index = 2      parent = LineElement:
   s              = 0.5     UniversalParams
   s_downstream   = 3.5      kind            = "Drift"
                             name            = "d"
                             L               = 3.0
                             tracking_method = SciBmadStandard(
                               radiation_damping_on      = false,
                               radiation_fluctuations_on = false,
                               ibs_damping_on            = false,
                               ibs_fluctuations_on       = false,
                             )




 LineElement:
  BeamlineParams          InheritParams
   beamline_index = 4      parent = LineElement:
   s              = 4.0     UniversalParams
   s_downstream   = 7.0      kind            = "Drift"
                             name            = "d"
                             L               = 3.0
                             tracking_method = SciBmadStandard(
                               radiation_damping_on      = false,
                               radiation_fluctuations_on = false,
                               ibs_damping_on            = false,
                               ibs_fluctuations_on       = false,
                             )

Finally, elements in a beamline allow one to "get" parameters that may only be defined when said element is in a beamline. We showed the s and s_downstream, but another example would be the unnormalized magnetic field, if the normalied magnetic field is stored as an independent variable:

ele = Quadrupole(Kn1=2, L=2)
bl = Beamline([ele], p_over_q_ref=3)
println(bl.line[1].Bn1) # Returns Kn1 * p_over_q_ref = 2 * 3

6.0

The last parameter "set" will always define what the independent variable is. So if we then set the unnormalized quadrupole strength Bn1, that will be the independent variable:

ele.Bn1 = 10
println(bl.line[1].Kn1) # Returns Bn1 / p_over_q_ref = 10 / 3

3.3333333

Now, if we then change the reference energy of the beamline, Bn1 will remain constant but Kn1 will change:

bl.p_over_q_ref = 4
println(bl.line[1].Bn1) # == 10
println(bl.line[1].Kn1) # Now equals 10 / 4

10.0
2.5

Deferred Expressions

Earlier we set qf.Kn1 = 0.36, and qd.Kn1 = -0.36. But what if we want to ensure that qd.Kn1 == -qf.Kn1 always? We can bake-in such an interdependence, common in particle accelerator parameters, using a "deferred expression" - an expression where evaluation is postponed until its result is actually needed, rather than immediately when it is defined.

To do this, let's first define a function that returns the current value of -qf.Kn1. We can do this without giving the function any explicit name using lambda/anonymous functions:

lambdafun = () -> -qf.Kn1
println("Before: ", lambdafun())
qf.Kn1 = 0.1
println("After: ", lambdafun())

Before: -0.36
After: -0.1

Here lambdafun takes no arguments (specified by the empty tuple ()) and returns -qf.Kn1. In the context of programming, lambdafun is specifically called a closure, because it "encloses" qf, and at the time of evaluation gets the Kn1 of that enclosed qf and negates its sign.

Now we just wrap this function in Beamlines's DefExpr type, and we can set any LineElement parameter to be such a deferred expression:

qd.Kn1 = DefExpr(lambdafun)
qd.Kn1

-0.1

Now if we change qf.Kn1, evaluation of qd.Kn1 will always be -qf.Kn1:

qf.Kn1 = 0.7
qd.Kn1

-0.7

Deferred expressions can also be manipulated like any other number:

a = 1
da = DefExpr(()->a)
b = 2
db = DefExpr(()->b)
dc = da + db
println(dc())
a = 4
println(dc())
dd = sin(dc)
println(dd())

3
6
-0.27941549819892586

One can really "go crazy" with deferred expressions if they want to. They can be infinitely nested, and you can write any function that the programming language allows, for example file I/O, or even control system gets/puts with a real accelerator for a digital twin.

DefExpr{T}

julia> a = 0.36;

julia> qf = Quadrupole(Kn1L=DefExpr(()->a)); # captures variable a

julia> qf.Kn1L
0.36

julia> a = 0.7
0.7

julia> qf.Kn1L
0.7

julia> b::Float64 = 0.2;

julia> qf = Quadrupole(Kn1L=DefExpr(()->b)); # captures variable b with known type

julia> qf.Kn1L
0.2

julia> b = 0.4
0.4

julia> qf.Kn1L
0.4

julia> c = 64;

julia> cd = DefExpr(()->sin(c));

julia> c = pi;

julia> cd = DefExpr(()->sin(c));

julia> cd()
0.0

julia> c = pi/2;

julia> cd()
1.0

julia> dd = cd + 5;

julia> dd()
6.0

julia> d = DefExpr(c -> c.a + c.b);

julia> c1 = Context(a = 1, b = 2);

julia> d(c1)
3

Contexts

While DefExprs can wrap variables in the given scope as shown in the previous section, it can be useful and convenient to have a contained place where all control variables exist; this is the purpose of the Context. Contexts contain variables that can be optionally used when evaluating DefExprs that are defined with a single input argument of type Context. This is best shown with an example:

julia> c1 = Context(a = 1);
julia> c2 = Context(a = 2);
julia> d = DefExpr(c -> c.a); # one-argument lambda function
julia> d(c1)1
julia> d(c2)2
julia> c1.a = 3; # Can mutate the state of the variables
julia> d(c1)3

Contexts can be pushed on/popped from a global stack of contexts GLOBAL_CONTEXTS. In this case, when referencing a variable from a context, if it does not exist in that given context, then the first instance of that variable from the top of the GLOBAL_CONTEXTS stack will be used:

julia> c1 = Context(a = 1);
julia> push!(GLOBAL_CONTEXTS, c1);
julia> c2 = Context(b = 2);
julia> c2.a # `a` does not in exist in `c2`, get from `GLOBAL_CONTEXTS`1
julia> push!(GLOBAL_CONTEXTS, c2);
julia> c3 = Context();
julia> c3.a1
julia> c3.b2

If the variable does not exist anywhere, then an error will be thrown.

For DefExprs with a context input argument, a context must be provided in order to evaluate the deferred expression, or else an error will be thrown. For general, interactive use, contexts can be pushed on/popped from a global stack of contexts GLOBAL_CONTEXTS, which will be used if no context is provided:

julia> c1 = Context(a = 1);
julia> d = DefExpr(c -> c.a); # one-argument lambda function
julia> push!(GLOBAL_CONTEXTS, c1)DataStructures.Stack{Context}(Deque [Context[Context{Any}(Dict{Symbol, Any}(:a => 1, :b => 2, :c => 3)), Context{Any}(Dict{Symbol, Any}(:a => 1)), Context{Any}(Dict{Symbol, Any}(:b => 2)), Context{Any}(Dict{Symbol, Any}(:a => 1))]])
julia> d() # Finds `a` from the `GLOBAL_CONTEXTS` stack1

All Beamlines have a context property to store a context, which all containing LineElement parameters defined use when getting properties at the element-level:

julia> c1 = Context(Kn1=0.36);
julia> qf = Quadrupole(Kn1=DefExpr(c -> c.Kn1), L=0.5);
julia> bl = Beamline([qf], context=c1);
julia> bl[qf][1].Kn1 # Index the beamline with the `qf` to get all child `qf`s0.36

Context{T}

julia> c1 = Context(a = 1, b = 2, c = 3);

julia> push!(GLOBAL_CONTEXTS, c1);

julia> c2 = Context(a = 4);

julia> d = DefExpr(c -> c.a + c.b);

julia> d() # Uses c1.a and c1.b from the GLOBAL_CONTEXTS
3

julia> d(c2) # Uses c2.a, but falls-back to c1.b from GLOBAL_CONTEXTS
6

julia> push!(GLOBAL_CONTEXTS, c2);

julia> d() # Now uses c2.a and c1.b both from GLOBAL_CONTEXTS
6

julia> pop!(GLOBAL_CONTEXTS);

julia> d() 
3

julia> c1 = Context(Kn1=0.36);

julia> qf = Quadrupole(Kn1=DefExpr(c -> c.Kn1), L=0.5);

julia> bl = Beamline([qf], context=c1);

julia> bl[qf][1].Kn1
0.36

Parameters

Beamlines.jl supports a continually-growing list of parameters to define accelerator elements. To see a full list of the parameters you can set, look at the docstring for the LineElement type. This can be retrieved in a Julia session using Doc.docs(LineElement).

LineElement

  AlignmentParams              ApertureParams               BMultipoleParams
   x_offset                     x1_limit                     KnX
   y_offset                     x2_limit                     KsX
   z_offset                     y1_limit                     BnX
   x_rot                        y2_limit                     BsX
   y_rot                        aperture_shape               KnXL
   tilt                         aperture_at                  KsXL
                                aperture_shifts_with_body    BnXL
                                aperture_active              BsXL
                                                             tiltX

  BeamlineParams               BendParams                   EMultipoleParams
   beamline                     g_ref                        EnX
   beamline_index               tilt_ref                     EsX
   s                            e1                           EnXL
   s_downstream                 e2                           EsXL
   line                         edge1_int                    etiltX
   context                      edge2_int
   branch
   branch_index

  FourPotentialParams          InitialBeamlineParams        MapParams
   four_potential               species_ref                  transport_map
   four_potential_params        E_ref                        transport_map_params
   four_potential_normalized    pc_ref
                                p_over_q_ref
                                dE_ref
                                dpc_ref
                                dp_over_q_ref

  MetaParams                   PatchParams                  RFParams
   alias                        dt                           rf_frequency
   label                        dx                           harmon
   description                  dy                           phi0
                                dz                           zero_phase
                                dx_rot                       traveling_wave
                                dy_rot                       is_crabcavity
                                dz_rot

  UniversalParams
   kind
   name
   L
   tracking_method

ele = LineElement()
ele.L = 1   # Length 
ele.Kn1 = 2 # Normal quadrupole strength
ele.rf_frequency = 1e6 

For a more detailed description of each parameter, see the docstrings for individual parameter groups. These are shown in the Parameter Groups section of the documentation.

Polymorphism/Differentiability

To enable full auto-differentiability of all accelerator parameters, Beamlines.jl is fully polymorphic. Full polymorphism this means that you can set any parameter to be any type that you want. For auto-differentiability, a special number type that propagates the partial derivative(s) with actual value must be used in-place of the regular 64-bit floating point numbers – polymorphism allows that. Differentiable codes in accelerator physics are not new: an early example of such a differentiable code is the Polymorphic Tracking Code (PTC), written in Fortran back in the 90s.

As an example, let's see how to compute the derivative of the total length of the beamline w.r.t. a particular element length. We will use the GTPSA.jl package to do so.

using GTPSA
d1 = Descriptor(1, 1) # 1 variable, 1st order

@elements begin
    qf = Quadrupole(Kn1=0.36, L=0.5)
    d = Drift(L=1.2)
    qd = Quadrupole(Kn1=-0.36, L=0.5)
end

fodo = Beamline([qf, d, qd, d], species_ref=Species("electron"), E_ref=18e9);
println(fodo.line[end].s_downstream)

3.4000000000000004

Now we just need to update L to be a differential-algebra variable,

ΔL = vars(d1)[1] # get the first differential

d.L += ΔL
println(fodo.line[end].s_downstream)

GTPSA.TPS64{GTPSA.Dynamic}:
Descriptor(NV=1, MO=1)
 COEFFICIENT             ORDER   EXPONENTS
 3.4000000000000004E+00    0     0
 2.0000000000000000E+00    1     1

This output shows that the total length of fodo is equal to $3.4 + 2L$, which is exactly what we'd expect given that there are two drifts.

Here we just showed the length, but **any*** accelerator parameters defined in Beamlines.jl can be set to any number type (fully polymorphic), so that derivatives can be computed using any automatic-differentiation package in Julia.

After computing derivatives, e.g. during an optimization, one might want to restore all number types back to their primitive values ( Float64, Float32, etc). This can be done using the scalarize! function:

scalarize!(fodo)
println(fodo.line[end].s_downstream)

3.4000000000000004

scalarize!(ele::LineElement)

scalarize!(bl::Beamline)

scalarize!(branch::Branch)