Srw

SRW stands for Synchrotron Radiation Workshop.

It uses a Fourier Optics approach to compute the wavefronts.
The method is described in the publication
"Accurate and Efficient Computation of Synchrotron Radiation in the Near Field Region"
http://cds.cern.ch/record/858849

Basically, one computes first the orbit of the electron.
Then the electric field of the radiation is given by (eq. 1 in above paper)

(1)
\begin{align} \vec E(k,\vec X) = iek \int_{-\infty}^{\infty} [\vec \beta - \frac{\vec n}{R}[1+\frac{i}{kR}]]e^{ik(c\tau + R)}d\tau \end{align}

Here, $\vec \beta(\tau)$ is the velocity vector. $\vec n(\tau)$ points from the electron to the observation point. R is the distance between electron and observation point.
k is the wavenumber $k=\frac{2\pi}{\lambda} = \frac{\omega}{c}$
SRW assumes that kR is large.
Compiling and testing SRW Library and its Python binding on Linux:

1. Download and compile fftw-2.1.5 library as required for SRW (assuming that absolute path to SRW directory is "SRW_Dev"):
cd SRW_Dev/ext_lib
tar -zxvf fftw-2.1.5.tar.gz
cd fftw-2.1.5
configure —enable-float —with-pic
Manually (using editor) add -fPIC option to CFLAGS in Makefile
make
make install
cp fftw/.libs/libfftw.a SRW_Dev/ext_lib/

2. Compile SRWLib library and Python binding:
cd SRW_Dev/cpp/gcc
Make sure that Python 3.2 or higher (or Python 2.7) is installed
In Makefile, specify correct PYFLAGS path to Python include file and library
rm libsrw.a
make all
cp srwlpy.so ../../env/work/srw_python/

3. Check examples:
Make sure that path to Python 3.2 (or 2.7) is added to the PATH variable and "srw_python" to PYTHONPATH variable:
export PATH="$PATH:<absolute path to Python 3.2>" # this is not necessary if you install python using the distro's package manager export PYTHONPAH="$PYTHONPATH:SRW_Dev/env/work/srw_python/" # temporarely solution
or
echo "export PYTHONPATH=\$PYTHONPATH:SRW_Dev/env/work/srw_python/" » ~/.bashrc # permanent solution for a single user
Setting up PYTHONPATH allows to import srwlpy module from any directory.
cd SRW_Dev/env/work/srw_python
python SRWLIB_ExampleXX.py
-

from future import print_function #Python 2.7 compatibility
from srwlib import *
#import os
#import sys

#**Predefined Input Parameters:

1. as defined in glossary (Gaussian Beam)

ElectronEnergy = 6.04
ElectronCurrent = 0.2
ElectronBeamSizeH = 395e-6
ElectronBeamSizeV = 9.9e-6
ElectronBeamDivergenceH = 10.5e-6
ElectronBeamDivergenceV = 3.9e-6

1. as defined in glossary (Insertion Device)

PeriodID = 32e-3
N = 50
#case 1 E_III = 15816
#Kv = 1.452
#outFile = 'srw_id28_case1.spec'
#case 3 E_III = 21747
#Kv = 0.994
#outFile = 'srw_id28_case3.spec'
#case 4 E_V = 21747
#Kv = 1.7263
#outFile = 'srw_id28_case4.spec'
#case 5 E_V = 23750
Kv = 1.6
outFile = 'srw_id28_case5.spec'

#case 2 E_II = 15816
#Kv = 0.46
#PeriodID = 17.6e-3
#N = 91
#outFile = 'srw_id28_case2.spec'

1. as defined in glossary (driftspace) = distance from undulator to slit

d = 27.0

1. as defined in glossary (slit)

gapH = 2.4e-3
gapV = 1e-3
#others (scan, output)
PhotonEnergyMin = 4000.0
PhotonEnergyMax = 500000.0
PhotonEnergyPoints = 10000
Nmax = 91

#derived
B0 = Kv/0.934/(PeriodID*1e2)

#**End Predefined Input Parameters:

#print('SRWLIB Python Example # 6:')
print('Running SRW (SRWLIB Python)')
print('Calculating spectral flux of undulator radiation by finite-emittance electron beam collected through a finite aperture and power density distribution of this radiation (integrated over all photon energies)')

#***Undulator
harmB = SRWLMagFldH() #magnetic field harmonic
harmB.n = 1 #harmonic number
harmB.h_or_v = 'v' #magnetic field plane: horzontal ('h') or vertical ('v')
harmB.B = B0 #magnetic field amplitude [T]
und = SRWLMagFldU([harmB])
und.per = PeriodID #period length [m]
und.nPer = N #number of periods (will be rounded to integer)
magFldCnt = SRWLMagFldC([und], array('d', [0]), array('d', [0]), array('d', [0])) #Container of all magnetic field elements

#***Electron Beam
eBeam = SRWLPartBeam()
eBeam.Iavg = ElectronCurrent #average current [A]
eBeam.partStatMom1.x = 0. #initial transverse positions [m]
eBeam.partStatMom1.y = 0.
eBeam.partStatMom1.z = 0. #initial longitudinal positions (set in the middle of undulator)
eBeam.partStatMom1.xp = 0 #initial relative transverse velocities
eBeam.partStatMom1.yp = 0
eBeam.partStatMom1.gamma = ElectronEnergy/0.51099890221e-03 #relative energy
sigEperE = 0.00089 #relative RMS energy spread
sigX = ElectronBeamSizeH #horizontal RMS size of e-beam [m]
sigXp = ElectronBeamDivergenceH #horizontal RMS angular divergence [rad]
sigY = ElectronBeamSizeV #vertical RMS size of e-beam [m]
sigYp = ElectronBeamDivergenceV #vertical RMS angular divergence [rad]
#2nd order stat. moments:
eBeam.arStatMom2[0] = sigX*sigX #<(x-<x>)^2>
eBeam.arStatMom2[1] = 0 #<(x-<x>)(x'-<x'>)>
eBeam.arStatMom2[2] = sigXp*sigXp #<(x'-<x'>)^2>
eBeam.arStatMom2[3] = sigY*sigY #<(y-<y>)^2>
eBeam.arStatMom2[4] = 0 #<(y-<y>)(y'-<y'>)>
eBeam.arStatMom2[5] = sigYp*sigYp #<(y'-<y'>)^2>
eBeam.arStatMom2[10] = sigEperE*sigEperE #<(E-<E>)^2>/<E>^2

#***Precision Parameters
arPrecF = [0]*5 #for spectral flux vs photon energy
arPrecF[0] = 1 #initial UR harmonic to take into account
arPrecF[1] = Nmax #final UR harmonic to take into account
arPrecF[2] = 1.5 #longitudinal integration precision parameter
arPrecF[3] = 1.5 #azimuthal integration precision parameter
arPrecF[4] = 1 #calculate flux (1) or flux per unit surface (2)

arPrecP = [0]*5 #for power density
arPrecP[0] = 1.5 #precision factor
arPrecP[1] = 1 #power density computation method (1- "near field", 2- "far field")
arPrecP[2] = 0 #initial longitudinal position (effective if arPrecP[2] < arPrecP[3])
arPrecP[3] = 0 #final longitudinal position (effective if arPrecP[2] < arPrecP[3])
arPrecP[4] = 20000 #number of points for (intermediate) trajectory calculation

#***UR Stokes Parameters (mesh) for Spectral Flux
stkF = SRWLStokes() #for spectral flux vs photon energy
#srio stkF.allocate(10000, 1, 1) #numbers of points vs photon energy, horizontal and vertical positions
stkF.allocate(PhotonEnergyPoints, 1, 1) #numbers of points vs photon energy, horizontal and vertical positions
stkF.mesh.zStart = d #longitudinal position [m] at which UR has to be calculated
stkF.mesh.eStart = PhotonEnergyMin #initial photon energy [eV]
stkF.mesh.eFin = PhotonEnergyMax #final photon energy [eV]
stkF.mesh.xStart = -gapH/2 #initial horizontal position [m]
stkF.mesh.xFin = gapH/2 #final horizontal position [m]
stkF.mesh.yStart = -gapV/2 #initial vertical position [m]
stkF.mesh.yFin = gapV/2 #final vertical position [m]

stkP = SRWLStokes() #for power density
stkP.allocate(1, 101, 101) #numbers of points vs horizontal and vertical positions (photon energy is not taken into account)
stkP.mesh.zStart = d #longitudinal position [m] at which power density has to be calculated
stkP.mesh.xStart = -gapH/2 #initial horizontal position [m]
stkP.mesh.xFin = gapH/2 #final horizontal position [m]
stkP.mesh.yStart = -gapV/2 #initial vertical position [m]
stkP.mesh.yFin = gapV/2 #final vertical position [m]

#sys.exit(0)

#**Calculation (SRWLIB function calls)
print(' Performing Spectral Flux (Stokes parameters) calculation … ', end='')
srwl.CalcStokesUR(stkF, eBeam, und, arPrecF)
print('done')

partTraj = SRWLPrtTrj() #defining auxiliary trajectory structure
partTraj.partInitCond = eBeam.partStatMom1
partTraj.allocate(20001)
partTraj.ctStart = -1.6 #Start Time for the calculation
partTraj.ctEnd = 1.6
partTraj = srwl.CalcPartTraj(partTraj, magFldCnt, [1])

#print(' Performing Power Density calculation (from trajectory) … ', end='')
#srwl.CalcPowDenSR(stkP, eBeam, partTraj, 0, arPrecP)
#print('done')

print(' Performing Power Density calculation (from field) … ', end='')
srwl.CalcPowDenSR(stkP, eBeam, 0, magFldCnt, arPrecP)
print('done')

#**Saving results

1. open output file

fs = open(outFile, 'wb')

1. dump inputs

header = header+ '#U ElectronCurrent = ' + repr(ElectronCurrent ) + '\n'
header = header+ '#U ElectronBeamSizeH = ' + repr(ElectronBeamSizeH ) + '\n'
header = header+ '#U ElectronBeamSizeV = ' + repr(ElectronBeamSizeV ) + '\n'
header = header+ '#U ElectronBeamDivergenceH = ' + repr(ElectronBeamDivergenceH ) + '\n'
header = header+ '#U ElectronBeamDivergenceV = ' + repr(ElectronBeamDivergenceV ) + '\n'
header = header+ '#U PeriodID = ' + repr(PeriodID ) + '\n'
header = header+ '#U N = ' + repr(N ) + '\n'
header = header+ '#U Kv = ' + repr(Kv ) + '\n'
header = header+ '#U PhotonEnergyMin = ' + repr(PhotonEnergyMin ) + '\n'
header = header+ '#U PhotonEnergyMax = ' + repr(PhotonEnergyMax ) + '\n'
header = header+ '#U PhotonEnergyPoints = ' + repr(PhotonEnergyPoints ) + '\n'
header = header+ '#U d = ' + repr(d ) + '\n'
header = header+ '#U gapH = ' + repr(gapH ) + '\n'
header = header+ '#U gapV = ' + repr(gapV ) + '\n'
header = header+ '#U B0 = ' + repr(B0 ) + '\n'

#

1. write flux

#
header="#N 2 \n#L PhotonEnergy[eV] Flux[phot/sec/0.1%bw] \n"
for i in range(stkF.mesh.ne):
ener = stkF.mesh.eStart+i*(stkF.mesh.eFin-stkF.mesh.eStart)/(stkF.mesh.ne-1)
fs.write(' ' + repr(ener) + ' ' + repr(stkF.arS[i]) + '\n')

#

1. write power density (mesh scan)

#
header="#N 3 \n#L V[mm] H[mm] PowerDensity[W/mm^2] \n"

for i in range(stkP.mesh.nx):
for j in range(stkP.mesh.ny):
xx = stkP.mesh.xStart + i*(stkP.mesh.xFin-stkP.mesh.xStart)/(stkP.mesh.nx-1)
yy = stkP.mesh.yStart + j*(stkP.mesh.yFin-stkP.mesh.yStart)/(stkP.mesh.ny-1)
ij = i*stkP.mesh.nx + j
fs.write(repr(xx) + ' ' + repr(yy) + ' ' + repr(stkP.arS[ij]) + '\n')

#

1. write trajectory

#