PEP (or rather PEPO.py) is a program that calculated predictions for spectral line emission. Calculations are based on the escape probability formalism, and a couple of alternative forms of the escape probability formula are included. The input parameters include the kinetic temperature and the specifications of density and column density grids. The calculations for the different combinations of density and column density are performed in parallel and preferably on a GPU. With a modern GPU, the calculation of 100×100 grids points should take of the order of one second.
References:
The name of the initialisation file is give for PEP as a command line parameter. One can add at the end of the command line also another parameter “-d”, to perform all calculations in double precision. That is usually not necessary, although does not seem to be much slower either, even on GPUs.
The initialisation file is a text file, where each line starts with a keyword that may be followed by arguments. The keywords are listed below, and an example is given in example section . The lines can appear in arbitrary order, and lines starting with the comment character ‘#’ are ignored.
mol name
Molecular data are to be read from a file name.dat.
These follow the format of the files in the Leiden Atomic and Molecular
Database (LAMDA) and can be downloaded directly from there
(https://home.strw.leidenuniv.nl/~moldata/
)
tkin T
The value of kinetic temperatue T_{kin}=T for the calculations.
tlim limit
Limit included energy levels to those with energy in Kelvin units below this
threshold value, E<limit.
density min max N
Specifies a logarithmic density grid from density min to density max, with
N grid points.
colden min max N
Specifies a logarithmic column density grid from column density min to column density max,
with N grid points.
fwhm dv
Specifies the FWHM line width dv in units of km s^{-1}.
gpu flag
With flag equal to 1, asks PEP to use a GPU for the calculations.
The default (keyword not given or flag equal to zero) is to use the CPU.
cabu x1 x2 ...
Specifies relative abundances for the collisional partners (in case there are
several in the LAMDA file). Thus, the density of each collisional partner is
assumed to be the total density (as given with the keyword density) multiplied
with the corresponding x factor.
output filename
Specifies the name for the output file.
alfa val
A floating point number used to dampen the change in the escape probability
values beta over iterations. Typically not needed, but if calculation does not
converge due to oscillations, values 0<val<1 damped these.
tolerance atol rtol maxiter
Convergence criteria in terms of change per iteration for absolute values of
the level populations (atol), relative change in level populations
(rtol) and maximum number of iterations allowed (maxiter).
escape tag
Specifies the form of the escape probability formula. The argument tag can
be a string or a number: 0 = lvg, 1 = simple, 2 = slab, 3 =
sphere. These are described in the article, apart from the option simple
that has the formula beta=(1-exp(-tau))/tau.
sdevice string
Can be used to select a specific OpenCL device, based on the given string
being part of its name. One can check the device names with the clinfo
utility that is found on many software repositories.
The output from a PEP run is a single binary file that contains the computed excitation temperatures, optical depths, and line intensities (radiation temperatures). The file can be read with the following piece of Python code, which explains the exact format of the file. The name of the output file (as specified in the ini file with the keyword output) is in this example output.dat.
import numpy as np
fp = open("output.dat", 'rb')
# three 32-bit integers:
# number of column densities, densities, transitions, energy levels
NCD, NRHO, NTR, NL = np.fromfile(fp, np.int32, 4)
# NTR 32-bit floats: the frequencies of the transitions
F = np.fromfile(fp, np.float32, NTR)
# NRHO 32-bit floats: the density values
RHO = np.fromfile(fp, np.float32, NRHO)
# NCD 32-bit floats: the column density values
CD = np.fromfile(fp, np.float32, NCD)
# excitation temperatures TEX[NCD, NRHO, NTR] as np.float32 values
# the last Python-array index runs fastest in the file (row-order)
TEX = np.fromfile(fp, np.float32, NCD*NRHO*NTR).reshape(NCD, NRHO, NTR)
# optical depths TAU[NCD, NRHO, NTR] as np.float32 values
TAU = np.fromfile(fp, np.float32, NCD*NRHO*NTR).reshape(NCD, NRHO, NTR)
# radiation temperatures TRAD[NCD, NRHO, NTR] as np.float32 values
TRAD = np.fromfile(fp, np.float32, NCD*NRHO*NTR).reshape(NCD, NRHO, NTR)
# level populations NU[NCD, NRHO, NL]
Ni = np.fromfile(fp, np.float32, NCD*NRHO*NL).reshape(NCD, NRHO, NL)
fp.close()
The energy levels and the transitions are in the same order as in the input LAMDA file.
Below is an example of the PEP input file.
mol co # molecule
tkin 20.0 # kinetic temperature [K]
tlim 100.00 # energy limit [K]
density 1.0e+03 1.1e+06 100 # density grid [cm-3]
colden 1.0e+15 1.1e+18 100 # column density grid [cm-2]
fwhm 1.0 # line FWHM [km/s]
cabu 1.0 0.0 0.0 0.0 # coll. partner abundances
output demo.output # output file name
alfa 1.0 # damping (alfa=1.0 == no damping)
tolerance 1.0e-07 1.0e-05 100 # abs. and rel. tolerances, max iter.
escape lvg # form of escape prob. formula
gpu 1 # choose gpu
sdevice 4070 # choose OpenCL device by name
If the file is saved as a file demo.ini, PEP can be run from the command line as
> python PEPO.py demo.ini
This should then produce the output file my.output, which can be read as described above in the section on outputs .
The GitHub page on PEP contains an example script demo_plot.py for plotting the PEP results. In this case this produces the following figure.