DIES

DIES

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:

  • Juvela M. (2024) DIES: parallel dust radiative transfer program with the immediate re-emission method, submitted

Topics:

Input files

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.

Ini-file
  • 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 .T and spectra to .spe. The default prefix is “prefix”.

  • 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).

Dust file

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.

Cloud model

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.

Example ini-file

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 

Output files

Dust temperatures

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].

Emission spectra

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.

Full example

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.