PEP

PEP

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:

  • Juvela M.: PEP: parallel estimation of spectral line parameters with the escape-probability method, 2024, submitted

Topics:

Initialisation file

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 Tkin=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.

Outputs

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.

Example

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.