imcrt

ice Monte Carlo Radiative Transfer.

Dependencies

A working installation of MATLAB. The test suite runs on R2025b with no toolboxes. Developed on R2020b and Octave 7.2.0. Older releases and GNU Octave are nominally supported but untested. The project has no dependencies other than the code included here.

Install

Run Setup.m. If running in Octave, check .octaverc.

Usage

Run mcrt_verify to verify model accuracy. The reflect case compares hemispherical and angular reflectance and transmittance with van de Hulst's tabulated solutions to the transfer equation (Vol. 2, Table 35). The fluence case runs the internal-fluence problem from Wang et al. (1995).

The examples directory includes code needed to reproduce the detector interference simulations reported in Cooper et al 2021, the work this code was developed for. If you wanted to investigate the influence of an instrument on optical measurements, that code would be a good place to start (e.g. see rodintersect.m).

For general use, there is a library of "inherent optical properties" (scattering and absorption coefficients) for water ice Ih saved in dat/mie_iops_dE.mat. The prefix mie_ refers to the Mie scattering formulas used to compute the scattering coefficients. The suffix _dE refers to the "delta Eddington" approximation used to compute the extinction coefficients: $c=\sqrt{3ab_e}$ with extinction coefficient $c$, absorption coefficient $a$, and effective (or 'reduced') scattering coefficient $b_e$. Note that this equation is identical to the diffusion approximation, where $c$ is sometimes called the "propagation coefficient". It's inverse $1/c$ is the "transport length". See doc/tc-2020-53-supplement.pdf for more details on how these values enter into the Monte Carlo model.

dat/ssa_iops_dE.mat contains the same values of absorption coefficient, but values of scattering and extinction coefficient computed with the "specific surface area" approximation, also called the "geometric optics" approximation. This approximation is valid for scatterers about the same size or slightly larger than the interacting wavelengths. The effective particle radii are saved in the libary.

Verification

mcrt_verify runs one case at a time as a set of seeded simulations (8 for reflect, 4 for fluence) and prints a PASS/FAIL table. The reflect case compares hemispherical and angular reflectance and transmittance with van de Hulst's tabulated solutions (Vol. 2, Table 35). The fluence case runs the internal-fluence problem from Wang et al. (1995) with self-consistency checks, since no numeric reference is in the repository. Call mcrt_verify('fluence') to pick a case; the default is reflect, and the function returns the verdict table and the runs. Every run carries its own standard errors for the reflectance and transmittance outputs in RT.se; with fewer than four runs the verdict uses those instead of the spread over runs.

For the full multi-run check, run the driver from the repo root (MATLAB only: it uses datetime and the JVM's SHA-256). Run the dry run at scale = 1e-2 first. The driver checkpoints every run and writes a dated PASS/FAIL report and one figure per case under verify_out/:

Setup; addpath('tests/verify'); vdhverify('verify_out', 1e-2)
Setup; addpath('tests/verify'); vdhverify('verify_out')

The cases are in tests/verify/verifycases.m. To add a case, add an entry there, add a verdict function for its reference, and add a branch to mcrt_verify. The fast test suite is runtests('tests'); tests/README.md describes it and the on-demand perf suite.

Limitations

The model is a plane-parallel, homogeneous slab with a vertical pencil source. Scattering follows the Henyey-Greenstein phase function. Lengths are in cm.

• There is no refractive-index boundary, so the model has no Fresnel reflection, no refraction, and no specular term. Direct reflectance is therefore zero.
• The tally options default to wmin 1e-4, wrr 10, R 2 cm, dr 0.001 cm, and da pi/60, which gives 30 angular bins over the hemisphere. Set any of them using name-value input options in mcrt.
• The roulette threshold comparison (1e-4 against 1e-5) was measured only for the van de Hulst case at albedo 0.9, where no packet reaches either threshold. The fluence case at albedo 0.999 was not compared. Roulette is unbiased, so the threshold changes only the variance and the run time.
• The radial bin beyond R pools every packet at r > R and is normalized as one more ring of width dr. Its per-area density is meaningless when much light leaves past R.
• The depth bin beyond Z stays empty. A packet that crosses Z is scored as transmitted before any absorption is tallied.
• RT.se gives per-run standard errors for the reflectance and transmittance outputs only. Absorption and fluence errors come from the spread over multiple seeded runs.
• The fluence case is verified by reproducing Fig. 4 of Wang et al. 1995. The repository holds no numeric fluence reference. Only the reflect case has a tabulated reference.
• chgdir takes its on-axis branch only below a polar angle of 1e-12 rad. A direction between 1e-12 and 1e-6 rad loses about eps/theta^2 of its unit norm. Such an angle occurs about once in 1e12 scatters for near-isotropic directions and about once in 2e6 scatters at g = 0.999. The verification and paper cases use g up to 0.9.

tests/reports/Post-publication corrections.md states the corrections made to src/mcrt.m after publication and their sizes. tests/reports/impact-report.md quantifies them.

How do I cite this?

If you find this model useful, please consider citing the software release, and/or the following paper.

Cite the version DOI of the release you ran, so the citation resolves to that exact code. Each release lists its version DOI in its GitHub release notes and on its Zenodo record. Note: Zenodo mints that DOI after the tag exists, so a downloaded release archive carries only the concept DOI 10.5281/zenodo.4579072, which resolves to all versions. CITATION.cff on the repository branches carries both. CHANGELOG.md records what changed between releases.

Cooper, M.G., Smith, L.C., Rennermalm, A.K., Tedesco, M., Muthyala, R., Leidman, S.Z., Moustafa, S.E., Fayne, J.V., 2021. Spectral attenuation coefficients from measurements of light transmission in bare ice on the Greenland Ice Sheet. The Cryosphere 15, 1931–1953. https://doi.org/10.5194/tc-15-1931-2021

@article{cooper_2021_TC,
  title = {Spectral Attenuation Coefficients from Measurements of Light Transmission in Bare Ice on the {{Greenland Ice Sheet}}},
  author = {Cooper, M. G. and Smith, Laurence C. and Rennermalm, {\AA}sa K. and Tedesco, Marco and Muthyala, Rohi and Leidman, Sasha Z. and Moustafa, Samiah E. and Fayne, Jessica V.},
  year = {2021},
  month = apr,
  journal = {The Cryosphere},
  volume = {15},
  number = {4},
  pages = {1931--1953},
  publisher = {{Copernicus GmbH}},
  issn = {1994-0416},
  doi = {10.5194/tc-15-1931-2021},
  langid = {english}
}

More details

The model roughly implements the one described in Wang et al. (1995). It is not meant to be an exact implementation. The model simulates transfer through a uniform slab. See doc/tc-2020-53-supplement.pdf and the two references for technical descriptions.

Wang, L., Jacques, S. L. and Zheng, L.: MCML?Monte Carlo modeling of light transport in multi layered tissues, Computer Methods and Programs in Biomedicine, 47(2), 131?146, https://doi.org/10.1016/0169 2607(95)01640 F, 1995.