A simple three-dimensional dust continuum radiative transfer code for different computation methods comparison. Monte Carlo, Quasi Monte-Carlo and Directions Grid Enumeration methods are available. The methods are described in Shulman (2018) and new paper in preparation.
This code is based on mcpolar program by Kenneth Wood (http://www-star.st-and.ac.uk/~kw25/research/montecarlo/montecarlo.html). The code was translated from FORTRAN to C++, refactored and improved. Now it realizes two radiative transfer techniques: Monte Carlo method and directions grid enumeration method (DGEM). DGEM uses precalculated directions of the photons propagation instead of the random ones to speed up the calculations process.
The physics and method detail will be presented in a paper which is in preparation now.
The provided makefile allows compiling the program on Linux.
Cmake configuration can be used to build the program and for IDEs.
Commands for the program compilation for different platforms may be found in .github/workflows configuration files.
PLOT3.plt is a gnuplot script for plotting the resulting images.
In utils folder there are utilities to compare images, add them and obtain a mean image.
Utilities are described in section Utils.
The first program argument is a configuration file.
If the argument is not provided, the program use parameters.json file.
Example parameters.json is provided in the repository root.
More examples are available in regression_tests folder.
The configuration file has the following sections:
Common settings
- fMonteCarlo — is True or False. Set False for using DGEM or set True for using Monte Carlo method
- fWriteScatterings — is True or False. The default value is True. Set False to avoid writing single and double scatterings to separate files.
- taumin — is the minimum optical depth which gives scatterings
- nscat — number of scatterings considered
- SphereSurfaceDirectionsLevel — 2-4 HEALPix grid (Gorski et. al. 2005) for direct photons emission points on the sphere star surface
Monte Carlo parameters
- nphotons — total number of photon packets for Monte Carlo method
- MonteCarloGenerator — Monte Carlo random number generator description
- inverseSphereSourceOrder — For sphere sources: determine the direction of photon propagation before the emission point.
DGEM parameters
- fUseHEALPixGrid — defines the sphere grid type. Set True to use HEALPix grid (Gorski et. al. 2005). Otherwise, a grid based on a partition of a regular icosahedron is used
- PrimaryDirectionsLevel — the primary directions grid parameter. If fUseHEALPixGrid is True the primary directions number is 12 * PrimaryDirectionsLevel ^ 2. Otherwise, it is 5 * 4 ^ PrimaryDirectionsLevel
- SecondaryDirectionsLevel — the secondary directions grid parameter. If fUseHEALPixGrid is True the primary directions number is 12 * SecondaryDirectionsLevel ^ 2. Otherwise, it is primary directions number is 5 * 4 ^ SecondaryDirectionsLevel
- NumOfPrimaryScatterings — number of scatterings in every direction in the primary grid. This parameter should be zero to use optimized one-parameter DGEM. One-parameter DGEM is a recommended mode.
- NumOfSecondaryScatterings — number of scatterings in every direction in secondary grid
- MonteCarloStart — number of scatterings of the photon package after which Monte Carlo method is used instead of DGEM (1 is recommended)
- defaultStarRadius — a default star radius for point sourced used for the inner step computation in one-parameter DGEM method. This parameter is optional. If it is not provided, the solar radius is used.
- dgemBinType — experimental option for one-parameter DGEM implementation. Possible values are Point, Line, and HexLines.
"MonteCarloGenerator" section has the following parameters:
- type — random number generator type
- seed — initialization parameter for pseudo random number generators
- inputRandomFile — a file with the initial state of the random number generator. If it is provided, iseed option is ignored.
- outputRandomFile — a file to store final state of the random number generator
- numberOfPoints — a total number of points for Hammersley (1960) set. The default value is nphotons.
There ane nine supported generator types. Four pseudo-random number generators:
- MinimumStandard — STL implementation following Park (1993).
- MersenneTwister — STL implementation following Matsumoto and Nishimura (2000).
- Ranlux48 — STL implementation following after Lüscher (1994) and James (1994).
- LEcuyer — custom implementation based on L'Ecuyer (1988). The default generator in our program.
Four quasi-random number generators with our implementation:
- Halton — described in Halton (1960) and Halton and Smith (1964)
- Faure — following Faure (1981)
- Sobol — Suggested by Sobol' (1967). The implementation uses Bratley and Fox (1988) and Joe and (2008) results.
- Niederreiter — Niederreiter (1988) base 2 generator implemented after Bratley, Fox, and Niederreiter(1992)
And Hammersley — Hammersley (1960) set.
DGEM supports two types of dust now.
- kappa — The extinction opacity. It is the common dust parameter.
The second Dust description part is "white" or "mie" section.
"white" section corresponds to the dust with Henyey and Greenstein (1941) phase function and White (1979) approximation for polarization functions. It is described by the following parameters:
- albedo — Single scattering albedo
- hgg — Henyey-Greenstein phase function anisotropy
- pl — Peak linear polarization
- pc — Peak circular polarization
- sc — Skew Factor
"mie" section describes the dust with tabulated scattering matrix coefficients. They are usually calculated based on Mie theory. There are two parameters:
- albedo — Single scattering albedo
- tableFile — the path of a file with elements of scattering matrix for Stokes vector (I_r, I_l, U, V). Each line of the file contains a scattering angle in degrees and four coefficients: p1, p2, p3, and p4. Scattering angles from 0 to 180 degrees should be listed in the table in increasing order.
The geometry of the scattering matter is described as a set of geometric shapes. On the top level, there are should be flared disk, sphere envelope, fractal cloud or max/sum list. All angles of this section are measured in degrees, all distances are supposed to be in astronomy units.
The density of the flared disk is described by the following formulae:
where the scale height
and the radial coordinate in the disk plane
.
The model parameters are:
- rInner — the inner disk radius
- rOuter — the outer disk radius
- rho0 — the density at the disk midplane at a radius r0
- h0 — the disk scale height at a radius r0
- r0 — the radius, where h0 and rho0 are defined
- alpha — the radial density exponent
- beta — the flaring power
One can add the disk wind suggested by Safier (1993a, 1993b) to the flared disk. The density of the disk with the wind is the maximum of the wind and disk densities. The density of the Safier wind is
where χ = z / r is the dimensionless height above the disk plane and ρ0 is the wind density on the disk surface at distance 1 AU from the star.
The function η(χ) can be obtained by solving the gas-dynamic equations. ξ'0 and ψ0 are wind model parameters, which are defined in Safier papers (1993a, 1993b) for a list of models.
The parameters of the wind are:
- model — the model of the wind from Safier papers. Should be B, C, D, I, E, F or G
- mOut — the mass outflow rate in solar masses per year
- mStar — the stellar mass in solar masses
- h0 — the dimensionless height from which the wind begins
- rMin — the inner radius of the wind formation region
- rMax — the outer radius of the wind formation region
The conical disk wind (Kurosawa, 2006) may be added to the disk. the wind density is described depending on (ω, l) coordinates. Here ω is the distance from the star to a wind streamline on the equatorial plane, and l is the distance from the equatorial plane to the point along the streamline. The wind density is described by the equation:
The relationship between Cartesian coordinates (x, y, z) and wind parameters (ω, l and others) is expressed by the following relations:
Here D is the distance from the wind 'source' point to a point with coordinates (x, y, z).
δ is an angle between z axis and a streamline.
The local mass outflow rate for ω is obtained based on the dependency
and the total mass outflow rate for the wind.
The wind poloidal velocity is expressed by the equation
Here Rs is the wind scale length, β is the wind acceleration rate, cs is the sound speed in the local disk circular velocity units, vcirc is the disk circular velocity at radius ω, vterminal is the wind terminal velocity.
In sum, the wind parameters are:
- d — the distance from the star to a wind 'source' point
- p — the mass outflow rate on radius dependency power
- mOut — the mass-loss rate in solar masses per year
- mStar — the stellar mass in solar masses
- rInner — the inner radius of the wind formation region in AU
- rOuter — the outer radius of the wind formation region in AU
- rScale — the wind scale length in AU
- windAccelerationRate — the wind acceleration rate
- terminalV — the wind terminal velocity in km/s
- soundSpeed — the sound speed at the wind launching point on the disc in the disk circular velocity units
One can add humps based on the Gaussian function to model disk perturbations. There are two types of humps: a round hump and an azimuthal hump. Both humps may be applied to the flared disk and Safier wind. In the case of the disk hump, it is applied to the disk scale height h. In the case of the wind hump it is applied to the density ρ. Only one hump is allowed for the disk (or wind). The hump centre is located on the x axis of the disk.
The round hump has a shape
The azimuthal hump has a bit more complicated shape
In this equations v is the value the hump changes and other values are hump model parameters:
- h — the hump relative height
- r — the distance from the disk axis to the centre of the hump
- sigma2 — the variance of the hump mass distribution along the disk plane (for the round hump) or the radius (for azimuthal hump)
- sigma2azimuthal — the variance of the hump mass distribution along the azimuth for the azimuthal hump
- sigma2azimuthalBackward — the optional variance of the hump mass distribution along the azimuth for the azimuthal hump. When sigma2azimuthalBackward is provided it is used for negative values of arctan(y, x) and sigma2azimuthal is used only for positive values of arctan(y, x).
The density of the sphere envelope is described by the equation
where the radius is
.
The model parameters are:
- rInner — the inner sphere radius
- rOuter — the outer sphere radius
- rho0 — the density at a radius r0
- r0 — the radius, where rho0 is defined
- alpha — the radial density exponent
The fractal cloud suggested by Elmegreen (1997) is a clumpy dust cloud, which is obtained by the following algorithm:
- Consider a cube space with size 2max, consisting of n3 cubical cells.
- Place dotsN points randomly in the cube.
- For every point build a smaller cub with a center in the point. The size of new cubes is max/Δ, where Δ = dotsN1/dCube. dCube is the fractal dimension.
- Place other dotsN points randomly in every small cube. Do not shift any points outside of the considered big cube and use them in the next steps.
- Repeat steps 3–4 twice more. The total number of created points is dotsN4.
- Shift all points outside the big cube to within it by changing point coordinates per 2max.
- Set the density in every cell proportional to the number of dots in this cell
The model parameters are:
- n — the number of cells along each direction of the cube.
- max — the distance from the cube centre to a cube border along each axis
- dCube — the fractal dimension
- rho0 — the density per one dot
- dotsN — the number of initial dots
- seed — the seed for a random number generator
The list of other geometry shapes (including other lists). The density of the sum list is the sum of densities of all elements. The density of the max list is the maximum density of all elements.
The Flared Disk and the Sphere Envelope may contain additional Translation section. It describes a rotation and a translation of the shape in the global coordinate system. We use the Euler angles for rotation. So there are six parameters: intrinsicRotation, nutation, precession, x, y, z. Intrinsic rotation, nutation and precession are measured in degrees. x, y, z are in AU. All transformations are applied in the global coordinate system, thus the order of operations is:
- Intrinsic Rotation (rotation around z-axis)
- Nutation (rotation around x-axis)
- Precession (rotation around z-axis)
- (x, y, z) vector translation
Every value may be omitted (it will be treated as zero).
Two types of grids are supported: a regular cartesian grid and an unstructured tetrahedral grid.
The cartesian grid has following parameters xmax, ymax, zmax, nx, ny, and nz. The studied area is [-xmax : xmax] x [-ymax : ymax] x [-zmax : zmax] and consists of nx * ny * nz cells. nx, ny, and nz are the numbers of cells along each axis.
The tetrahedral grid is based on the Delaunay triangulation. The area is from -max to max along all axes. The nodes of the grid are described in the nodesFile and elements are listed in elementsFile. Instead of these two files the grid may be presented in the form of one binary file gridBinFile. If all three files exist, grid from nodesFile and elementsFile will be saved in gridBinFile. The grid may be constructed in a separate programs, e.g. gmsh (Geuzaine and Remacle, 2009)
The first line of the nodesFile is the number of nodes. All other lines contain four values: node number and x, y, z coordinates.
The first line of the elementsFile is the number of nodes. All other lines contain five unused values and four node numbers of the tetrahedral vertices.
Is a list of sources with 4 parameters:
- x, y, z — source coordinates
- l — source luminosity
- r — optional star radius in AU. The default value is 0. When r > 0 the star is modelled as a sphere source, otherwise as a point source.
"stars": [
{
"x": 0.0,
"y": 0.0,
"z": 0.0,
"l": 1.0,
"r": 0.0047
}
]We use spherical coordinates [θ, φ] to describe the direction towards each observer. θ is a zenith distance (an angle between positive z-axis direction and observers direction), measured from 0 to 180°. φ is an azimuth, measured counterclockwise in xy plane from the positive x-axis direction. The azimuth may be in the range from -180° to 180° or from 0 to 360°.
There are four additional parameters to describe the observer. We specify them once for all observers.
- rimage is a radius of the visible area in astronomy units. Both axes of the image will be from -rimage to rimage.
- rmask defines the radius in astronomy units of the circular area in the center of the image, which will be covered with a mask (allows you to make an analog of a coronagraph). Optional parameter: 0 is the default value.
- nx and ny are optional parameters describing the image size in pixels. The default size is 200x200.
There are three ways to specify observer directions:
One should specify the list of coordinates φ and θ for each observer position.
All observers are evenly distributed on the circle with a constant θ. In this configuration, one should specify the number of observers numberOfObservers and θ. The first observer has φ equal to zero.
All observers are evenly distributed on the meridian with a constant φ. In this configuration, one should specify the number of observers numberOfObservers and φ. If numberOfObservers is the odd one of the observers has θ equal to 90°. This configuration does not include poles.
You can combine all three ways in one problem. The JSON configuration may be specified in the following way:
"observers": {
"rimage": 800.0,
"rmask": 0.1,
"manual": [
{
"phi": 45.0,
"theta": 45.0
},
{
"phi": 60.0,
"theta": 120.0
}
],
"parallel": {
"numberOfObservers" : 10,
"theta" : 90.0
},
"meridian": {
"numberOfObservers" : 2,
"phi" : 0.0
}
}It may be important to check that the scattering matter has the geometry you expect in some cases. For this purpose, the configuration file contains densitySlice optional section. It allows dumping a slice of the scattering medium along the z-axis into a file. The slice parameters are:
- filename — the output file name
- phi — longitude of the slice in degrees
- radiusMax — the maximum radius value. The radius varies from 0 to radiusMax.
- heightMax — the maximum height value. The height varies from 0 to heightMax.
The effective height is the height (z value) at which the optical thickness of the matter, when moving downward along the z-axis from the infinity, reaches 1. This value is convenient for understanding the geometric shape of the scattering matter. This debug output is described in effectiveHeight optional section. The effective height parameters are:
- filename — the output file name
- radiusMax — the maximum module of x and y coordinates
- heightMax — the maximum height value. Integration along the z axis starts from this height
- dRadius — the step along x and y axes
- dHeight — the step for integration along the z axis
- printCoordinates — boolean flag. If it is true each output line contains three values: x, y and z. Otherwise, only the z values are displayed in the table form.
There are three utilities: imdiff, immean, and imsum.
They can be compiled both with makefile makeutils and Cmake.
Computes image difference statistics and produces three files:
dif.dat with the image pixel difference,
dif_abs.dat with the absolute pixel difference,
and dif_rel.dat with the relative pixel difference.
In the output statistics the main value is refsum.
It is the difference norm for the images.
Usage: ./imdiff [OPTION] FILE1 FILE2
Compare image files FILE1 and FILE2, compute difference norm, and create difference files.
Options:
-h, --help display this help and exit.
-c, --compute just compute difference norm and do not create files.
Computes the sum of images. May be used to combine image files for different scatterings.
Usage: ./imsum [OPTION] FILES -o SUMFILE
Compute sum of images from FILES and store the resulting image into SUMFILE.
Options:
-h, --help display this help and exit.
Computes the mean of images. May be used to combine results of different Monte Carlo simulations with equal photon numbers.
Usage: ./immean [OPTION] FILES -o MEANFILE
Computes mean of images from FILES and stores the resulting image into MEANFILE.
Options:
-h, --help display this help and exit.
- nlohmann/json to parse a configuration file
- catch2 for unit testing
- Bratley P. and Fox B.L., 1992. Algorithm 659: Implementing Sobol's Quasirandom Sequence Generator. ACM Trans. Math. Softw. 14, 88–100. doi:10.1145/42288.214372
- Bratley P., Fox B.L., and Niederreiter H., 1992. Implementation and Tests of Low-Discrepancy Sequences. ACM Trans. Model. Comput. Simul. 2, 195–213. doi:10.1145/146382.146385
- Elmegreen B.G., 1997. Intercloud structure in a turbulent fractal interstellar medium. Astrophys. J. 477, 196–203. doi:10.1086/303705
- Faure H., 1981. Discrépances de suites associées à un système de numération (en dimension un). Bulletin de la S. M. F. 109, 143–182. doi:10.24033/bsmf.1935
- Fox B.L., 1986. Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators. ACM Trans. Math. Softw. 12, 362–376. doi:10.1145/22721.356187
- Geuzaine C. and Remacle J.-F., 2009. Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering, 79(11), 1309-1331. doi:10.1002/nme.2579
- Gorski K.M., Hivon E., Banday A.J., Wandelt B.D., Hansen F.K., Reinecke M., and Bartelmann M., 2005. HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere, Astrophys. J., 622, 759–771. doi:10.1086/427976
- Halton J., 1960. On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals. Numerische Mathematik 2, 84-90. doi:10.1007/BF01386213
- Halton J. and Smith G. B., 1964. Algorithm 247: Radical-inverse quasi-random point sequence. Commun. ACM 7, 701–702. doi:10.1145/355588.365104
- Hammersley J., 1960. Monte Carlo methods for solving multivariable problems. Proceedings of the New York Academy of Science 86, 844-874. doi:10.1111/j.1749-6632.1960.tb42846.x
- Henyey L.G. and Greenstein J.L., 1941. Diffuse radiation in the Galaxy. Astrophys. J. 93, 70–83. doi:10.1086/144246
- James F., 1994. RANLUX: A Fortran implementation of the high-quality pseudorandom number generator of Lüscher. Computer Physics Communications 79, 111-114. doi:10.1016/0010-4655(94)90233-X
- Joe S. and Kuo F.Y., 2008. Constructing Sobol Sequences with Better Two-Dimensional Projections. SIAM Journal on Scientific Computing 30, 2635–2654. doi:10.1137/070709359
- Kurosawa R., Harries T. J., and Symington N. H., 2006. On the formation of Hα line emission around classical T Tauri stars. MNRAS 370, 580–596. doi:10.1111/j.1365-2966.2006.10527.x
- L'Ecuyer P., 1988. Efficient and Portable Combined Random Number Generators. Commun. ACM 31, 742–749,774. doi:10.1145/62959.62969
- Lüscher M., 1994. A portable high-quality random number generator for lattice field theory simulations. Computer Physics Communications 79, 100–110. doi:10.1016/0010-4655(94)90232-1
- Matsumoto M. and Nishimura T., "Dynamic Creation of Pseudorandom Number Generators", Monte Carlo and Quasi-Monte Carlo Methods 1998, Springer, 2000, 56–69. doi:10.1007/978-3-642-59657-5_3
- Niederreiter H., 1988. Low-discrepancy and low-dispersion sequences. Journal of Number Theory 30, 51–70. doi:10.1016/0022-314X(88)90025-X
- Park S.K., Miller K.W., and Stockmeyer P.K., (1993). "Technical Correspondence: Response". Communications of the ACM. 36, 108–110. doi:10.1145/159544.376068
- Safier P. N., 1993a. Centrifugally Driven Winds from Protostellar Disks. I. Wind Model and Thermal Structure. Astrophys. J. 408, 115. doi:10.1086/172574
- Safier P. N., 1993b. Centrifugally Driven Winds from Protostellar Disks. II. Forbidden-Line Emission in T Tauri Stars. Astrophys. J. 408, 148. doi:10.1086/172575
- Shulman S.G., 2018. Three-dimensional heuristic radiation transfer method based on enumeration using the directions grid. Astronomy and Computing 24, 104–116. doi:10.1016/j.ascom.2018.06.002
- Sobol’ I.M., 1967. Distribution of points in a cube and approximate evaluation of integrals. U.S.S.R Comput. Maths. Math. Phys. 7, 86–112. doi:10.1016/0041-5553(67)90144-9
- White R.L., 1979. Polarization in reflection nebulae. I. Scattering properties of interstellar grains. Astrophys. J. 229, 954–961. doi:10.1086/157029