Simple particle sources and coordinates

This example shows how to create particles and what are the relevant coordinates an user should be aware of.

The full example, sources_and_coordinates.py, is located at examples subdirectory.

Importing piol

Piol is imported into a Python program by

import piol

The examples add a parent path to sys.path so that the examples can be run without installing the library first.

import sys
sys.path.append('..')
import piol

Defining the reference particle

The reference particle is an imaginary particle which travels through the optical system at the optical axis. In matrix formalism all the coordinates for the reference particle are zero. This can be understood as a transfer matrix is essentially a collection of Taylor polynomes up to required order. The coordinates of the other particles are relative to this reference particle. PIOL needs to know the energy, mass and charge for the reference particle to calculate the transfer matrices and required fields in the optical elements.

PIOL uses SI units for length and time but the energies are given in eV, masses in atomic mass units and charges in multiples of the elementary charge.

See documentation for piol.referenceparticle.ReferenceParticle.

# Creating a reference particle of 50 MeV, 100 u, 18 e.
refple = piol.ReferenceParticle(energy=50e6,mass=100,charge=18)

Creating particles

Most of the simulations requires creation of particles which are then tracked through the optical system. There are multiple different source distributions available in PIOL. These are

• GeneralParticleSource,

• KVSource,

• WaterbagSource and

• GaussianSource.

The GeneralParticleSource is a source where the coordinates are totally independent on each other. The distributions can be individually chosen for each coordinates.

# Creating N=10 particles
N = 10
Source = piol.GeneralParticleSource(
    N,
    ('x',piol.Uniform(-1e-3,1e-3)),            # x from [-1,1] mm, [x] = m
    ('a',piol.Gaussian(mean=0, stddev=50e-3)), # [a] = rad     
    ('y',piol.Uniform(-2e-3,3e-3)),            # y from [-2,3] mm, [y] = m
    ('b',piol.Gaussian(mean=0, stddev=30e-3)), # [b] = rad
    ('t',0),                                   # [t] = s
    ('m',np.random.choice( np.arange(refple.mass-3, refple.mass+4 ), N ) ),
    ('K',refple.energy),('q',refple.charge) )

In PIOL the particle sources are elements as for example a dipole. This means that the when the source element is executed in the optical system, it adds the particles into the system in that moment. Therefore, to see the created particles we need also an element stack, which represents the optical system, and a driver.

Creating the optical system and a driver

The class ElementStack is a special element which which can contain other elements, also other ElementStack elements.

The elements do not contain a state. An instance of class Driver is needed to store the tracked particles.

# Creating a storage for actual elements and appending all elements to it
system = piol.ElementStack(name='example')
system.elements.append( Source )

# Creating a driver which contains all the coordinates and in general
# the state of the calculation.
driver = piol.Driver()

Running the system

All the elements have method transfer which takes driver as an argument. After the call the driver includes the coordinates of the particles in its member particles. It is a dictionary of numpy arrays, the coordinate names acting as keys. The coordinates x, a, …, q are those created by the source. There is a special coordinate id which is created by the driver and is unique integer tag for each particle.

Driver has also another member, traj, which is a list of coordinate dictionaries at each profile plane.

The following lines runs the system and prints all particle coordinates after the system.

# Running the system (only one element -- the source)
system.transfer(driver)

# Printing all coordinates after the optical system
for key in driver.particles:
    print( "{:13s} ".format(key), end='' )
print()
for index in range( len( driver.particles['id'] ) ):
    for key in driver.particles:
        print( "{:<13g} ".format( driver.particles[key][index] ), end='' )
    print()

The result can be like:

jhsaren@myyra:~/git/piol/examples$ python3 sources_and_coordinates.py
id            x             a             y             b             t             m             K             q
0             -0.000653651  -0.0151447    0.00096545    0.0450925     0             102           5e+07         18
1             0.000218771   -0.0261737    0.00237408    -0.00452792   0             101           5e+07         18
2             0.000664871   0.0778851     0.00216692    0.025384      0             98            5e+07         18
3             0.000166958   -0.0549145    -0.00125974   0.0239548     0             100           5e+07         18
4             0.000324334   -0.0116778    -0.000909095  0.0220662     0             99            5e+07         18
5             7.22529e-05   0.0750293     -0.00124698   -0.0405254    0             98            5e+07         18
6             0.00074122    -0.0139119    0.000578136   0.00708552    0             101           5e+07         18
7             -0.000855707  0.0620997     -6.78905e-06  0.00603737    0             97            5e+07         18
8             -0.000470933  0.0305785     0.00158259    -0.0185863    0             103           5e+07         18
9             -0.000774972  0.0960833     0.00141956    0.0226045     0             102           5e+07         18

The full example

# This example is documented in
# ../docs/source/example_source_and_coordinates.rst Changing the line
# numbers will mess the example documentation. So, check the
# documentation if changing anything here.

import numpy as np
# Next line can be removed if piol is installed in the system
import sys
sys.path.append('..')
import piol

# -- line 12 --
# Creating a reference particle of 50 MeV, 100 u, 18 e.
refple = piol.ReferenceParticle(energy=50e6,mass=100,charge=18)

# -- line 16 --

                                
# -- line 19 -- 
# Creating N=10 particles
N = 10
Source = piol.GeneralParticleSource(
    N,
    ('x',piol.Uniform(-1e-3,1e-3)),            # x from [-1,1] mm, [x] = m
    ('a',piol.Gaussian(mean=0, stddev=50e-3)), # [a] = rad     
    ('y',piol.Uniform(-2e-3,3e-3)),            # y from [-2,3] mm, [y] = m
    ('b',piol.Gaussian(mean=0, stddev=30e-3)), # [b] = rad
    ('t',0),                                   # [t] = s
    ('m',np.random.choice( np.arange(refple.mass-3, refple.mass+4 ), N ) ),
    ('K',refple.energy),('q',refple.charge) )
# -- line 31 --


# -- line 34 --
# Creating a storage for actual elements and appending all elements to it
system = piol.ElementStack(name='example')
system.elements.append( Source )

# Creating a driver which contains all the coordinates and in general
# the state of the calculation.
driver = piol.Driver()
# -- line 42 --


# -- line 45 --
# Running the system (only one element -- the source)
system.transfer(driver)

# Printing all coordinates after the optical system
for key in driver.particles:
    print( "{:13s} ".format(key), end='' )
print()
for index in range( len( driver.particles['id'] ) ):
    for key in driver.particles:
        print( "{:<13g} ".format( driver.particles[key][index] ), end='' )
    print()
# -- line 57 --