DIES is a Monte Carlo program for one-dimensional, spherical models, to calculate equilibrium dust temperatures and the resulting dust emission spectra. DIES uses the immediate re-emission method, and cannot be used for problems where the stochastic heating of the grains is important. The Python program DIES.py assumes constant dust properties througout the model, but the alternative script DIES_MD.py allows dust properties to be set cell by cell.
References:
The initialisation file (‘ini-file’) is a text file where each line contains a key word, which is followed by one ore more arguments. Comment lines start with the # character.
background <filename>
Specifies a text file that contains the intensity spectrum of the external
radiation field. Each lines should contain two numbers, the frequency [Hz] and
the intensity value [erg/s/cm2/Hz/sr], with lines in the order of increasing
frequency.
batch <int>
Used to overcome rounding errors in runs with very high number of photon
packages. For integer arguments int>1, the run is split into that many batches.
bgpackets <int>
Specifies the number of photon packages to be generated for the background
radiation field.
cloud <filename>
Specifies a text file with the cloud model (discretised density vs. radius),
as described in the section on cloud model
. This is an obligatory
keyword.
density <float>
Optional keyword that is used to scale the density values read from the cloud file
with the constant float factor.
dust <filename>
DIES.py: Specifies the file with the dust optical properties. The format is the same as
for the CRT and SOC programs and is described in the separate section on dust
.
This is an obligatory keyword.
DIES_MD.py: The argument is the name of a text file, which lists (with
one entry per line) the names of the actual dust files (again, in the same
format as for the CRT and SOC programmes). Therefore, the number of
lines in that text file should be the same as the number of cells in the model.
forced <int>
With argument int>0, instructs DIES to use the method of “forced first
interaction” (i.e. each photon package is guaranteed to scatter or be absorbed
at least once).
gpu <int>
With gpu=0 the run is performed on the CPU, with gpu>0 on a GPU.
gridlength <float>
Optional keyword used to sets the value of the outer radius [pc] of the model cloud,
overriding the values read from the cloud file (see keyword cloud).
offsets <int>
Specifies the number of impact parameters (from the model centre to model
surface) for which emission spectra are calculated. The default value is 10.
platform <int>
Selects a specific OpenCL platform for the calculations. By default, platforms
are tested in order, until the first matching the gpu parameter is found. The
keyword sdevice may be a preferable way to select a specific device.
pointsource <filename> <k> <R>
Specifies a radiation source at the centre of the model. The first argument is the
name of a file with the spectral luminosity (text file with two columns,
frequency [Hz] and spectral luminosity [erg/s/Hz]). The second argument is a
scaling factor that is applied to the luminosity values read from the file, and the
third argument is the source radius [pc] (i.e. actually not a point
source).
prefix <prefix>
Specifies a prefix for the output files. Temperatures will be automatically
saved to
pspackets <int>
Specifies the number of photon packages to be generated from the source at the
centre of the model.
sdevice <string>
Specifies the selected OpenCL device by using a string that appears somewhere
in the name of the device (names can be checked for example by using the
clinfo program that can be found on GitHub or in many software repositories).
seed <float>
Seed for the random generators. With float>0, run is always repeated with the
same random numbers. With float<0, a random seed is generated for each run.
weightbg <float>
Weighting of the initial directions of photon packages from the external
radiation field. Value float=0.5 is the normal (unweighted case) and smaller
values cause more photon packages to be sent towards the model centre.
weightemit <float>
Weighting of the directions of reemitted photon packages. The weight function
is a Henyey-Greenstein function with the float argument being the asymmetry
parameter g. The directions are generated with respect to the direction
towards the model centre (i.e. float=g>0 directing emission preferentially
towards the model centre).
weightemitps <float>
Same as weightemit but the direction measured with respect to the radial
direction away from the model centre.
weightsca <float>
Weighting of the directions of scattered photon packages. The weight function
is a Henyey-Greenstein function with float=g (the asymmetry parameter) and with
respect to the direction towards the model centre.
weightscaps <float>
Same as weightsca but with the direction measured with respect to the
radial direction away from the model centre.
weightstepabs <float>
Weighting of the free paths for absorption, modifying the default probability
exp(-tau) to exp(-float*tau). The value float=9.0 has a special meaning,
implying the use of density-dependent weighting of the free paths (see
kerner_DIES.c for details).
weightstepsca <float>
Weighting of the free paths for scattering, modifying the default probability
exp(-tau) to exp(-float*tau).
The format of the dust files is the same as for the CRT and SOC programs without stochastic heating. The dust properties are defined by a single text file. The first two lines contain nominally the number ratio f between dust grains and H atoms and the grain radius a [cm], one value per line. The remaining file has four values per line: frequency [Hz], the value of the asymmetry parameter g, the efficiency Q(abs) for absorptions, and the efficiency Q(sca) for scattering. Lines starting with the comment character ‘#’ are ignored. For DIES.py, the dust file is specified by the keyword dust, for DIES_MD.py the dust files are listed in the text file that is specified using the keyword dust.
The format has some historical reasons, but in practice dust models consist of a distribution of grain sizes, not of a single grain size. However, all that matters is that the resulting values f×pi×a^2*Q are the correct optical depths (for absorptions and scatterings separately) per Hydrogen atom. This assuming that also the densities in the cloud file are correspondingly given as the number density of Hydrogen atoms.
The density distribution of the one-dimensional model is described with a text file. The first line contains just the number of cells. The remaining lines contain each two numbers, the outer radius of the cell (shell) [pc] and the value of the number density [cm^-3]. Optionally, the size of the model can be rescaled using the keyword gridlength and the density values can be rescaled using the keyword density. The densities are conventionally specified as the number density of Hydrogen atoms, in which case the dust file (keyword dust) must specify optical depths per Hydrogen atom, as described in the section Dust file . Other conventions are possible, as long as the descriptions of the cloud (density) and the dust (opacities) are mutually consistent.
Below is an example ini file that contains every keyword but some commented out. One should start without any of the weighting options, and try those only if they are in tests found to be useful. Only the weighting option forced should be usually safe to be included, although even that is beneficial only for optically thin models.
cloud my.cloud
background isrf.dat 1.0 # intensity of isotropic background
pointsource ps.txt 1.0 4.0e-5 # source luminosity
prefix test # prefix for output files
offsets 128 # number of output spectra
dust COM.dust # dust model
bgpackets 1000000 # external field photon packages (pkg)
pspackets 1000000 # central source photon pkg
seed -1.0 # random seed
# gpu 1 # use gpu
# platform 0 # select OpenCL device by number
sdevice 3090 # -"- by string "3090" in the name
# weightbg 0.2 # weight direction of bg packages
# weightsca 0.5 # weight direction of scattered bg pkg
# weightemit 0.5 # weight direction of reemitted bg pkg
# weightemitps 0.4 # weight direction of reemitted ps pkg
# weightscaps 0.4 # weight direction of scattered ps pkg
# weightstepsca 0.6 # weight free path for scatterings
# weightstepabs 0.6 # weight free path for absorptions
# weightstepabs 9.0 # same with density-dependent weighting
Temperatures will be automatically saved to a file prefix.T, where the start of the filename is specified with the prefix keyword in the ini-file. The file is a plain text file with two columns: the outer radius of a cell [pc] and the dust temperature [K].
The emission spectra will be saved automatically to a binary file with the name prefix.spe, where the prefix is specified in the ini-file with the prefix keyword. The file starts with two 32-bit integers, the number of frequencies and the number of offsets from the model centre. This is followed by a vector of 32-bit floating point numbers that lists the frequencies. The remaining data are the intensity values [Jy/sr] for each offset and each frequency, with the frequency being the faster running index.
The file could be read for example using the following piece of Python code
fp = open('IRE.spe')
NF, OFFS = fromfile(fp, int32, 2) # number of frequencies and offsets
F = fromfile(fp, float32, NF) # vector of frequencies
RES = fromfile(fp, float32).reshape(OFFS, NF) # intensities [Jy/sr]
The file does not directly contain the values (in pc) of the offsets, but these are equidistant between the centre and the surface of the projected model cloud.
The GitHub directory PEP contains an example script that creates a simple model cloud and point source luminosities. Using existing files (included in the directory) for the intensity of an isotropic background and for the dust properties, the script also runs the DIES program and plots its results. In that example case, the outputs should look like in the following figure.