mocassin(1)

📦 version 2.02.73.2, in debian bookworm, debian trixie, ubuntu jammy, ubuntu noble, ubuntu resolute, ubuntu stonking

mocassin 💫 (section 1 of the manual, bestie — mocassin man page, we're so back)

🌸 NAME 🌸

mocassin — Monte Carlo Simulations of Ionised Nebulae ✨

🦋 SYNOPSIS 🦋

mocassin, mocassinWarm, mocassinOutput, mocassinPlot

okay so yes, that's literally four programs in one little synopsis line. it's giving quadruple threat 💅

👑 DESCRIPTION 👑

so mocassin is a fully 3D or 2D photoionisation and dust radiative transfer code, and the way it works is it uses a Monte Carlo approach to move radiation through media of arbitrary geometry and density distribution. arbitrary!! no neat little spheres required, she's flexible like that 🧚

it was originally developed for modelling photoionised regions — think HII regions and planetary nebulae 🌌 — and since then it has expanded and been applied to a whole bunch of astrophysical problems, no cap. that includes:

modelling clumpy dusty supernova envelopes 💥
star forming galaxies
protoplanetary disks
and inner shell fluorescence emission in the photospheres of stars and disk atmospheres

and the code can deal with arbitrary Cartesian grids of variable resolution. it has successfully been used to model complex density fields from SPH calculations, and it can deal with ionising radiation extending from the Lyman edge to the X-ray. from Lyman edge to X-ray. the range is unhinged and I respect it 😭🔥

also, and this part matters so much: the dust and gas microphysics is fully coupled — both in the radiation transfer and in the thermal balance. both!! she's not doing half the job and calling it a day 🫶

🎀 the four icons themselves 🎀

mocassin — this is the main MOCASSIN driver, the one you run to start a new simulation ✨
mocassinWarm — this one resumes an interrupted simulation 💖 (your sim fell over? we're back, bestie)
mocassinOutput — runs the output routines using the current grid files in the output subdirectory 📂
mocassinPlot — uses the current grid files in the output/ subdirectory to create 3d-emission maps 🌈🗺️

fr fr, that's the whole squad. four binaries, one vibe, zero half-measures 💅

🌈 How to run the benchmark cases 🌈

ok bestie, this is the part where we actually run stuff. buckle up 🎀

✨ Pure photoionisation benchmarks ✨

💖 The benchmark problems 💖

so once you've downloaded, unpacked and compiled mocassin successfully — slay — you can first of all try to run one or more of the available benchmark problems. these were designed at a series of workshops on photoionised plasma codes held at Meudon (France) and Lexington (USA). 💅

the latest version of these benchmarks is included in Pequignot et al. (2001) in Spectroscopic Challenges of Photoionised Plasmas, ASP Conference Series Vol 247; G. Ferland and D.W. Savina eds (http://adsabs.harvard.edu/abs/2001ASPC..247.....F). should those conference proceedings not be available to you, no stress — the same benchmarks are also reported in Ercolano et al., 2003, MNRAS, 340, 1136 (http://cdsads.u-strasbg.fr/abs/2003MNRAS.340.1136E). that's also where you'll find some more general info about mocassin. (there's also a link to this paper on the publications page, http://mocassin.nebulousresearch.org/publications.) 📚

ok real talk: the mocassin results included in the Pequignot et al. (2001) (http://cdsads.u-strasbg.fr/abs/2001ASPC..247..533P) paper were only from a very early version of the code. the best results are those listed in the Ercolano et al., 2003, MNRAS, 340, 1136 (http://cdsads.u-strasbg.fr/abs/2003MNRAS.340.1136E) paper. periodt. 👑

🌸 Input files for benchmarks 🌸

copy the input file of the benchmark you wish to run (you should have found these in a subdirectory called benchmarks, also included in the tar ball) into the mocassin.X.XX/input subdirectory and rename it input.in or link the input.in file to the benchmark file. 💫

e.g. for the Meudon standard HII region (Tstar=40kK) (assume mocassin is in ~/mocassin.X.XX/ ) type

cp ~/mocassin.X.XX/benchmarks/gas/HII40/meudonHII40.in ~/mocassin.X.XX/input/input.in

ln -s ~/mocassin.X.XX/benchmarks/gas/HII40/meudonHII40.in ~/mocassin.X.XX/input.in

also copy the respective nebular abundance file into the same directory (no need to rename!). in the case of the Meudon standard HII region benchmark this will be abunHII40.in (you can check the name of the abundance file, which is given in the input.in file with keyword nebComposition)

cp ~/mocassin.X.XX/benchmarks/gas/HII40/abunHII40.in ~/mocassin.X.XX/input/abunHII40.in

🦋 Run a mocassin pure photoionisation benchmark 🦋

the command you need to use to actually run mocassin will depend on how your system is set up, and whether you are required to run your models through a queueing system — in which case you'll probably need some sort of shell script to do that (ask your network administrator if unsure). 🙈

if you want to run mocassin outside any queueing system then you can try one of the following commands and hope that one works 😭 if not, ask your network administrator about how to interactively run MPI programs on your system.

for the SUN:

mprun -np N ./mocassin   (where N=number of processors you wish to use)

for the IBM

./mocassin -procs N (where N=number of processors you wish to use)

for the SGI, Beowulf and single Linux PCs using the standard MPI as above:

mpirun -np N ./mocassin (where N=number of processors you wish to use)

👑 Checking that everything is going OK 👑

if the output keyword is included in the input file, at the end of each iteration mocassin will write out a number of output files in the output/ directory. some of them are written out regardless of whether the output keyword is included — these are grid files that are needed by the mocassinWarm driver (should you decide to stop the simulation and re-start it from the same iteration); see the writeGrid keyword to see how to control the output grid files. also mocassin will keep you updated on what he is doing by outputting stuff to the screen (if that annoys you, just pipe it to a file e.g. ./mocassin -procs N > log; and then you can check the log at any time you wish 💅).

after each iteration it will tell you the percentage of grid-cells converged so far, as well as the number of 'dark cells' (i.e. those cells not reached by any of the photon packets) and also will provide a summary of eventual convergence problems. generally this is the last bit of info written before the beginning of the new iteration, so, if you have decided to pipe the screen output into a file called log, then you can easily check the progress so far by typing tail -f log, or grep Summary log, or more simply, the convergence history is summarised in a file called output/summary.out. ✨

for the benchmark cases you can start having an idea of whether mocassin is working properly by waiting until it's about 60-70% converged and then checking the output/lineFlux.out file. a description of this file is given in output files. also check the output/temperature.out and the output/ionratio.out files (also described in output files). 🎀

💜 Completing the benchmark runs 💜

the benchmark input files included are set up to run a simulation of 13*13*13 grid cells (this is sufficient to reproduce the benchmark results; see Ercolano et al., 2003 (http://cdsads.u-strasbg.fr/abs/2003MNRAS.340.1136E)) and a given starting number of energy packets (e.g. this can be set by the keyword nPhotons 100000); one iteration of 100000 packet should run in a few minutes on most systems. 🕐

the benchmarks input files make use of the autoPackets option, which automatically increases the number of energy packets to be used as soon as a convergence plateau is reached. 🔥

this could also be done manually, by excluding the autoPackets option from the input file and starting off with about 1000000 packets. however 1000000 photons are not enough to reach full convergence in this case and you will have to stop the simulation when the percentage convergence reaches 30-40%, then edit the grid3.out file changing 1000000 to 3000000, and restart the simulation using the previously compiled mocassinWarm driver, e.g.

mprun -np N ./mocassinWarm

this will basically tell mocassin to increase the number of energy packets to be used in the simulation from 1000000 to 3000000. it will, of course, now take longer to complete each iteration, but the convergence percentage should increase quite quickly. you can stop the simulation when you reach convergence > 90%. for simplicity, until you familiarise yourself with the program, it is advised that you keep the autoPackets option as set in the input files, in which case it will not be needed to manually stop and restart the simulation. 🫶

💀 Comparing with the benchmark tables 💀

the description of the output files is given in output files. the output files contain all the info (and much more) you need to compare your simulation with the benchmark tables. remember that mocassin employs a statistical method, so do NOT expect to see exactly the same figures as those quoted in the benchmark paper, as some of the differences will be due to the statistical error and also to the convergence level/limit employed. furthermore some of the atomic data has been updated since the publication of the paper. 😭

🌈 Pure Dust Benchmarks 🌈

Version 2.0 was designed for the modelling of regions where dust grains and photoionised gas coexist in the same volume; however extreme care has been taken so that the code could also be run efficiently when only one or the other process is included. in this section we will see briefly how to run dust-only models by attempting to reproduce pure dust 1D and 2D models as shown by Ercolano, Barlow and Storey (2005, MNRAS, 362, 1038) (http://cdsads.u-strasbg.fr/abs/2005MNRAS.362.1038E). 💖

🧚 1D Benchmarks 🧚

a set of spherically symmetric benchmark models and solutions are described by Ivezic et al. (1997, MNRAS, 291, 121) (http://cdsads.u-strasbg.fr/abs/1997MNRAS.291..121I). the input files used by mocassin for some of these benchmarks are included in the directory ~/mocassin.X.XX/benchmarks/dust/1D. 🌸

copy the input file of the benchmark you wish to run into the mocassin.X.XX/input subdirectory and rename it input.in or link the input.in file to the benchmark file. then run the code (see Section 3.1.3) and finally compare the output files (SED.out, dustGrid.out, see output files) with the results published by Ercolano, Barlow and Storey (2005, MNRAS, 362, 1038) (http://cdsads.u-strasbg.fr/abs/2005MNRAS.362.1038E). ✨

🍓 2D Benchmarks 🍓

a set of 2D disk benchmark models and solutions are described by Pascucci et al. (2004, A&A, 417, 793) (http://cdsads.u-strasbg.fr/abs/2004A%26A...417..793P). sample input files used by mocassin for some of these benchmarks are included in the directory ~/mocassin.X.XX/benchmarks/dust/2D. 🦋

as for the 1D benchmarks you will have to copy the input files of the benchmark you wish to run into the mocassin.X.XX/input subdirectory and rename it input.in or link the input.in file to the benchmark file. then run the code (see Section 2.1.3) and finally compare the output files (SED.out, dustGrid.out, see Section 4) with the results published by Ercolano, Barlow and Storey (2005, MNRAS, 362, 1038) (http://cdsads.u-strasbg.fr/abs/2005MNRAS.362.1038E). 💜

🌟 Running things other than benchmark cases 🌟

once you have convinced yourself that mocassin is behaving as it should, you can run your own simulations. if you are not confident that this is the case, please contact me (http://mocassin.nebulousresearch.org/contact) before going any further. 🙈

okay bestie, the very first step of all first steps: you define your model ✨ That happens through a whole lil' family of keywords you drop into input.in. And listen — YOU SHOULD NEVER HAVE TO CHANGE ANYTHING IN THE SOURCE CODE (I hope). Like fr fr. If you find you do have to, please let me know 🫶. The rest of this section lists and describes every keyword you can put in the input file. How many you use, and what order you list them in? Entirely up to you, no rules, do your thing 💅. Miss out a required, fundamental keyword and the simulation stops and tells you exactly what you left out — king behavior honestly. But heads up: for most of these keywords, default values are defined, so make sure the default is actually what you want before you decide to leave a given keyword out 👀. Defaults for each keyword are given below, wrapped in square brackets. The fundamental keywords — the ones that MUST be in every single input file — are marked by an asterisk inside those square brackets and have no default value at all.

🌸 List of keywords 🌸

`autoPackets` real1 real2 integer3

The number of energy packets used in the next iteration gets automatically boosted by a factor of real2 whenever a convergence plateau hits — and a plateau is defined by real1, i.e. when the convergence level increase is less than the value you specified in real1. Example: autoPackets 0.10 5. 1e8 will cause the total number of energy packets to increase by a factor of 5 whenever the convergence level increase from one iteration to the other is less than 10%. This only happens when the convergence level is less than 85% and the maximum number of packets defined by integer3 (1e8 in that example) hasn't been reached yet 🚀.

Default value: .false.

`continuumCube` real1 real2

This keyword builds a 3D cube of the escaped continuum radiation, in the direction given by the inclination keyword and over all directions 🌈. If no inclination is specified in the input, the continuum cube will include packets escaping in all directions. The continuum band wavelength limits are defined by real1 and real2 (given in μm). A negative value, or just omitting this keyword entirely, results in no continuum cube being written out. periodt.

Default value: -1.

`contShape` `'string'` [optional: real]

The shape of the ionising continuum 💫. The default is 'blackbody', in which case the ionising stellar continuum is approximated by the Planck function for the stellar temperature defined by the TStellar keyword. If 'string' == 'powerlaw' then this must be followed by a real number indicating the index of the power law distribution such that Fnu ∝ nu^(-index). E.g. contShape powerlaw 1. will generate an input spectrum following a nu^(-1) distribution. Note that the final result will be normalised to the luminosity entered. Also, the power law keyword is not compatible with the LPhot keyword — it only works with the LStar keyword 🚫.

If a stellar atmosphere data file is to be used, 'string' must specify the path of the external file containing the data. For example, contShape NLTE140lg65 tells the program to look for the nLTe140lg65 data file in the current directory. The stellar atmosphere files must be in a format consisting of two columns: the first listing the wavelength points in units of [A], and the second containing the corresponding wavelength-dependent stellar Eddington fluxes in units of [erg/cm^2/s/A/sr]. A set of stellar atmosphere flux tables has been compiled by Dr T. Rauch in a mocassin-friendly format and are available from http://astro.uni-tuebingen.de/~rauch/ (but please manually remove the header from the flux tables 🫠).

Default value: blackbody

`convLimit` real

This is the convergence limit for the variation of a given parameter, in each grid cell, from one Monte Carlo iteration to the next (so 0.05 means changes of 5% maximum). For a pure-gas (or gas+dust) simulation the criterion is based on the rate of change of neutral hydrogen. For dust-only, the criterion is based on the rate of change of the dust temperatures 🌡️.

Other convergence criteria can also be used — at the moment, that would require a simple editing of some source modules. If you'd like to use a different convergence criterion, please email me (http://mocassin.nebulousresearch.org/contact) and I can do the editing for you 💌.

Default value: 0.05

`debug`

Logical switch that enables debugging mode 🐛. When this keyword is included, mocassin will calculate a number of extra quantities (see Section 5.), which will, of course, slow the process down and also require more memory.

Default value: .false.

`densityFile` `'string'`

The density structure of the nebula can be defined cell by cell using an external density file 📄. mocassin knows a density file is to be used when densityFile 'string' is included in the input file, where 'string' contains the name and path of the file where the data is stored. This file must consist of four, five or six columns:

the first three columns contain the x-, y- and z- coordinates of the grid cell in [cm]
the fourth column contains the value of the hydrogen density by number in [cm^{-3}] at that particular grid cell

The x, y and z axes don't have to be equally spaced — irregular grids are perfectly acceptable to mocassin, and the extent of each axis can vary too (as long as this is consistent with the values given in the nx, ny and nz fields). The fifth column is optional. If the multiChemistry keyword is specified, the fifth column must contain an integer number in the range [1, Ncomponents] which indicates what component this cell belongs to (so that mocassin can assign the chemical abundances for this component).

It's also possible to specify nx, ny and nz from within 'string2' — the first row of the file should then be

# nx ny nz

and the keywords can be omitted from input.in.

Default value: none

`densityLaw` real...

This keyword is usually followed by a set of parameters which are to be fed into the density law routine, included in the grid_mod module. Any density law can be specified by editing the code in the setGrid subroutine. If the nebula is homogeneous, this keyword must be omitted and the Hdensity keyword included instead. And ngl, if neither of the two keywords is included, and an external density file isn't specified with the densityFile keyword, the nebular density distribution is left undefined and the simulation is halted with an error message being produced 💀.

Default value: 0. 0. etc..

`diffuseSource` real1 real2 `'string1'` integer1 integer2

You'd use this keyword in the case of a non-central source — like the heating by radioactive decay of 56Co in supernova remnants ☢️. real1 is the total luminosity of the diffuse source in [1e36 erg/s], and real2 is the temperature of the diffuse source. string1 is the shape of the source spectrum, with the same syntax as contShape.

integer1 is the number of diffuse photons to use at the start of the model (same use as nPhotons). integer2 is the index of the grid in which the diffuse source emits, which lets you control the extent of the diffuse source. For example, a medium with clumps embedded could be represented by a mother grid (index 0) and subgrids (indices 1..n). To restrict emission to the inter-clump medium only, integer2 would be set to 0. If there are no subgrids, set integer2 to 0.

Currently, for "legacy" reasons, a value of TStellar must still be specified (i.e. TStellar=0.) but nPhotons, LStar and contShape can be removed from the input.in file.

🚨 diffuseSource should not be used at the same time as symmetricXYZ.

Default value: *

`dustFile` string1 string2

names dust data files —

string1 = grain species file
string2 = grain size info file

Default value: "none", "none"

`dustMass` real1

By default, mocassin calculates the total dust mass from the user-specified number densities, dust species and grain size distributions. If this keyword is specified, the specified number densities are scaled when read in, such that the total dust mass is equal to real1. real1 has units of [Msun] ✨.

Default value: .false.

`echo` real1 real2

If you need to account for light travel time effects — for example if the source luminosity is changing — then this keyword will cause the results in SED.out to be integrated only over grid cells lying between ellipsoids corresponding to light travel times of real1 and real2 ⏳. real1 and real2 are in light days, and real1 must be less than real2. The ellipsoids open out along the z-axis, and so you should also specify 'inclination 1 0.0 -1' in input.in when using this keyword. A review of the geometries of light echoes can be found in Sugerman, 2003, AJ, 126, 1939 (http://cdsads.u-strasbg.fr/abs/2003AJ....126.1939S), and a case of mocassin modelling including light travel time considerations is described by Wesson et al, 2010, MNRAS, 403, 474 (http://cdsads.u-strasbg.fr/abs/2010MNRAS.403..474W).

Default value: .false.

`edges` real1 real2 real3

Defines the grid edges (only to be used with automatic grids). real1 real2 and real3 are respectively the x, y and z edges in cm 📐.

Default value: -1., -1., -1. *must be given if automatic grids are used

`fillingFactor` real1

real1 can assume all values from 0. to 1. to define the gas volume (and/or dust) filling factor 🎈.

Default value: 1.

`getEquivalentTau`

Logical switch to enable equivalent optical depth calculations (see Ercolano et al. 2007 (http://cdsads.u-strasbg.fr/abs/2007MNRAS.375..753E)) — useful for diffuseSource and multiPhotoSources cases. The last column in the output file equivalentTau.out gives the source SED in the same units as SED.out, which may be useful to you as well 💖.

Default value: .false.

`getTau`

…

Logical switch to enable optical depth calculations and the output that comes with them — fair warning bestie, this can be lowkey time consuming if your grids are big 💀

Default value: .false.

Hdensity real — sets a constant hydrogen density, by number, all across the nebular region. So like, Hdensity 300 will set the hydrogen density by number to the constant value of 300 cm-3. Slay. ✨

Default value: 0.

inclination integer real_1, real_2 ..... real_{n1}, real_{n2} — this one controls the viewing angles at which the SED gets calculated as it'll show up in the SED.out file. The number of viewing angles comes first (that's the integer; in this version integer <= 2, so no more than 2 angles okay?) and then the θ and φ inclination in radians. Wanna turn off the φ dependence? Set real_2 to -1. (Fun fact: θ=0. when the line of sight coincides with the z-axis.) 🦋

Default value: 0

inputNe — says that the density distribution of the grid is in terms of electron densities instead of H density. What this does is make the code calculate, at each iteration, the values of H density from the local Ne values by taking into account the local ionisation structure. It's giving ✨smart✨.

Default value: .false.

isotropicScattering — this keyword activates the logical switch that turns off unisotropic scattering (which is implemented via the Heyney-Greenstein phase function — with <cos θ> calculated from the dielectric constants via MIE theory, iykyk).

Default value: .false.

LPhot real — the number of hydrogen-ionizing photons the source emits per unit time. In the literature this is generally called Q(H0) and it's in units of [E36 sec-1]. If you give this, then the stellar luminosity LStar gets automatically derived from it. Periodt.

Default value: * if LStar not given

LStar real — the stellar luminosity, in units of [E36 erg sec-1]. If you hand this in as an input, then the number of hydrogen-ionizing photons, Q(H0), gets automatically derived from it and from the input spectrum. Iconic behaviour tbh. 👑

Default value: * if LPhot not given

maxIterateMC integer1 real1 — integer1 is the maximum number of Monte Carlo iterations to be performed in the simulation. real2 is the minimum convergence level to be achieved before the end of a simulation. The program will however stop after integer1 iteration even if real1% of convergence has yet to be reached. Not me being stubborn but okay 💅

Default value: 30 95.

MdMg string1 real1 / string2 — dust to gas ratio by mass, fr fr. If string1 = 'constant' then it must be followed by real1, holding the value of MdMg to be applied homogeneously to all cells in the grid. If string1=file then it must be followed by string2, the name of the file defining the MdMg at each location. Note that MdMg, MdMh and Ndust are mutually exclusive.

That file has to consist of four columns, with the first three columns containing the x-, y-, and z- coordinates of the grid cell in [cm] and the fourth column containing the dust to gas mass ratio at that particular grid cell. The x, y and z axis do not have to be equally spaced — irregular grids are perfectly acceptable by mocassin, and the extent of each axis can vary too (as long as this is consistent with the values given in the nx, ny and nz fields).

It's also possible to specify nx, ny and nz from within 'string2' — the first row of the file should then be

# nx ny nz

and the keywords can be omitted from input.in. 🎀

Default value: .false. 0./'none'

MdMh string1 real1 / string2 — dust to Hydrogen ratio by mass. If string1 = 'constant' then it must be followed by real1, containing the value of MdMh to be applied homogeneously to all cells in the grid. If string1=file then it must be followed by string2, the name of the file defining the MdMh at each location. Again, note that MdMg, MdMh and Ndust are mutually exclusive.

This file must consist of four columns, with the first three columns containing the x-, y-, and z- coordinates of the grid cell in [cm] and the fourth column containing the dust to hydrogen mass ratio at the particular grid cell. The x, y and z axis do not have to be equally spaced — irregular grids are perfectly acceptable by mocassin, and the extent of each axis can vary (as long as this is consistent with the values given in the nx, ny and nz fields).

Same glow-up as above: you can specify nx, ny and nz from within 'string2' — the first row of the file should then be

# nx ny nz

and the keywords can be omitted from input.in. 💖

Default value: .false. 0./'none'

multiChemistry integer1 string(1) string(2) ..... string( integer1) — use this keyword when you're doing a chemically inhomogeneous model. The integer1 value defines the number of components to be used, and string(1), string(2) ... string(integer1) are the names of the files describing the abundances in each component. When the multiChemistry keyword is included, the density distribution MUST be specified via the densityFile keyword. The fifth column of the density file must contain the index for the abundance file describing the chemical composition at each location. No skipping this step bestie, it's mandatory. 😤

Default value: .false.

multiGrids integer1 string1 — defines a multiple grid environment. integer1 is the total number of grids to be used (mother + subgrids) and string1 is the name of the file containing the subgrid information. Go check out the Running multiple spatial grids section in this document for the whole tea. 🫶

Default value: .false. 'none'

multiPhotoSources string1 — this keyword defines multiple (or single) ionising sources, and they can be placed at literally any locations. 'string1' is the name of the file containing the central star parameters. This file must contain, in the first line, the number of sources to be included, and then each source should be specified on successive lines giving (in this order) Teff [K], L* [E36erg/s], ContShape, x-, y-, z-position. Heads up: the location of the source must be given normalised to the length of the relative axis. There's an example included in the example/ directory. 🌈
nbins integer — the total number of points to be used in the frequency mesh.

Default value: 600

Ndust string1 real1 / string2 — number density [cm-3] of dust grains. If string1='constant' then it must be followed by real1, containing the value of Ndust to be applied homogeneously to all cells in the grid. If string1=file then it must be followed by string2, the name of the file defining Ndust at each location. Note that MdMg and Ndust are mutually exclusive. This file must consist of four columns, with the first three columns containing the x-, y-, and z- coordinates of the grid cell in [cm] and the fourth column containing the number density of dust in [cm-3] at the particular grid cell. The x, y and z axis do not have to be equally spaced; irregular grids are perfectly acceptable by mocassin, and the extent of each axis can vary too (as long as this is consistent with the values given in the nx, ny and nz fields).

And yep, you guessed it: it's possible to specify nx, ny and nz from within 'string2' — the first row of the file should then be

# nx ny nz

and the keywords can be omitted from input.in. ✨

Default value: .false. 0./'none'

nebComposition 'string' — specifies the path of the nebular composition data file. If you wanna use the default 'solar' composition (defined in the file solar.dat) you can just omit this keyword. But heads up: solar.dat includes all elements with Z<30, which will result in a more memory expensive simulation. So it's advised to set to zero those elements which are not needed in the simulation. Otherwise, the nebular composition can be specified in the user-defined 'string' file to be found in the current directory. All composition input files must be in a format consisting of one column containing the abundances by number and relative to hydrogen for the first thirty elements in order of ascending atomic number. The abundances of elements which are not to be included in the simulations must be set to zero (this will automatically exclude them by flagging them out throughout the program). And if the multiChemistry keyword is specified, the nebComposition keyword must be omitted. 🙈

Default value: * if gas is present and no multiChemistry

NeStart real — initial guess for the nebular electron density.

Default value: 0.

noPhotoelectric — when this keyword is included in the input file, all procedures associated with the calculation of the grain charges and photoelectric yield get switched off, as well as gas-grain collision processes. Bye bestie 👋

Default value: .false.

nPhotons integer — the number of energy packets to be used in the Monte Carlo simulation, and it has to be specified for each model. Non-negotiable. 💀

Default value: *

nStages integer — the number of ionisation stages to be used in the model. Max allowed is currently 10. If you have the atomic data necessary and would like to use more than 10 ionisation stages, please contact me (http://mocassin.nebulousresearch.org/contact) — or, if you're confident, you can edit the data/fileNames.dat and include the new data files to the pool. Main character energy tbh. 👑

Default value: 7

nuMax real — high limit of the frequency mesh, in units of [Ryd].

Default value: 15.

nuMin real — low limit of the frequency mesh, in units of [Ryd].

Default value: 0.

nx integer — number of axial points in the x-direction. Can also be specified on the first row of the densityFile, Ndust, MdMg and MdMh files, in the form of a row containing

# nx ny nz

The keywords nx, ny and nz can then be omitted. So flexible, we love that for us. 🦋

Default value: 30

ny integer

okay so you also get to pick how many axial points you want in the y-direction, same deal as before. And lowkey you can sneak it onto the first row of your densityFile, Ndust, MdMg and MdMh files instead — like, as a row that literally goes

# nx ny nz

and then you just ✨omit✨ the nx, ny and nz keywords entirely. no cap.

Default value: 30

nz integer — the number of axial points in the z-direction, bestie. yes, z. we're going 3D. You can also slide it onto the first row of your densityFile, Ndust, MdMg and MdMh files, as a row containing

# nx ny nz

and then the nx, ny and nz keywords can be left out. bye girl.

Default value: 30

output — the moment you add this keyword, the output files get written at the end of every single iteration. If you leave it out, no output gets created at all — but ngl that's fine, because if the grid files are being created you can just recover the output files with the mocassinOutput driver. icon behaviour.

Default value: .false.

planeIonization real1 real2 — this one's for when your ionisation source is a whole plane and not a lil point source. real1 has to hold the value of the incident ionizing flux above υ = real2 [Ryd] (constant at each point on the plane), in units of [phot/s/cm2]. And if real2 = 0. then real1 is assumed to be the bolometric flux.

heads up: when you use this keyword, your density distribution must be defined through the densityFile option. no exceptions.

The ionizing plane is the x-z plane. The energy packets get reflected when they hit the y-z and the y-x planes, and they can escape from the x-z planes. (This can be easily changed to suit. Please contact me (http://mocassin.nebulousresearch.org/contact) if your work requires it to be different.)

Default value: .false.

quantumHeatGrain real1 real2 — this keyword activates the temperature spiking routines (aka quantum grain heating), and it is only valid for simulations that include dust. real1 is the limiting size of grain radii [μm] that will be allowed to spike (i.e. grains with a < real1 will spike). real2 is the minimum convergence level for the spiking to occur. Please read the notes on Grain temperature spiking procedures. 💅

Default value: 1.e-3, 99.

quantumHeatGrainParameters real1 integer1 logical1 — okay so this gives you controls over some internal parameters in the quantum grain heating procedures. Do not touch these unless you really, truly know what you're doing, bestie, I'm serious. real1 is the max temperature to be considered in the grain temperature mesh of the spiking. integer1 is the number of temperature (and enthalpy) bins considered. logical1 switches on and off the writing out of a file containing all the probability functions for all the grains in every cell.

⚠️ that resulting file can be really gigantic, so set this to .true. only for debugging purposes, and use it with care. Please read the notes on Grain temperature spiking procedures contained in this manual.

Default value: 700., 300, .false.

Rin real — inner radius of the ionised region, in units of [cm].

Default value: *

Rout real — outer radius of the ionised region in units of [cm].

Default value: 0.

recombinationLines — if this keyword is in the mix, the recombination line intensity of astrophysically abundant ions gets computed and appended to the lineFlux.out file. 🎀

Default value: .false.

resLineTransfer real — real tells you at what level of grid convergence the resonance line photons escape fractions should be calculated. You wanna include this when both dust and gas are present, fr fr.

Default value: 101.

slit real1 real2 — this makes the results in lineFlux.out, temperature.out and ionratio.out get integrated over only those cells that fall under the projection of a slit aligned along the z-axis of the grid. The slit x and y dimensions (in [cm]) are defined by real1 and real2 respectively. If real1 and real2 have value 0. or if they're just omitted, then no slit, and the results get integrated over the whole active volume. simple as.

Default value: 0., 0.

symmetricXYZ — when the nebula you're modelling shows axial symmetry in the x-, y- and z-directions, you can use this keyword to enable the symmetric grid procedures. What that does is put the ionizing source in a corner of the grid instead of the centre, which means only one eighth of the nebula has to be computed. efficiency queen 👑

not me warning you again: symmetricXYZ should not be used in models illuminated by a diffuse source.

Default value: .false.

talk — this switch enables the verbose version of the program. she talks a lot. we love her for it.

Default value: .false.

TeStart real — initial guess for the nebular temperature.

Default value: 10000.

traceHeating — a logical variable to switch on the thermal balance channel tracing. When you include this in the input file, a file called thermalBalance.out gets written to the output/ directory. Just be aware that depending on the size of your grid, this may be quite a large file and time-consuming in the I/O phase. 😭

Default value: .false.

TStellar real — temperature in [K] of the ionizing stellar source. your star's whole vibe, temperature-wise.

Default value: *

writeGrid real — real indicates the minimum grid convergence percentage after which the grid files will be written out.

Default value: 0.

📁 Input and Output Files

Input files

The source files live in a subdirectory called source/. mocassin looks for the input.in file from a subdirectory called input/. All the atomic data files should be sitting in a subdirectory named data/. Most of the atomic data files should not need to be changed at all — unless you decide to update some of them, in which case (under the GPL agreement) you should also email me (http://mocassin.nebulousresearch.org/contact) with the changes so I can include them in the public version of the code. 🫶

The dust optical data library and other dust related data files are contained in a subdirectory named dustData/.

The user's input files may be a combination of the following files, depending on the processes included in a given simulation:

input.in — this is the main input file where you get to specify all the keywords that define your simulation. Some example input files are given for the Meudon/Lexington benchmark cases (see pure photoionisation benchmarks).
gas abundances file — the nebular abundances file, which should have the name specified by the user in the nebComposition field of the input.in file. Some sample files are given for the Meudon/Lexington benchmarks.
gas density file — the nebular density structure file, which should have the name specified by the user in the densityFile field of the input.in file. The format of this file is given in Section 4 (see densityFile).
stellar atmosphere file — the stellar atmosphere file, which should have the name specified by the user in the contShape field of the input.in file. The format of this file is given in Section 4 (see contShape).
dust number density file — this should contain the dust number density distribution across the grid. Its name and path should be specified in the input.in file by the keyword Ndust.
dust to hydrogen or dust to gas ratio — this contains the dust to hydrogen or dust to gas ratio distribution across the grid. Its name and path should be specified in the input.in file by the keyword MdMh or MdMg, respectively.
dust species and grain size distribution files — the names of these two files must be specified in the input.in file following the keyword dustFile. The dust species file should contain a first line specifying how many species are to be included, and then successive lines containing the names of the optical data (n,k or Qs) file and the relative abundance of the species.

e.g. for a pure silicate model (using the Draine and Lee 1984 data (http://cdsads.u-strasbg.fr/abs/1984ApJ...285...89D)):

1
'dustData/sil-dl.nk' 1.0

The grain size distribution file should contain a first line specifying how many grain sizes are to be included, the rest of the file should consist of three columns : index, radius (in um), weight. Grain size distribution files can be created using the makeGrainSizeDistribution.f90 program included in the accessories/ subdirectory.

e.g. for a single grain size:

1 size
1 0.16 1.0

Input files for multigrid simulations are described in [[running multiple spatial grids.

plot.in — the plot.in file is used by the mocassinPlot driver in order to create 3D grids of line emission. This file must be placed in (or linked to) the input/ subdirectory. The plot.out and the grid4.out files are written out to output/ and can then be used to create emission line maps by integration along any given line of sight. 🌈

Monochromatic grids are created using the mono keyword, and individual lines using the line keyword.

For example:

mono
line 2        4861.   4861.
line 93       4686.   4686.
line 1529     5007.   5007.
line 1540     4363.   4363.
line 2407     6733.   6733.
line 2408     6718.   6718.
line 929      6583.   6583.

The integer following the keyword line is the line index number as given in the lineFlux.out file.

NOTE: the line indices will be different for different simulations, because they depend on which elements are included and on the number of ionisation stages accounted for.

The 2nd and 3rd indices contain the central wavelength of the line (these are redundant for monochromatic plots, however they must be included).

Output files

the output files, and bestie there are so many 💖

okay so mocassin produces several output files at various times during the simulation — it's literally posting updates the whole time, no cap. They all get contained in a subdirectory named mocassin.X.XX/output/. The files ionratio.out, lineFlux.out, temperature.out, (tau.out), ionDen.out and SED.out are all produced by the output_mod module. ✨

The plot.out and grid4.out files are produced by the mocassinPlot file.

`ionratio.out` 👑

ionratio.out contains the volume averaged ionic fractions. Different authors in the past have used slightly different definitions of this quantity in their models, so pls refer to Ercolano et al. (2003) (http://cdsads.u-strasbg.fr/abs/2003MNRAS.340.1136E) for further information on the description used by mocassin. iykyk, definitions are a whole thing.

The first two columns of the ionratio.out file give the atomic number of the element and its ionisation stage (1 for neutral, 2 for singly ionised etc.), and the third column gives the required quantity. If a multiChemistry model is being run, then the results will be given for each individual component. 🧚

`lineFlux.out` 🌸

The file lineFlux.out contains the volume integrated intensities of all the emission lines calculated by mocassin. These are all given relative to Hβ, which is in absolute units. slay.

The first two columns give the element and ion number, these are followed by mocassin codes for the levels of the transition; these are followed by the wavelength in [A]. The wavelength column is followed by the analytical and Monte Carlo line intensities relative to Hβ, which is given in absolute units at the top of each region. Finally the last column gives the ratio of the two previous columns.

NOTE that the Monte Carlo line intensities are only calculated if the debug keyword is included in the input file. In normal mode only the intensities calculated using the formal solution (which are in general more accurate) will be available. so yeah, turn on debug or you simply don't get them 🙈

If a multiChemistry model is being run, then the results will be given for each individual component.

`temperature.out` 💜

The file temperature.out contains the mean electronic temperatures weighted by the ionic species (see Ercolano et al., 2003 (http://cdsads.u-strasbg.fr/abs/2003MNRAS.340.1136E) for definition). This file has the same structure as the ionratio.out file. twin behaviour, honestly iconic.

If a multiChemistry model is being run, then the results will be given for each individual component.

`equivalentTau.out` 🦋

Please see Section 2.3 of Ercolano, Barlow and Sugerman (2007) (http://cdsads.u-strasbg.fr/abs/2007MNRAS.375..753E) for a description of this quantity. The first column contains Energy [Ryd], the second column contains wavelength [μm], the third column contains equivalent tau [see paper] and the fourth column contains Fλ0 [see paper]. Note that in the case of diffuse illumination this is the only meaningful quantity — see discussion in paper.

`tauNu.out` ⭐

The file tauNu.out contains the frequency dependent optical depth tau(nu) [which includes ALL opacity sources] at the edge of the grid in the three positive axial directions starting from the origin of the axes. ALL of them, bestie, every single opacity source. 💅

`tau.out` 💀

THIS FILE IS NOT PRODUCED ANY MORE IN VERSIONS >= 2.02.49 PLEASE CONTACT ME (http://mocassin.nebulousresearch.org/contact) IF YOU NEED IT.

The file tau.out contains the run of the optical depth from the centre of the nebula to the outer edge along the three axial directions. The optical depths are calculated at the neutral hydrogen ionisation threshold, nu = 1.0 Ryd (13.6 eV), at the neutral Helium ionisation threshold, nu = 1.8 Ryd, and at the singly ionised Helium ionisation threshold, nu = 4.0 Ryd. The first column of the file gives the distance in [cm] from the centre of the nebula and the second column gives the optical depth from the centre to that point. This file is not always calculated correctly if the internal switches are not set up properly. Please use equivalentTau. See description above. fr fr, just use equivalentTau.

`ionDen.out` 🍓

The file ionDen.out contains the ionic fractions at each grid cell. The first three columns give the x-, y- and z-axis indices of the cell, the fourth and fifth columns give the atomic number and the ionisation stage of the element (as above, 1 for atom, 2 for singly ionised etc.) and, finally, the sixth column gives the corresponding ionic fraction.

`SED.out` 🌈

This file contains the emerging spectral energy distribution from the grid. As indicated in the files' header, column 1 and column 2 contain the frequency [Ryd] and wavelength [μm] grid, column 3 contains the direction averaged SED per unit direction (must multiply by π to obtain the total overall directions); the following columns contain the SED emerging in any given line of sight as requested by the user in the input.in file with the inclination keyword. ✨

The grid files and `photoSource.out` 🎀

As we have already mentioned elsewhere, grid files are also written out after each iteration by routines contained in the grid_mod module. These are needed by the warm start driver (mocassinWarm) to re-initialise an interrupted simulation. These files are formatted such that they can be written out and read back in quickly and therefore they may not be very clear to the human eye. However, most of the information they contain is also given in a more intelligible form in the other output files listed above — dust temperatures are an exception as discussed in plotting dust temperatures. Gas-only simulations will result in 4 grid files being written out: grid0.out, grid1.out, grid2.out and grid3.out.

The first line of the grid0.out file gives the number of grids included (mother + subgrids), on the next line are the x-, y- and z-axes points in the mother grid, followed by the outer radius in [cm]. The next few lines list the x-axis points, then the y-axis points and, finally, the z-axis points. The rest of the file contains the convergence info for each grid cell. The active index of the cell in the first column, whether it has converged (1=yes, 0=no), and whether it is a black cell (i.e. if it could not be reached by any photons, 1=black, 0=normal). Cells that have a 0 index are inactive cells; cells with a negative index are cells that have been replaced by a subgrid, whose index is equal to the absolute value of the negative mother-cell index. In the case of multiGrid simulations, the file will loop around all subgrids included. 🔁

The grid1.out file contains the electron temperatures, electron densities and hydrogen densities for each grid cell in each grid. As for the grid0.out file, this information is given for each grid cell, with the last index varying the fastest (i.e. (1,1,1), (1,1,2), etc.).

The file grid2.out contains the ionic fractions at each grid cell for the ions included in the simulation. These are given in order of increasing atomic number and ionisation stage, with each element occupying one line. The grid cell indices vary in the same fashion as in the grid1.out file.

The file grid3.out contains a list of the specified simulation parameters in a fixed order (the keywords are indicated on the right of each value).

NOTE that it is not possible any more to change nuStepSize from the input.in file. If you wish to change this parameter (you shouldn't need to), this is defined in the set_input_mod module.

Dust-only simulations only produce grid0.out, grid3.out and dustGrid.out. 🫶

The file dustGrid.out contains the dust number density at each cell followed by the grain temperatures for each grain size of each grain species. For each cell Ndust is written on one line and this is followed by n_size+1 lines of n_species+1 columns containing the individual grain temperatures for each size and species, where n_size=number of size bins and n_species=number of grain species included. The average temperature for the grain mixture weighted by the size distribution and the species abundances is given at (0,0).

Dust+gas simulations will produce all the files above. ✨

grid4.out is written out by mocassinPlot and it just contains the volume of each gridcell in [e45 cm3], which is needed for visualisation purposes.

The ionising source(s) input parameters are written out to a file named photoSource.out.

`plot.out` 🎀

The mocassinPlot driver produces also an output file, containing the luminosities of each individual grid volume element in the required emission lines. This file, named plot.out, is written in a format which should be easily readable into a data visualisation software, such as IDL or PDL. A grid4.out is also written out, containing the volume of each gridcell in [e45 cm3].

🌈 Miscellaneous notes on mocassin 🌈

Analytical and Monte Carlo line fluxes 💫

The total luminosity of the nebula in various emission lines longward of the Lyman limit can be obtained by using two methods. 👯

The first method, which is only available to Monte Carlo codes, consists of summing up the number of energy packets in the given line over all the grid cells. From this, the power emitted in that line can be readily obtained (see e.g. Ercolano et al., 2003 (http://cdsads.u-strasbg.fr/abs/2003MNRAS.340.1136E)).
The second method consists of using the values of the local electron temperatures and ionic abundances given by the final model solution to obtain the line emissivities for each grid cell. The luminosity of the nebula in any given line can then be calculated easily by integrating the emissivity of the required line over the volume of the nebula.

Both routes, same destination. periodt. 💅

A side-by-side look at the results from the two methods described above gives you a lil vibe check on the statistical accuracy you actually got out of the simulation, because the two methods only agree when enough energy packets were used to get good statistics for every single line. But bestie, lowkey, the second method (the formal solution) is the one giving the most accurate results in general, especially for weak lines — aka the ones that only emit a few photons. 🫶

Working out the line emissivities with the first method is straightforward, sure, but it needs extra book-keeping, and that gets expensive for bigger simulations. So for that exact reason, this calculation only happens when the keyword debug is in the input file. Otherwise, only the speedier (and more accurate) formal solution gets employed. Periodt.

🦋 Running multiple spatial grids 🦋

WARNING: This section is out of date — new help is in the process of being written as we speak — almost ;-) . Please email me (http://mocassin.nebulousresearch.org/contact) if you want to use multiple grids and need some help

From Version 2.00 onwards you can choose to run simulations with multiple spatial sampling. The classic example of when you'd need this is modelling knots embedded in a gaseous nebula. The density enhancement inside the knot needs a finer spatial grid than whatever is fine enough for the main nebula. In those cases the subgrid describing the knot can be modelled simultaneously and self consistently by mocassin's multigrid routines. The radiation field gets safely transferred from mother to sub grids with (hopefully) no error being introduced along the way. The overheads involved only reflect the eventual introduction of extra grid cells. 💅

A typical example of how such a model gets set up is briefly described below. All files mentioned should be in the examples/multigrid subdirectory.

The example below shows the input.in file for a plane parallel model consisting of a cubic blob embedded into a cubic grid:

autoPackets 0.20 2. 1000000000
output
contShape  blackbody
nebComposition 'examples/abunHII40.in'
maxIterateMC  30 95.
nPhotons 10000000
nx 10
ny 10
nz 10
planeIonization 3.0
Rin 0.
Rout 2.1e19
TStellar 45000.
writeGrid 5.
convLimit 0.03
densityFile 'examples/cubenew.dat'
multiGrids 2 'examples/subGrids.in'

The density distribution of the mother grid (the one holding the main nebula) is provided provided by an external file using the densityFile keyword as described in Section 3.1; e.g.

densityFile 'examples/cube.dat'

mocassin knows it's gonna have to deal with more than one grid (2 in this case) thanks to the line

multiGrids 2 'examples/subGrids.in'

The integer after the multiGrids keyword is the total number of grids including the mother grid; the file example/subGrids.in holds all the info about the position of the subgrids and their specification. This file gets created automatically by using the makeSubGrids.f90 program included in the accessories/ directory. ✨

The makeSubGrids.f90 program reads the makeSubGrids.in input file, which has the following form

1
cube.dat 10 10 10 1
cubenew.dat
1 9 9 9 knot_den.in
7.e18 14.e18 7.e18 14.e18 7.e18 14.e18
100. -1

line 1: number of subgrids to be included
line 2: name of mother grid file; its x,y,z dimensions; multiChemistry? (1=no; >1=yes)
line 3: name of the modified mother grid file (to be created by makeSubGrids.f90)
line 4: knot index; its x,y,z extent; name of knot file (to be created by makeSubGrids.f90)
line 5: xmin,xmax,ymin,ymax,zmin,zmax in cm, defining the region occupied by knot 1
line 6: H number density in [cm-3] of knot 1, abundance index (negative if no multiChemistry)

The maximum number of subgrids you can include in a simulation is controlled by the parameter maxGrids in the source/constants_mod.f90 module; in the example above only one subgrid is included, and you could add more by simply repeating lines 4,5 and 6 for each subgrid. 🌸

The knot_den.in file above has the following format

0.0000000E+00  0.0000000E+00  0.0000000E+00   100.0000
0.0000000E+00  0.0000000E+00  0.1250000       100.0000
0.0000000E+00  0.0000000E+00  0.2500000       100.0000
0.0000000E+00  0.0000000E+00  0.3750000       100.0000
0.0000000E+00  0.0000000E+00  0.5000000       100.0000
0.0000000E+00  0.0000000E+00  0.6250000       100.0000
0.0000000E+00  0.0000000E+00  0.7500000       100.0000
0.0000000E+00  0.0000000E+00  0.8750000       100.0000
0.0000000E+00  0.0000000E+00   1.000000       100.0000
0.0000000E+00  0.1250000      0.0000000E+00   100.0000
0.0000000E+00  0.1250000      0.1250000       100.0000
etc ........

Columns 1, 2 and 3 contain the x, y, and z positions in the subgrid, normalised to 1. The fourth column contains the h number density. (when running a multiChemistry model a fifth column would appear; containing the abundance file index).

It's worth noticing that since the x y and z positions in the knot_den.in files are given normalised to unity, we can include knots with the same geometry and density at several grid locations without needing to create a density file for each knot. This is done by specifying the same filename for the knots in the makeSubgrids.in file. 💖

The makeSubgrids.f90 file will also process the cube.dat file (which is the original mother grid density file provided by the user) turning off all the cells that are to be replaced by subgrids. The new file will be called cubenew.dat in the above example, and this is the file that must be specified in the input.in file of mocassin.

The makeSubgrids.f90 file can only create subgrids of homogeneous density and chemical indices. It should be very easy for you to customise the grids you get to include your own 'favourite' density/chemistry distribution. 🎀

The output from multigrid models gets summarised in the usual files, and you get an overall result as well as a grid by grid result.

Christophe Morisset of UNAM is currently in the process of creating an IDL visualisation tool that can deal with multiGrids and reconstruction, for a quick and efficient analysis of the results. The plan is to distribute the new tool as soon as it becomes available. 👑

🧚 Running models that include dust and gas. 🧚

The gas and dust radiative transfer routines in mocassin are now fully integrated. So it is now possible to run mocassin in its original gas-only mode, as well as dust-only and of course dust+gas. The basics on how to run the code for dust-only or gas-only cases have already been given in pure photoionisation benchmarks and pure dust benchmarks; here we're gonna concentrate on an example where both dust and gas are present. 💫

Below is the input.in file used for the dust and gas model of NGC 3918 as described by Ercolano, Barlow and Storey (2005, MNRAS, 362, 1038) (http://cdsads.u-strasbg.fr/abs/2005MNRAS.362.1038E).

autoPackets 20. 2. 500000000
densityFile 'ngc3918/ngc3918den.dat'
symmetricXYZ
multiChemistry 2 'ngc3918/ngc3918.dat' 'ngc3918/ngc3918.dat
contShape 'ngc3918/nLTe140lg65'
maxIterateMC  30 95.
nPhotons 1000000
nx 23
ny 23
nz 23
nbins 700
LStar 27.64
nuMax 23.7
nuMin 3.1e-4
Rin 0.
Rout 3.27142e17
TStellar 140000.
MdMh constant 0.0011
dustFile 'ngc3918/primary_grainspecies.dat'
'ngc3918/primary_grainsizes.dat'
writeGrid 50.
convLimit 0.03
resLinesTransfer 90.

In summary, it should be clear from the example above that the only keywords that differentiate the file above from the input of a gas-only model are:

MdMh (dust to hydrogen mass ratio)
dustFile (files containing the grain size distribution and species information)
resLinesTransfer (which tells at what level of grid convergence the resonance line photons escape fractions should be calculated).

Note that the first two keywords that define how much, what type and what size grains are to be used are also needed to run dust-only models (although MdMh could be substituted by Ndust or by MdMg). resLinesTransfer is only needed if there is gas also in the simulation (which would then be emitting resonance lines capable of heating the grains). 🔥

💜 The accessories/ subdirectory 💜

A number of useful (or not) FORTRAN and IDL programs are included in the accessories/ subdirectory. Some words of warning: the programs are very basic and poorly commented, as they were developed for personal use. Anyone is welcome to use them at their own risk!! The IDL programs, in particular, are specific to the simulation they were developed for, and some editing will be necessary to customise them to the user's needs. 😭

🌈 Plotting dust temperatures 🌈

Depending on how complex the dust model employed in a given model is, and how geometrically complex the grid is, it may be quite challenging to explore the dust temperature distribution calculated by mocassin and written out to dustGrid.out. 😵‍💫

Sometimes the only way is to plot out the results or create 3D maps using a visualisation program such as IDL, PDL etc..

The accessories/ subroutine ships with an example — dustTemperatures2.pro — showing how you'd visualise a grid like that using IDL. It was written for a grid holding 2 grain species and 20 grain sizes, and because that particular simulation was a multiChemistry dust and gas one, the results are also split up by sector. 🎀

🌡️ Grain temperature spiking routines 💥

The grain temperature spiking routines in mocassin are built on the Guhathakurta & Draine (1989 ApJ 345, 230) method (http://cdsads.u-strasbg.fr/abs/1989ApJ...345..230G). Bestie, this method is so powerful — it lets you compute the time-dependent grain temperatures from quantum heating without it taking forever. ✨

For the limitations of the method, also go read Siebenmorgen et al. (1992 A&A 266, 501) (http://cdsads.u-strasbg.fr/abs/1992A%26A...266..501S). 🫶

The temperature spikes only affect the output SED from dust grains. So in the case of gas+grain simulations it's advisable not to switch these time-consuming procedures on until the model has almost converged — keep a high value of real2, see quantumHeatGrain. 😌 In gas-only or dust-only models though, it's worth having the procedures working right from the start, because the convergence criterion is then based on dust temperatures — so you literally have to take this effect into account from the very beginning or you'll get convergence fluctuations. 🙈 And ngl, it's in general not worth running the quantum heating routines on large grains that are unlikely to spike. For a full discussion of when quantum heating routines actually must be considered, again: Siebenmorgen et al. (1992 A&A 266, 501) (http://cdsads.u-strasbg.fr/abs/1992A%26A...266..501S).

At present the grain temperature spiking routines are only implemented for carbonaceous or silicate grains. mocassin will expect to be told what type of grain he's dealing with when calculating the spiking, and that's done by adding a -capital- 'S' or 'C' at the beginning of the species label in the optical constants file. E.g. for the Draine & Laor (1993) (http://cdsads.u-strasbg.fr/abs/1993ApJ...402..441L) silicates data (dustData/sil-dlaor.nk):

nk
Ssil_dl  1400. 3.3  0.588 20.077
0.10000E-02 0.99956E+00 0.97380E-04
0.10120E-02 0.99955E+00 0.10160E-03
0.10230E-02 0.99954E+00 0.10610E-03
0.10350E-02 0.99953E+00 0.11060E-03
0.10470E-02 0.99952E+00 0.11520E-03
0.10590E-02 0.99951E+00 0.12000E-03
.........

Please email me (http://mocassin.nebulousresearch.org/contact) if you'd like to include T-spiking effects for other species. 💌

🚧 Limitations and future development 🛠️

mocassin's principal limitation is imposed by the computer power available. Periodt. 💅 The great volume of data that has to be handled in a three-dimensional simulation implies you need a system with multi-processing capabilities to accelerate the computational time. But the fast development of Beowulf Linux clusters is making parallel computing more affordable — and that's also a reason why the MPI formalism was chosen (as opposed to openMP), because it allows information to be passed from one processor to another and therefore it doesn't necessarily require a system with shared memory facilities. Such systems, which include the Silicon Graphics Origin 2000 machine used for this work, are generally much more expensive than Beowulf clusters. 😭

Pure-dust models are less computationally expensive than gas ones, or dust-and-gas ones — and reasonably sized grids can feasibly be run on single processor machines. 👑

mocassin was designed for modelling the photoionised region of planetary nebulae and H II regions, and it does not, at present, include the high energy physical processes needed for stuff like modelling AGNs. However the inclusion of processes such as inner shell photoionisation and Compton heating is straightforward, and that's intended to be one of the developments of the near future. ✨ Moreover, current efforts are geared toward including a self-consistent treatment for PDR processes — and that requires incorporating a chemical network. 🧪

🐛 BUGS 🐛

No known bugs. iconic behaviour, honestly 🌈

👑 AUTHOR 👑

Barbara Ercolano

And that's a wrap, bestie — this is mocassin version 2.02.70, dated 31 Dec 2015, from the mocassin(1) manual. Go forth and model some nebulae 💖✨🦋

📦 where this page comes from

debian bookworm: mocassin 2.02.73.2-1 · GPL-2.0-or-later
debian trixie: mocassin 2.02.73.2-1 · GPL-2.0-or-later
ubuntu jammy: mocassin 2.02.73.2-1 · GPL-2.0-or-later
ubuntu noble: mocassin 2.02.73.2-1 · GPL-2.0-or-later
ubuntu resolute: mocassin 2.02.73.2-1build1 · GPL-2.0-or-later
ubuntu stonking: mocassin 2.02.73.2-1build1 · GPL-2.0-or-later
💖 girlified by deepseek-ai/DeepSeek-V4.1-Flash on 2026-10-13: a language model's rewrite of the original, under its terms

mocassin(1)

mocassin 💫 (section 1 of the manual, bestie — mocassin man page, we're so back)

🌸 NAME 🌸

🦋 SYNOPSIS 🦋

👑 DESCRIPTION 👑

🎀 the four icons themselves 🎀

🌈 How to run the benchmark cases 🌈

✨ Pure photoionisation benchmarks ✨

💖 The benchmark problems 💖

🌸 Input files for benchmarks 🌸

🦋 Run a mocassin pure photoionisation benchmark 🦋

👑 Checking that everything is going OK 👑

💜 Completing the benchmark runs 💜

💀 Comparing with the benchmark tables 💀

🌈 Pure Dust Benchmarks 🌈

🧚 1D Benchmarks 🧚

🍓 2D Benchmarks 🍓

🌟 Running things other than benchmark cases 🌟

🌸 List of keywords 🌸

autoPackets real1 real2 integer3

continuumCube real1 real2

contShape 'string' [optional: real]

convLimit real

debug

densityFile 'string'

densityLaw real...

diffuseSource real1 real2 'string1' integer1 integer2

dustFile string1 string2

dustMass real1

echo real1 real2

edges real1 real2 real3

fillingFactor real1

getEquivalentTau

getTau

📁 Input and Output Files

Input files

Output files

the output files, and bestie there are so many 💖

ionratio.out 👑

lineFlux.out 🌸

temperature.out 💜

equivalentTau.out 🦋

tauNu.out ⭐

tau.out 💀

ionDen.out 🍓

SED.out 🌈

The grid files and photoSource.out 🎀

plot.out 🎀

🌈 Miscellaneous notes on mocassin 🌈

Analytical and Monte Carlo line fluxes 💫

🦋 Running multiple spatial grids 🦋

🧚 Running models that include dust and gas. 🧚

💜 The accessories/ subdirectory 💜

🌈 Plotting dust temperatures 🌈

🌡️ Grain temperature spiking routines 💥

🚧 Limitations and future development 🛠️

🐛 BUGS 🐛

👑 AUTHOR 👑

`autoPackets` real1 real2 integer3

`continuumCube` real1 real2

`contShape` `'string'` [optional: real]

`convLimit` real

`debug`

`densityFile` `'string'`

`densityLaw` real...

`diffuseSource` real1 real2 `'string1'` integer1 integer2

`dustFile` string1 string2

`dustMass` real1

`echo` real1 real2

`edges` real1 real2 real3

`fillingFactor` real1

`getEquivalentTau`

`getTau`

`ionratio.out` 👑

`lineFlux.out` 🌸

`temperature.out` 💜

`equivalentTau.out` 🦋

`tauNu.out` ⭐

`tau.out` 💀

`ionDen.out` 🍓

`SED.out` 🌈

The grid files and `photoSource.out` 🎀

`plot.out` 🎀