Mcrx config and output file format

Mcrx (for Monte Carlo Radiative Xfer) is the main program that does the dust scattering. It takes a grid file generated by sfrhist and produces an output file containing images of the object as well as a large amount of supplemental information.

Mcrx takes a configuration filename as its only argument. The configuration file consists of keyword-value pairs which are used to control the behavior of the program. String values are shell-expanded, so environment variables like $HOME and ~ can be used.

Important Note: No checking that required keywords are supplied is done before they are accessed. This means that the code may crash well into the run because it doesn't find some required keyword. While suboptimal, the overhead in maintaining a list of the keywords was considered worse. Also note that some keywords have default values. This is nice in that it doesn't force the user to supply all keywords, but it also means that if you mistype a keyword name, the default will silently be used. At the end of the run, a list of keywords whose values were never accessed is printed. You are strongly encouraged to check that this list doesn't contain any surprises.

Config file format

The required keywords (with example values) are, organized by category:

File parameters:

input_file grid_000.fits

The input file (created by makegrid).
output_file mcrx_000.fits

The name of the output FITS file where the results are saved. (If this file exists, mcrx tries to resume the run where it was interrupted.)
compress_images

If true, CFITSIO tile compression is used to compress the image cubes. There are issues with this because the compression is destructive. Recommend you not use it.
write_images_parallel

If true, several processes are spawned to do the tile compression in parallel.
mcrx_data_dir /home/patrik/sunrise/src/

This specifies the directory containing the interpolator file for mapping the mass of a spherical particle into a cubic cell. This file is in the /src subdirectory when you check it out with svn.

Emergence parameters

(definitions of the cameras forming images from the rays leaving the simulation volume):

cameradistance 1000 / kpc Camera distance from origin

The distance of the camera from the origin determines the perspective. The camera should be further away from the origin than any part of the simulation, otherwise strange things will happen!
camerafov 200 / kpc Full field of view of cameras

This is specifies the linear field of view of the cameras at the origin.
npixels 300 / .md Camera image dimension

The number of pixels in the images generated (always the same in both dimensions).

The positions of the cameras can be specified in several ways:

ntheta 3 / .md Number of cameras in theta direction
nphi 4 / .md Number of cameras in phi direction
exclude_south_pole false

Determines the number of viewpoints by putting cameras in a regular grid in spherical coordinates. Note that (cos theta) and phi, not theta and phi, are sampled uniformly so that the cameras are distributed uniformly in solid angle (see this URL for an explanation). The example given puts cameras in the six axes directions.
camera_position 10.0 10.0 10.0 / kpc
camera_direction
camera_up

Alternative way of specifying an arbitrary camera position explicitly.
camera_positions

Alternative way of specifying camera positions by reading a file containing lines of theta, phi in degrees. You can also specify the complete camera positions by putting 10 numbers per line, which encode position, direction, up-vector and field of view (in radians).

Cameras can also use different projection types. "Normal" cameras use a rectilinear projection, ie just projecting the points onto a tangential plane. It is also possible to create cameras to use an Hammer-Aitoff equal-area projection to make all-sky maps. (In this case, the field of view is ignored, but must still be specified.) This is controlled by setting the camera_projection keyword. It is currently not possible to have different cameras use different projections.

camera_projection rectilinear

Determines the projection used for the cameras. Valid options are "rectilinear" (default) or "aitoff".

Scatter parameters

(determines the scattering of rays):

dust_grain_file /u7/patrik/dust_data/kext_albedo_WD_MW_3.1_60

This parameter the name of the file containing the scattering characteristics of the dust. (Note that if you are also calculating dust emission, you must instead use the "grain_model" keyword in the next section.)
dust_to_gas_ratio .0 / .md M_dust/M_gas for solar Z

The constant used to convert the gas mass in the grid cells to a dust mass. This is normally 0 as we use gas-phase metals to track dust, not the total gas density.
dust_to_metal_ratio 0.4

The constant used to convert the mass of gas-phase metals to dust mass.
use_multiphase_model

If true, the Springel & Hernquist multiphase model is used to determine how much mass is in diffuse vs. clumpy phases. See JGC10 for more information.
multiphase_t0_star 1.1e9

The t0_star parameter from the S&H multiphase model.
multiphase_rho_ph 1.71e7

The threshold density parameter from the S&H multiphase model.
wavelength_file mcrx.wavelengths

.OBSOLETE The name of a text file containing the wavelengths for which the scattering should be done. Not used for polychromatic calculations, the wavelength scale is then determined by the wavelengths in the stellar SED file. Only used for v2.

Dust emission parameters:

grain_model wd01 or wd01_Brent_PAH

The name of the dust grain model used, i.e. size distribution, composition, opacities, etc. If it's set to wd01_Brent_PAH, then the template_pah_fraction keyword must be defined. Note that this is mutually exclusive with specifying dust_grain_file.
wd01_parameter_set MW3.1_60

If you are using any of the wd01 grain models, this parameter sets the particular parameter set, using the labels of WD01 and DL07. The possible labels are "MW3.1_" with one of 00, 10, 20, 30, 40, 50, or 60 appended, "LMCavg_" with one of 00, 10, or 20 appended, and "SMCbar". "MW3.1_60" is what used to be the implicit value. The DL07 models are labeled e.g. DL07_MW3.1_60. Note that you have to set the version of the PAH opacity curves you want, see next keyword.
use_dl07_opacities true

If true, the Draine & Li 2007 PAH opacity curves will be used, otherwise the Li & Draine 2001 ones. If you are using the DL07 grain models, you should set this to true in order for the model to be self-consistent.
grain_data_directory

The directory where the Draine & Lee grain opacity files are located.
grain_min_emission_wavelength

.OBSOLETE The smallest wavelength for which grain emission is calculated. Excluding wavelengths speeds up the calculation. If absent, the entire wavelength range is used. Using this is a bad idea because it screws up energy cons.
template_pah_fraction

For the wd01_Brent_PAH model, this determines the fraction of the PAH grains that are emitting according to the MAPPINGS PAH template and which fraction is emitting like thermal equilibrium grains.
ir_template_file

.OBSOLETE The name of the IR SED template file used with v 2.x of Sunrise.

Accuracy parameters:

nrays_nonscatter 1000000
nrays_scatter 1000000
nrays_aux
nrays_intensity (v3.03+)
nrays_ir

The number of Monte Carlo rays to shoot during the scatter-free, scattering, aux quantity, radiation intensity (for determining dust temperature) and final dust emission phases, respectively. (Note the meaning of nrays_ir in the presence of the integrate_ir keyword below.)
i_min 0.01
i_max 10

These keywords describe the minimum and maximum ray intensities that are maintained by the Monte Carlo algorithm. The maximum limit is enforced by the "stratification" of the forced scatterings, which sets the polychromatic reference wavelength and the maximum optical depth a ray will travel before a scattering is forced in the interval (see the paper, currently in prep). When the intensity in a ray falls below the minimum, the ray is subject to Russian Roulette and possibly terminated. (In Sunrise v3-, these were called i_last_scatter and i_split, and had slightly different functionality.)
i_forced
n_forced
n_forced_ir

.OBSOLETE The minimum intensity for rays to use forced scattering, and the maximum number of scatterings for rays to use forced scattering. These keywords are no longer used in Sunrise v4, since forcing scatterings is done automatically.
n_scatter_min 0
n_scatter_max inf

These two keywords can be used to only select radiation that has undergone a specified range of scatterings, for experimental purposes. The range is inclusive, ie setting both to 0 disables scatterings and only selects direct, attenuated light. Setting min to 1 and max unset selects only light that's scattered at least once, etc.
track_intensities

.OBSOLETE If this is false, no intensity information will be saved. This lessens memory requirements and speeds up the run if you are only interested in scattered and attenuated starlight. It defaults to true. (3.02 and earlier. The way to accomplish the same in 3.03+ is to set nrays_intensity to 0.)
intensity_subsampling

.OBSOLETE The factor by which the intensity vectors are subsampled for the purpose of calculating dust temperature. The number of wavelengths in the SED must be a multiple of this number. (3.02 and earlier.)
reference_wavelength
ir_reference_wavelength

.OBSOLETE The reference wavelength used for shooting, for the scatter and ir stages, respectively. Experience indicates that using 0.9um seems to be a good value for the scatter stage, while something around 50um works for the ir. These keywords are no longer used in v4+, since reference wavelength is selected automatically by the stratification algorithm.
ir_equilibrium_tolerance 0.1

Specifying this keyword causes mcrx to calculate the dust temperature including dust self-absorption, which requires iterating to convergence. This tolerance is enforced for all individual cells except those that fall below the luminosity percentile threshold. Note that this is a large change in meaning from its use in v3: the current meaning is that the luminosity in each cell must be converged to within this value, and the algorithm monitors the statistical uncertainty and adjusts the number of rays necessary to meet that criterion. Note that the previous meaning of the keyword applied to the integrated luminosity, so this is a much more severe criterion. A reasonable starting value is 0.1 (but even this may be overkill, do some testing in your case). This means the integrated luminosity, being the sum of all the cells, will be converged to a much higher degree than this. Also, because it will now actually measure the statistical uncertainty, setting this value lower and lower will subject you to the full consequences of the sqrt(N) Monte-Carlo convergence rate. Lowering it below .1 will very quickly increase the calculation time.
ir_luminosity_percentile 0.01

This keyword sets the percentile at which the convergence check ignores low-luminosity cells. Note the sense in which the value is defined: 0.01 will ignore the lowest-luminosity cells making up 1% of the total luminosity, and consequently include the highest-luminosity cells making up 99%. This keyword is necessary to avoid expending (a large) computational effort converging cells that have a low radiation intensity and thus poor signal/noise.
integrate_ir true

If this is true, the direct contribution to the dust emission will be determined by directly integrating the radiation transfer equation along the line of sight to each pixel in the cameras. Since this is completely deterministic, it gives noise-free images. The contribution from scattered light can't be done this way, so this will still be done based on the number of rays specified, but the scattered contribution at wavelengths where dust emission dominates is normally very small so this is quick.
use_grain_temp_lookup true

Setting this to true will make mcrx use an interpolation table for calculating the dust grain temperatures, which is a lot faster than actually doing the calculation. The default is false, but using the lookup is highly recommended for production runs. Note that this keyword was introduced in v3.02, before that it always used the interpolation table.
n_wavelengths_intensity

Determines the number of wavelengths to use in the intensity stage, when determining dust heating. This is independent of the number of wavelengths in the SED and can generally be done with much lower resolution without introducing problems. Lower number of wavelengths here will drastically lower the memory required for the dust temperature equilibrium calculation. The results with a value of 25 only differs by ~1% from that with full res, so that's the recommended value. (v3.03+)

Optional keywords

save_cell_seds false

.v4+ If true, the dust emission SEDs of all the cells will be saved to the HDU "CELLSEDS". Note that this takes n_cells*n_wavelengths*8 bytes, so for large grids and high wavelength resolutions this can make the output file significantly larger.
save_dust_temps false

.v4+ If true, the dust temperatures in the cells are saved to HDUs named "DUSTTEMPSn", where n is the dust species (currently this is always 0 but the code can handle more). The HDU is a 2D image of size n_cells*n_grains. The number of grains varies depending on the dust model, so to interpret this you need to know how the dust model works, and some grain types (for example those that just do emission based on a template) don't have a concept of temperature, so for these grains, the temperatures are always NaN. For example, using the wd01_Brent_PAH model, there are 6 different grains: First neutral and charged PAH template emitters (30 sizes each), then neutral and charged PAH thermal equilibrium emitters (30 sizes each), then graphite (81 sizes) and silicate (81 sizes) thermal equilibrium emitters. The total number of grains is thus 282, and the first 60 will always have NaN as temperature.
skip_postprocessing

If true, the run will terminate after the IR stage, and the postprocessing stages skipped.
postprocess_command_variable

If this keyword names a defined environment variable, it is used as a command which is executed instead of running the postprocessing stages.
n_threads 16

The number of threads over which the shooting and temperature calculation will be distributed. The default value is 1.
bind_threads false

If true, the threads will be bound to distinct processor cores. This will usually increase performance slightly, since it keeps cache data current.
thread_bank_size 16
Current Intel processors contain a feature called "Hyperthreading", whereby a processor core can execute more than one thread simultaneously. This makes each core appear to be two, and while it may be useful on some workloads, it is almost certain to cause a significant performance decrease in Sunrise. Ideally, hyperthreading should be disabled at the hardware level, but if that's not an option, you can use this keyword to select how mcrx binds threads to cores. The bank size works such that threads will be bound to cores 0-(banksize-1), 2*banksize-(3*banksize-1), etc. The formula is
```
core = thread_number%bank_size + thread_number/bank_size)*2*bank_size
```
Experimentation suggests that if you have a 2-socket, quad core machine, cores number 0-3 and 8-12 are real cores, while 4-7 and 13-15 are the "fake" hyperthreading cores. By setting 'n_threads=8andthread_bank_size=4`, you'll get the threads bound to the correct cores. However, this may vary on different hardware, so you'll have to experiment. If 2 threads get bound to the same core, you'll notice essentially a factor of 2 drop in performance. Note that the default value if this keyword is not specified is the value of n_threads, which makes each thread bound to the same-numbered core.
cpu_time_margin 180
wall_clock_margin 180

These parameters determine whether mcrx should try to monitor its resource limits and exit gracefully while saving state before being terminated by the system. The wall clock limit is taken from the environment variable WALL_CLOCK_LIMIT. The values are in seconds.
resubmit_command_variable RESUBMIT_COMMAND

If this parameter is set, mcrx will examine the specified environment variable and use it to resubmit itself to the batch queue as it exits because it is running out of time, as set above.
sed_file sed.txt

The name of file to which the integrated SED of the object (both non-attenuated, absorbed and attenuated in different directions) is written. The default is not to dump.
use_counters false

Specifies whether counters should print to standard out. This is nice if you are running interactively, but makes log files messy for batch jobs. Default is true.
use_hpm true

Specifies whether performance data should be collected. Default is false. This will only work if support for a performance monitoring library is enabled when compiling. Currently, the only library supported is PAPI. When this is enabled, performance data for each thread will be output to stdout after each shooting stage.
CCfits_verbose true

Specifies whether CCfits should be in verbose mode, which will print diagnostic messages if there's a problem. However, it will also print messages that look suspicious but are normal, so the default is false.
use_cuda false

If you have built Sunrise with CUDA support, this keyword will make it actually use the CUDA code. There's no harm in setting it to true if you don't have CUDA, it will just sliently fall back to using the CPU. The purpose is to make it possible to test that the CUDA calculation gives the same result as the CPU. (v3.02+)
mcrx_version v3.03

.OBSOLETE As of v3.04, this keyword is obsolete since Sunrise now saves the Mercurial version in the output files. This is a lot less error prone since it will always be correct. (Note that if you have uncommitted edits in the repository when compiling, the version string will have a "+" at the end. This of course means it's not possible to recover the exact code used, but I hope noone's doing production runs with uncommited edits...)
use_kinematics false

New in Version 4. Set to true in order to apply Doppler shifts to SEDs. This option requires that the wavelength spacing be logarithmic. Note that using very high spectral resolution over a large wavelength range will make the MCRX calculation extremely memory intensive.

Output file format

The output file is a (fairly large and complicated) FITS file. In order to propagate all the necessary input information, it contains a copy of all the HDU's in the input grid file except the "GRIDDATA" HDU.

pHDU

The pHDU has the following keywords set:

FILETYPE "TRANSFER"

The file contains the output from a radiative transfer calculation.
DATATYPE "GRID"

The data is on a grid (as opposed to a particle set).
TRANSFERCODE "MCRX"

Specifies the code which was used to do the radiative transfer calculation, as well as the HDU in which the configuration for the radiative transfer can be found. Currently, only MCRX is used here.

HDU "MCRX" (or whatever is a specified by TRANSFERCODE)

This HDU contains the keywords in the input configuration file.

HDU "LAMBDA"

This HDU is a table containing the wavelength vector used. The column "lambda" is the output wavelength vector, while "lambda_intensity" is the (shorter) logarithmic wavelength vector used for the radiation intensity/dust emission calculation.

HDU "SCATTERING_LAMBDAS"

For polychromatic (i.e., basically for v3+) runs the table is unused, the numbers of rays are saved as header keywords. When running monochromatic runs, this HDU contains the wavelengths for which the radiative transfer has been performed, as an ASCII table. The columns are:

entry long

The row in the LAMBDA HDU for the wavelengths
rlambda double

The requested wavelength (may be different from the actual wavelength used, it picks the closest entry).
normalization double

Image normalization (internal usage)
nrays long

Number of Monte Carlo rays completed

The first row, with entry = -1, is the non-scattering run. For this, only the normalization and nrays columns are used.

HDU('s) "CAMERAi-PARAMETERS"

These HDU's contain the parameters describing the set up of each of the cameras generating images of the simulation volume. They contain the following keywords:

CAMTYPE string "Pinhole"

The camera model used. "Pinhole" is the preferred model, because it has infinite focal depth. "Thin Lens" is another alternative, but is obsolete.
THETA float
PHI float

The spherical coordinates describing the position of the camera relative to the origin.
CAMERADIST float

The distance between the origin and the camera.
FOV float

The angular field of view of the camera.
LINEAR_FOV float

The linear field of view at the origin. (More user-friendly quantity than the angular fov.)
X/YSIZE int

The number of pixels in the camera image.

HDU('s) "CAMERAi-NONSCATTER"

These HDU's contain the images of the simulation volume generated without scattering. They are three-dimensional images (data cubes) with x and y dimensions given by the "npixels" configuration keyword, and z dimensions given by the number of the wavelengths in the SED's (number of rows in LAMBDA HDU and vector length of L_lambda column in GRIDDATA HDU in the input file). The pixel quantity is surface brightness as observed at the camera.

Note: In v3+, the image compression described here is typically not used, because the compression is destructive and can cause quantities to became negative. As a result, the v3 output files are commonly many gigabytes in size. Since these images are typically VERY large, the CFITSIO "image compression" convention is used. This means that every wavelength slice is compressed using the gzip algorithm and is stored as an entry in a variable-length table column. Hence, these images are NOT standard FITS image extensions, but rather binary tables. CFITSIO transparently decodes this and makes it look like an image, but software not using the CFITSIO library (e.g. IRAF) doesn't know what to do with it. This is a pain in the ass, but the compression gain is large enough that I think it's worth it.

The header contains a WCS describing the linear dimensions of the image, and also the following keywords:

IMUNIT string

The unit of the pixel values, as calculated from the SED unit along with length unit used in the grid.
SB_FACTR float

Conversion factor from internal units to surface brightness (which is what is saved).
normalization float

Image normalization (internal usage).
L_lambdai float

The luminosity of the simulation volume at wavelength specified by entry i, assuming emission is isotropic. (This is for comparison with the attenuated images.)

HDU('s) "CAMERAi-SCATTER"

These HDU's are only present when running monochromatic runs (i.e, v<=2) and contain the images of the simulation volume generated with scattering. They are three-dimensional images (data cubes) with x and y dimensions given by the "npixels" configuration keyword, and z dimensions given by the number of the wavelengths for which the radiative transfer calculations were done (number of rows in SCATTERING_LAMBDAS HDU, number of wavelengths in the file specified by the "wavelength_file" configuration keyword). The pixel quantity is surface brightness as observed at the camera. Unlike the non-scatter images, these are normal FITS image extensions. The header contains the same keywords as the non-scatter images, plus:

PIXEL_SR float

The solid angle subtended by the pixels on the "sky" of the image.

HDU('s) "CAMERAi-ATTENUATION"

These HDU's are only present when running monochromatic runs (v<=2) and contain images of the attenuation, i.e. the CAMERAi-SCATTER slice divided by the same-wavelength slice in CAMERAi-NONSCATTER. It's used when interpolating the final CAMERAi image. The header also contains the following keywords:

ATTNi float

The total attenuation at wavelength i, i.e. integrated luminosity with dust divided by without dust.

HDU('s) "CAMERAi"

These HDU's contain the images of the simulation volume including the effects of dust, with the full wavelength resolution of the "-NONSCATTER" images. For polychromatic runs, they are all calculated correctly using the polychromatic algorithm. For monochromatic runs, they are obtained by a pixel-by-pixel interpolation over wavelength of the attenuation, i.e. the "-SCATTER" image divided by the "-NONSCATTER" image. This is only valid if the direct flux dominates, which is the case most of the time in the galaxy simulations. If the pixel is dominated by scattered flux, a direct interpolation of the "-SCATTER" image is done, so in this case there will be no high-resolution features in the spectrum.

The following pertains only to Sunrise v2 (To save space, these images only contain the wavelengths in the interval spanned by the wavelengths in the SCATTERING_LAMBDAS table. Hence, the first wavelength in this image corresponds to the wavelength of the first entry in that table. These images are compressed in the same way as the "-NONSCATTER" images.)

HDU "INTEGRATED_QUANTITIES"

This HDU is created during postprocessing of the radiative transfer results. It is a binary table containing SED's integrated across the special dimensions of the images. The header also contains keywords describing bolometric quantities. The table columns are the following:

lambda double

The wavelength vector
L_lambda double

The total SED of the emission in the grid
L_lambda_absorbed double

The total SED of all the absorbed energy
L_lambda_nonscatteri double

The total SED of CAMERAi-NONSCATTER
L_lambda_scatteri double

The total SED of CAMERAi-SCATTER (v<=2 only)
L_lambdai double

The total SED of CAMERAi
L_lambda_iri double

The dust emission SED of CAMERAi
L_lambda_outi double

The sum of the L_lambdai and L_lambda_iri columns

The HDU also contains the following keywords:

L_bol_grid float

The bolometric luminosity in the grid.
L_bol_absorbed float

The bolometric luminosity absorbed in the grid.
L_nsi float

The bolometric luminosity in CAMERAi-NONSCATTER, assuming emission is isotropic.
L_scati float

The bolometric luminosity in CAMERAi-SCATTER, assuming emission is isotropic.

HDU "INTENSITY"

In v3+, this HDU contains a 2D image array of the stellar radiation intensity in the cells (number of cells by number of intensity wavelengths). This is used to calculate the dust temperature. The units are physical units of mean radiative intensity. The exact units are saved in the header keyword IMUNIT.

HDU "DUSTINTENSITY"

In v4, this HDU no longer exists, since both the stellar and dust radiation field is determined at the same time using the immediate reemission algorithm. In v3+, this HDU contains a 2D image array of the radiation intensity due to dust emission. The format is identical to that of the INTENSITY HDU, and the total radiation intensity in the cells is the sum of the two intensities.

HDU "DEPOSITION"

In v4, this HDU does not exist. In v3, this HDU is only written if "check_energy_conservation" is true, and is only used for testing purposes. It's like the INTENSITY HDU, except the quantity is energy deposition rate (which is equal to intensity divided by the mean free path to absorption). The exact units are saved in the header keyword IMUNIT.

In v<=2, this HDU contains a binary table with the deposited luminosity of the grid cells. The header also contains keywords with the quantities summed over all grid cells. The table columns are:

DEPOSITIONi double

Deposited luminosity in the cell at wavelength i
NDEPi double

Deposited luminosity in the cell at

wavelength i, derived from Ji
Ji double

Radiation intensity in the cell at

wavelength i (calculated using the path length method