Introduction

This page is the first of several on Monte Carlo techniques for solving the Radiative Transfer Equation. As an RTE solution technique, these pages could be considered as Level 2 material under the Radiative Transfer Theory chapter. However, because of the importance of Monte Carlo techniques in a wide range of scientific areas, they warrant a chapter unto themselves.

As used to solve the RTE, Monte Carlo techniques refer to algorithms that use probability theory and random numbers to simulate the fates of numerous photons propagating through a medium. Various averages over ensembles of large numbers of simulated photon trajectories give statistical estimates of radiances, irradiances, and other quantities of interest.

Monte Carlo techniques were developed in the 1940s for studies of neutron transport as needed for the design of nuclear weapons (Metropolis, 1949), (Eckhardt, 1987). The name Monte Carlo was at first a code name for this classified research. The name was well chosen because probability and statistics lie at the heart of both the simulation techniques and the gambling games played at the legendary Monte Carlo Casino in Monaco. Monte Carlo techniques are now highly developed and are used to solve many types of problems in physical and biological sciences; finance, economics, and business; engineering; computer graphics for movie production, and pure mathematics.

An essential feature of Monte Carlo simulation is that the known probability of occurrence of each separate event in a sequence of events is used to estimate the probability of the occurrence of the entire sequence. In the ray-tracing setting, the known probabilities that a light ray will travel a certain distance, be scattered through a certain angle, reflect off on a surface in a certain direction, etc., are used to estimate the probability that a ray emitted from a source at one location will travel through the medium and eventually be recorded by a detector at a different location.

The strengths of Monte Carlo techniques are that

• They are conceptually simple. The methods are based on a straightforward mimicry of nature, which in itself endows them with a certain elegance.
• They are very general. Monte Carlo simulations can be used to solve problems for any geometry (e.g., 3D volumes with imbedded objects), incident lighting, scattering phase functions, etc. It is relatively easy to include polarization and time dependence.
• They are instructive. The solution algorithms highlight the fundamental processes of absorption and scattering, and they make clear the connections between the ray-level and the energy-level formulations of radiative transfer theory.
• They are simple to program. The resulting computer code can be very simple, and the tracing of rays is well suited to parallel processing.

The weaknesses of Monte Carlo techniques are that

• They provide no insight into the underlying mathematical structure of radiative transfer theory. The simulations simply accumulate the results of tracing large numbers of rays, each of which is independent of the others.
• They can be computationally very inefficient. Monte Carlo simulation is a “brute force” technique. If care is not taken, much of the computational time can be expended tracing rays that never contribute to the solution, e.g., because they never intercept a simulated detector.
• They are not well suited for some types of problems. For example, computations of radiance at large optical depths can require unacceptably large amounts of computer time because the number of solar rays penetrating the ocean decreases exponentially with the optical depth. Likewise, the simulation of a point source and a point (or very small) detector is difficult. Monte Carlo techniques are based on tracking individual rays in the geometric optics limit of physical optics and thus cannot address wave phenomena such as diffraction.

This page discusses probability distributions and how they are sampled in Monte Carlo simulations. The next page discusses how light rays are tracked in Monte Carlo simulations. Subsequent pages discuss computational tricks for speeding up calculations and improving the accuracy of the statistical estimates, and the design of Monte Carlo simulations for particular problems. Those pages are placed in Level 2 because of their mathematical nature.

Probability Functions

As a ray travels through a medium, the distance it goes between interactions with the medium, whether it will be absorbed or scattered in a given interaction, the new direction it will travel after a scattering event, etc. are all random variables. In mathematics it is customary to let a capital letter, e.g. $X$, denote a random variable (such as the distance a ray travels) and to let a lower case letter, $x$ in this case, denote a value of $X$, e.g. $x=2.2m$ if $X$ is distance. Let $X$ be any such random variable that is defined over a range of values $x_1$ to $x_2$. If $X$ is distance traveled, $x_1=0$ and $x_2=\infty$; if $X$ is a polar scattering angle, $x_1=0$ and $x_2=\pi$ or 180 deg, and so on.

The probability density function (pdf) for $X$, denoted $p_X\left(x\right)$, is a non-negative function such that $p_X\left(x\right)dx$ is the probability (a number between 0 and 1) that $X$ has a value between $x$ and $x+dx$. A pdf must satisfy the normalization

$${\int \; <!--nolimits-->}_{x_1}^{x_2}{p}_X\left(x\right)dx=1.$$

(1)

That is, the probability is 1 that $X$ will have some value in its allowed domain. The cumulative distribution function (cdf) is a function $P_X\left(x\right)$ giving the probability that the random variable $X$ will have a numerical value less than or equal to $x$. The cdf is obtained from the corresponding pdf via

$$P_X\left(x\right)={\int \; <!--nolimits-->}_{x_1}^x{p}_X\left(x^{\prime }\right)d{x}^{\prime }.$$

Note that a pdf has units of 1/{units of x} and can have a magnitude greater than 1 for some values of $x$, whereas a cdf is nondimensional and grows monotonically from 0 at $x=x_1$ to 1 at $x=x_2$.

The mean or average value of $X$ is given by

$${\mu }_X={\int \; <!--nolimits-->}_{x_1}^{x_2}x{p}_X\left(x\right)dx,$$

and the variance of $X$ is given by

$${\sigma }_X^2={\int \; <!--nolimits-->}_{x_1}^{x_2}{\left(x-{\mu }_X\right)}^2{p}_X\left(x\right)dx={\int \; <!--nolimits-->}_{x_1}^{x_2}{x}^2{p}_X\left(x\right)dx-{\mu }_X^2.$$

The pdf for a random variable that can have values only between 0 and 1 is fundamental to Monte Carlo simulation. Let $\Re$ be a random variable drawn from the unit interval between 0 and 1 such that $\Re$ is equally likely to have any value $0\le 𝔯\le 1$ on the interval from 0 to 1. The pdf for $\Re$ is

$$p_{\Re }\left(𝔯\right)=\left\{\begin{array}{c}1if0\le 𝔯\le 1\hfill \\ 0if𝔯<0or𝔯>1.\hfill \end{array}\right.$$

(2)

$\Re$ is said to have a uniform 0 to 1 distribution, denoted by $\Re \sim U\left[0,1\right]$.

We wish to use a randomly drawn value of $\Re$, which is a known number $𝔯$, to determine a value for a random variable $X$. This is done by regarding going from $\Re$ to $X$ as a change of variables. Then the probablity that $\Re$ lies between $𝔯$ and $𝔯+d𝔯$ is the same as the probability that $X$ lies between $x$ and $x+dx$. Thus

$${\int \; <!--nolimits-->}_0^{𝔯}{p}_{\Re }\left({𝔯}^{\prime }\right)d{𝔯}^{\prime }={\int \; <!--nolimits-->}_{x_1}^x{p}_X\left(x^{\prime }\right)d{x}^{\prime }$$

Because $p_{\Re }\left(𝔯\right)$ is known from Eq. (2), the left-hand integral can be evaluated to obtain

$$𝔯={\int \; <!--nolimits-->}_{x_1}^x{p}_X\left(x^{\prime }\right)d{x}^{\prime }=P_X\left(x\right)$$

The fundamental principle of Monte Carlo simulation states that the equation $𝔯=P_X\left(x\right)$ uniquely determines $x$ in such a manner that $x$ lies in the interval $x$ to $x+dx$ with probability $p_X\left(x\right)dx$. That is, drawing a value $𝔯$ from a $U\left[0,1\right]$ distribution and then solving $𝔯=P_X\left(x\right)$ for $x$ gives a randomly determined value of $x$ that obeys the pdf for $X$.

The next sections illustrate how this principle is applied for specific examples of determining ray path lengths and scattering angles.

Determining Ray Path Lengths

Recall from the derivation of the RTE that radiance in a particular direction $\left(𝜃,\varphi \right)$ decays due to absorption and scattering out of the beam according to

$$\frac{dL\left(r,𝜃,\varphi \right)}{dr}=-c\left(r\right)L\left(r,𝜃,\varphi \right),$$

which integrates to give

$$L\left(r,𝜃,\varphi \right)=L\left(0,𝜃,\varphi \right)e^{-{\int \; <!--nolimits-->}_0^{r}c\left(r^{\prime }\right)d{r}^{\prime }},$$

where $c\left(r\right)$ is the beam attenuation coefficient and $r$ is the distance from some starting point. In terms of the optical path length $\tau ={\int \; <!--nolimits-->}_0^{r}c\left(r^{\prime }\right)d{r}^{\prime }$ this is

$$L\left(\tau ,𝜃,\varphi \right)=L\left(0,𝜃,\varphi \right)e^{-\tau }.$$

This experimentally established exponential decay of radiance can be explained in terms of the fate of individual rays if the probability of any particular ray being absorbed or scattered out of the incident direction between $\tau$ and $\tau +d\tau$ is

$$p_T\left(\tau \right)d\tau =e^{-\tau }d\tau$$

Note that this $p_T\left(\tau \right)$ satisfies the normalization condition (1) with $x_1=0$ and $x_2=\infty$. The corresponding cdf is $P_T\left(\tau \right)=1-e^{-\tau }$. Drawing a random number $𝔯$ from a $U\left[0,1\right]$ distribution and solving

$$𝔯=P_T\left(\tau \right)=1-e^{-\tau }$$

for $\tau$ gives

$$\tau =-ln\left(1-𝔯\right).$$

Because $1-\Re$ is also uniformly distributed on $\left[0,1\right]$, we can just as well use

$$\tau =-ln\left(𝔯\right)$$

to determine $\tau$. Optical distances randomly chosen in this manner, when applied to many rays, are consistent with the exponential decay of radiance with distance traveled. If the water is homogeneous, so that $c\left(r\right)$ does not depend on $r$, then $\tau =cr$ and the geometric distance a photon travels is given by

$$r=-\frac{1}{c}ln\left(𝔯\right).$$

Note that the average distance a ray travels is given by

$${\mu }_T={\int \; <!--nolimits-->}_0^{\infty }\tau e^{-\tau }d\tau =1$$

or, for homogeneous water,

$${\mu }_R=\frac{1}{c}.$$

The average distance a ray travels between an absorption or scattering interaction with the water is called the mean free path between interactions. Likewise, the standard deviation ${\sigma }_T$, which is the square root of the variance ${\sigma }_T^2$, of the optical distance traveled is also 1, or ${\sigma }_R=1∕c$. Thus rays travel on average one optical path length, or $1∕c$ meters, but with a large spread of values about that distance.

Determining Scattering Angles

Scattering is an inherently three-dimensional process and must be specified by both polar $\left(\psi \right)$ and azimuthal $\left(\chi \right)$ scattering angles. The scattering phase function $\tilde{\beta }\left({\psi }^{\prime },{\chi }^{\prime }\to \psi ,\chi \right)$ can be interpreted as a pdf for scattering from an incident direction $\left({\psi }^{\prime },{\chi }^{\prime }\right)$ to a final direction $\left(\psi ,\chi \right)$, per unit of solid angle. If we pick a spherical coordinate system centered on the incident direction $\left({\psi }^{\prime },{\chi }^{\prime }\right)$ and recall the expression for an element of solid angle in spherical coordinates, then we can write

$$\tilde{\beta }\left({\psi }^{\prime },{\chi }^{\prime }\to \psi ,\chi \right)d\Omega \left(\psi ,\chi \right)=\tilde{\beta }\left(\psi ,\chi \right)sin\psi d\psi d\chi .$$

Ocean water is usually well described as isotropic medium, which means that there are no optically preferred directions. (This is not the case, for example, in a cirrus cloud with non-randomly oriented ice crystals.) In that case, the polar and azimuthal scattering angles are independent, and we can write

$$\tilde{\beta }\left(\psi ,\chi \right)sin\psi d\psi d\chi =p_{\Psi }\left(\psi \right)d\psi p_X\left(\chi \right)d\chi$$

For an unpolarized beam, the azimuthal angle is equally likely to have any value from 0 to 360 deg, or 0 to $2\pi$ radians. Thus the pdf for azimuthal scattering is $p_X\left(\chi \right)=1∕\left(2\pi \right)$, the cdf is $P_X\left(\chi \right)=\chi ∕\left(2\pi \right)$, and $\chi$ is determined by $\chi =2\pi 𝔯$.

Using this $p_X\left(\chi \right)$ in the previous equation allows us to identify the pdf for the polar angle as

$$p_{\Psi }\left(\psi \right)=2\pi \tilde{\beta }\left(\psi \right)sin\psi$$

Recall from the discussion of the Volume Scattering Function that phase functions satisfy the normalization

$$2\pi {\int \; <!--nolimits-->}_0^{\pi }\tilde{\beta }\left(\psi \right)sin\psi d\psi =1,$$

so this function $p_{\Psi }\left(\psi \right)$ is indeed a pdf. Therefore, to determine the polar scattering angle, we draw a $U\left[0,1\right]$ random number $𝔯$ as always and solve

$$𝔯=2\pi {\int \; <!--nolimits-->}_0^{\psi }\tilde{\beta }\left({\psi }^{\prime }\right)sin{\psi }^{\prime }d{\psi }^{\prime }$$

(3)

for $\psi$.

In general, Eq. (3) must be solved numerically because of the complicated shape of most phase functions, or when the phase function is defined by tabulated data at a finite number of scattering angles and is fit with a spline (or other) function to generate a continuous function of scattering angle. However, a few commonly used phase functions allow this equation to be solved analytically. The simplest case, isotropic scattering, is instructive.

Isotropic Scattering. The phase function for isotropic scattering is $\tilde{\beta }\left(\psi \right)=\frac{1}{4\pi }$. Substituting this into Eq. (3) leads to

$$\psi ={cos}^{-1}\left(1-2𝔯\right)$$

for the determination of the scattering angle. This result may look peculiar until it is remembered that isotropic scattering means that scattering is equally likely into any element of solid angle, not equally likely at every scattering angle $\psi$. Figure 1 illustrates this important point. Scattering from a collimated beam is simulated two different ways. The scattering events occur at the center of a sphere, and the points plotted on the surface of the sphere show the scattering direction. The arrow at the “north pole” represents the direction of the unscattered beam and a scattering angle of 0. The blue lines are lines of constant scattering angle $\psi$, with the thick line at the “equator” being $\psi =90$ deg. In the left panel, the polar scattering angle $\psi$ was drawn from a $U\left[0,180\right]$ distribution, i.e., any value of $\psi$ between 0 and 180 deg was equally likely. The right panel drew $\psi$ from the ${cos}^{-1}\left(1-2𝔯\right)$ distribution. In both simulations the azimuthal scattering angle $\chi$ was drawn from a $U\left[0,360\right]$ distribution. It is visually clear from the left panel that the uniform distribution of $\psi$ values gives too many points near the north pole of the figure, i.e., too many points with scattering angles near 0. The right panel shows an even distribution of points per unit area of the surface of the sphere, i.e., per unit solid angle in any direction.

---

---

Table 1 displays various phase functions and the formulas obtained from solving $𝔯=P_{\Psi }\left(\psi \right)$ used to determine the corresponding scattering angles.

The widely used Fournier-Forand phase function is

$$\begin{array}{llll}\hfill {\tilde{\beta }}^{FF}\left(\psi \right)=& \frac{1}{4\pi {\left(1-\delta \right)}^2{\delta }^{\nu }}\left[\nu \left(1-\delta \right)-\left(1-{\delta }^{\nu }\right)+\left[\delta \left(1-{\delta }^{\nu }\right)-\nu \left(1-\delta \right)\right]{sin}^{-2}\left(\frac{\psi }{2}\right)\right]& \hfill & \\ \hfill +& \frac{1-{\delta }_{180}^{\nu }}{16\pi \left({\delta }_{180}-1\right){\delta }_{180}^{\nu }}\left(3{cos}^2\psi -1\right),& \hfill & \end{array}$$

where

$$\nu =\frac{3-\mu }{2}and\delta =\frac{4}{3{\left(n-1\right)}^2}{sin}^2\left(\frac{\psi }{2}\right),$$

Here $n$ is the real index of refraction of the particles, $\mu$ is the slope parameter of the hyperbolic distribution, and ${\delta }_{180}$ is $\delta$ evaluated at $\psi =180deg$. This phase function has an analytic cdf (Fournier and Jonasz, 1999)

$$\begin{array}{llll}\hfill P_{\Psi }^{FF}\left(\psi \right)=& \frac{1}{\left(1-\delta \right){\delta }^{\nu }}\left[\left(1-{\delta }^{\nu +1}\right)-\left(1-{\delta }^{\nu }\right){sin}^2\left(\psi ∕2\right)\right]& \hfill & \\ \hfill +& \frac{1}{8}\frac{1-{\delta }_{180}^{\nu }}{\left({\delta }_{180}-1\right){\delta }_{180}^{\nu }}cos\psi {sin}^2\psi .& \hfill \text{(4)}\end{array}$$ However, solving $𝔯=P_{\Psi }^{FF}\left(\psi \right)$ for $\psi$ (even if possible) would give a formula so complicated that it is numerically more efficient to use Eq. (3) for $P_{\Psi }^{FF}\left(\psi \right)$ to build up a table of $\psi$ vs. $P_{\Psi }^{FF}\left(\psi \right)$ values for closely spaced values of $\psi$ and the given $n$ and $\mu$ parameters, and then to interpolate within this table to obtain values of $\psi$ as $𝔯$ values are drawn in the course of a simulation. This is illustrated in Fig. 2, which shows a Fourier-Forand cdf for values of $n$ and $\mu$ that give a best fit to the Petzold average particle phase function phase function. The blue arrows show how drawing a value of $𝔯=0.7$ leads to a scattering angle of about 10 deg. When working with tabulated data for highly peaked phase functions, it is usually adequate to use linear interpolation between the tabulated values.

---

---