Monte Carlo project


\(\)Implement a Monte Carlo model to simulate shortwave radiative transfer in a homogeneous, plane-parallel cloud layer. Compare the results with those from a simple analytic solution. Perform numerical experiments designed to lend insight into radiative processes involving scattering.


There are many very different ways to numerically solve the radiative transfer equation (RTE) for monochromatic radiation in a scattering medium. One of the simplest and most flexible is the so-called Monte Carlo model, which is nothing more than a brute-force simulation of the random trajectories of individual photons. Basically, a photon is released from a specified point in a specified direction, and a random number generator is utilized to determine how far the photon will travel before encountering an extinction event. A second random number determines whether the photon is absorbed or scattered. If it is scattered, two more random numbers determine the direction into which it is scattered. It then proceeds to its next extinction event (again determined by a random number), and the process is repeated until the photon is either absorbed or else it emerges from the top or bottom of the cloud layer to be counted by a “virtual detector”.
The Monte Carlo method is fairly simple to understand, simple to implement, and it can be used for any imaginable cloud geometry. Properties of the radiation field can be obtained at any desired point in the model domain.
The principal drawback to random (or stochastic) methods like Monte Carlo is that one must accumulate statistics for a large number of photons in order to derive intensities that are of high precision. This used to mean enormous computing times for even moderate precision. Nowadays, computers are fast enough that Monte Carlo methods can be used for most problems (if desired) with good success, even if it is not necessarily the most efficient method for some geometries.

Project Description


Your assignment is to implement a Monte Carlo radiative transfer code in order to compute the reflected and transmitted intensities and fluxes associated with an externally illuminated cloud layer. To keep the problem simple, we will assume a plane-parallel cloud layer with specified total optical depth \(\tau\). The single-scatter albedo \(\tilde \omega\)  will be assumed constant throughout the cloud, and the scattering phase function \(p(\Theta)\) will be modeled via the widely used Henyey-Greenstein analytic phase function (see below). Thermal emission will be ignored. The incident illumination will be from a single direction specified as a solar zenith angle \(\theta_0\)  and azimuth \(\phi_0 = 0\). Photons emerging from the cloud layer top and bottom will be tabulated so that reflectivity and transmittance can be computed.


  • Writeup of methods and results (including applicable tables and plots) using the AMS preprint template for LaTeX.
  • Printout of complete Python source code.

Computational details


As before, you are asked to use Python to write your Monte Carlo code. Although Python doesn’t require this, I strongly recommend that you use this project as an opportunity to improve your comfort level with object-oriented programming (OOP).


You will most likely need to import a number of python modules.  Here are the ones I used:

import matplotlib.pyplot as plt
from numpy import log, cos, sin, exp, sqrt
from numpy.random import random as rand
from copy import copy

In my own program, most of the computational details are built into a single object class I called Photon.   This class has the following attributes and methods:

  • tau:  Present vertical position of the photon in the cloud layer, as measured from cloud top (tau = 0.)
  • k:   A 3-element tuple (k[0], k[1], k[2]) representing the unit vector describing the photon’s current direction of travel.
  • weight: Initialized to 1.0 for a new photon; reduced by the single-scatter albedo with each subsequent interaction.  Photon’s life ends when its weight falls below some specified threshold weightmin.
  • scattered: A boolean attribute, initialized to False, indicating whether the photon has ever been scattered.  Photons that exit the cloud with this attribute still set to False are considered to be part of the direct transmittance; otherwise they contribute to the diffuse transmittance.
  • __init__(self) : Initializes the various attributes of a newly instantiated photon.
  • move(): Uses a random number generator to determine the optical distance that a photon travels (see below) before its next extinction and updates its (vertical) position in the cloud layer.  Tests whether the photon has exited the cloud and, if so, calls either self.exittop() or self.exitbase(), as appropriate.
  • exittop(), exitbase(): Tabulates the weights of photons exiting either the top or base of the cloud and subsequently sets the weight to zero, so that the life of that photon effectively ends.
  • scatter(): Sets the photon’s attribute ‘scattered’ to True to indicate that the photon was scattered at least once before exiting the cloud.  Reduces the weight of the photon according to the single scatter albedo.  Uses random numbers to determine the new direction of the photon following scattering.  See below for the mathematics of computing the new direction unit vector k.

To see examples of how a class and its attributes and methods are implemented in Python, see the sample code I provided for the weighting function project and look at the Level and Layer classes.

In addition to defining the Photon class, I define the following global constants with representative values (which may be readily changed for different experiments):

SSA = 0.99
G = 0.85
PI = 3.1416
DEG = PI/180.0
TAUSTAR = 10.0
PHI0 = 0*DEG

Since most of the computational effort is embedded in the Photon class, the top-level program itself is only a few lines long and does the following:

  1. Initializes counters for reflected, diffusely, and directly transmitted photons to zero.
  2. Enters loop  that
    1. Instantiates a new photon
    2. Breaks out of the loop if the desired total number of photons has been exceeded
    3. Enters an inner loop that calls the photon.move() method repeatedly until the photon weight falls below WEIGHTMIN
  3. After exit from the outer loop, computes \(t_\mathrm{dir}\),  \(t_\mathrm{dif}\),  \(r\), and \( a\).
  4. Compares the above results with those computed from the two-stream method.

For running specific experiments, you may wish to enclose the entire sequence described above inside another outer loop that varies one model parameter and collects the results as lists that can be subsequently plotted in the same program rather than writing a separate plotting program.

Generating random variables obeying a prescribed PDF

A computer random number generator, such as the random() function in Python’s Numpy module, typically yields a uniformly distributed pseudorandom real number on the interval \([0, 1]\).  How do we use such a function to generate values following a prescribed PDF?

For the simple case that the desired random variable \(r\) is uniformly distributed on a different interval than [0,1], simple scaling and/or offsets can be used.  For example, in the scattering process, we need an azimuthal angle \(\Phi\) that is uniformly distributed on interval \([0, 2\pi)\), so we use \[ \Phi = 2\pi r .\]

For nonuniform PDFs, a little more math is involved.  Specifically, we use the fact that an arbitrary PDF  \(p(x)\) can be achieved by solving the following for \(x\) in terms of the random number \(r\):\[ r = \int_0^x p(x’) \; dx’ \]

The optical distance traveled

We know that the optical distance a photon travels without extinction is governed by Beer’s Law.  That is, the fraction transmitted beyond a specified optical path \(\tau\) is given by \( t = e^{-\tau}\).  It is less obvious (but easily proven) that this implies a PDF for optical distance traveled
\[ p(\tau) = e^{-\tau} . \] Using the method described for nonuniform PDFs, you will need to obtain an expression for the random optical distance traveled \(\tau_\mathrm{photon}\) in terms of the uniform random variable \(r\).

This will only give you the optical distance traveled along the path of the photon.  This will always be positive.   You will need to convert this to a vertical displacement (positive or negative) within the cloud layer.  Hint:  The third element 0f the unit vector \(\hat{\mathbf{k}}\) of the photon’s direction of travel will help you with this.

The Henyey-Greenstein phase function

The model phase function we will use is the Henyey-Greenstein phase function given by \[ p_\mathrm{HG}(\cos\Theta) = \frac{1-g^2}{(1+g^2 -2g\cos\Theta)^{3/2}} .\]

Again, by calling a random number \(r\) generating yielding a value on the interval [0, 1), you need to generate random values of \(\cos\Theta\) that obey the above phase function interpreted as a PDF.

Getting the photon’s new scattering direction

Before a photon is scattered, it has a direction of travel indicated by the unit vector \(\hat{\mathbf{k}}\).  This direction is expressed in what I call “model coordinates.”  That is,  the elements of  \(\hat{\mathbf{k}}\) correspond to the coordinate directions \((x,y,z)\).  When you use the random number generator to get scattering angles \((\Theta, \Phi)\), these are measured relative to the direction of travel of the photon before it scatters. It therefore takes some matrix algebra to transform the old \(\hat{\mathbf{k}}\) before scattering into the new \(\hat{\mathbf{k}}\) after scattering, given \((\Theta, \Phi)\).

This document gives the derivation of that transformation.   While you should read the first page-and-a-half to understand the reasoning, you only actually require the last half of the last page to implement the transformation in your Python code.

Putting it all together

The above information should be sufficient for you to write a working Monte Carlo code in Python that computes the direct and diffuse transmittance, reflectivity, and absorptivity for a plane-parallel cloud with optical thickness \(\tau^\star\), single scatter albedo \(\tilde \omega\), and asymmetry parameter \(g\).

In addition, add a function that uses the two-stream method (equations 13.25, 13.45, and 13.63-13.66) to compute the same quantities for comparison purposes.

For reference, my own complete Python program, not counting plotting functions, is only about 120 lines long, so it is not a complex program to write if you’re clever about it.

To test your model, run it for a relatively small number of photons at first (e.g., \(N=1000\)) with the following parameters:

SSA = 0.9
G = 0.85

Verify that you’re able to approximately reproduce the values predicted by Beer’s Law for \(t_\mathrm{dir}\) and by the two-stream method for all other variables.  The agreement should be good for the above parameter values, especially once you increase the number of photons to something much larger.

For reference, my program requires about 5 minutes per million photons using the above parameter values.  Time will increase as absorption is reduced and/or optical thickness is increased.

If something doesn’t seem to be working, try choosing values of the above input parameters that will help you isolate the part of the program that is giving erroneous results.


For all applicable experiments, plot the two-stream results together with your Monte Carlo results, using a dashed curve for the former and a solid curve for the latter.  In many cases, you should  plot the reflectivity, total transmittance, and absorptivity on a single plot as functions of whatever parameter is being varied.

Assume the following default values of parameter except where specified in experiments listed below:

SSA = 1.0
G = 0.85
  1. For an optical thickness of 1.0, determine the number of photons required to consistently obtain a direct transmittance that is correct to within 1% relative to the true value according to Beer’s Law.
  2. Determine the effects of varying the optical depth from 0.1 to 100 by factors of 2.
  3. Determine the effects of varying the solar zenith angle from 0 to 85 degrees.
  4. Determine the effects of setting the single scatter albedo to 1.0, 0.9999, 0.999, 0.99, and 0.9.
  5. Determine the effects of varying the asymmetry parameter g from 0.0 to 1.0.
  6. Using parameter values that you choose, test the validity of the similarity transformations described on pp. 416-417.  In particular, are you able to identify combinations of parameters for which the similarity transformations do not appear to yield radiatively equivalent clouds?

Here’s what I mean for the last task:

  1. For a given set of cloud parameters (tau-star, g, SSA), find the similarity-transformed values (tau-star’, SSA’, g=0).
  2. For both sets of values, run your MC model and determine reflectance and TOTAL transmittance (tdir + tdiff).  If the similarity transforms “work”, you should get similar results (r, t, and a) for both sets of values.
  3. Try out different combinations of parameters and see whether you can find any that lead to significant differences in the direct results vs. the similarity-transformed results.

Note that I’m deliberately leaving it up to you to define “significant” and to determine what ranges of parameter values you try out. This is effectively an open-ended science question, and I’m asking you to exercise some creativity and judgment.Note also that you don’t need to provide lots of figures or tabular results — it’s sufficient to provide one or two results that tell the essential story.


In your writeup, please address the following points:

  • Your findings with respect to each of the specific experiments.
  • The degree of consistency or inconsistency between the MC results and the two-stream results under varying conditions.