ECE580CLT: howtosample.html

Return to Monte Carlo Modeling page.

How to sample probability density functions with the Monte Carlo method.
A simple example using a uniform-irradiance flat-top circular beam.

The incident beam is a circular beam of uniform irradiance b with a total power of 1 W. The radius of the beam is a. Here is a step-by-step protocol for determining how to launch photons. This general method could be applied to a Gaussian beam or to some other beam with some specific irradiance pattern.


Specify the irradiance of your incident beam, E(r).

The irradiance E(r) [W/cm2] in this example is a flat-top beam of radius a. In other words, E(r) is uniform at a constant value we shall call b:

E(r) = b for r<=a
E(r) = 0 for r>a


Specify the probability density function, p(r).

The probability density p(r) is defined such that:

Therefore, letting E(r) = b and substituting for p(r) = E(r) 2 pi r in the above equation, we get:

So we now have p(r).


Integrate p(r) to yield probability distribution function, F(r).

The integration of p(r) from r = 0 to a particular r1 yields the probability distribution function, F(r1):


Equate F(r1) = rnd and solve for r1.

The key to the Monte Carlo sampling is to equate the F(r1) based on p(r) with a function F(rnd1) based on a random number rnd generated by a random number generator. But F(rnd1) = rnd1, so one is simply equating F(r1) = rnd1. Then one solves for r1 as a function of rnd1:


Rules for launching photons.

Now you can launch photons by specifying the launch position (x,y,z) and trajectory (ux, uy, uz), where (ux, uy, uz) are the cosines of the trajectory angle projected onto each of the three Cartesian axes. Because our example Monte Carlo program is scoring photon transport in cylindrical coordinates and hence is cylindrically symmetric, launching at a radial position r can be accomplished by setting y = 0 and setting x = r. Because our example is a collimated beam irradiating perpedicular to the tissue surface, the trajectory is pointed directly downward, so uz = 1 and ux and uy are zero. The photon launch is accomplished by the programming statements:

rnd = RandomGen(1,0,NULL); /* creates a random number */
x = a*sqrt(rnd);
y = 0;
z = 0;
ux = 0;
uy = 0;
uz = 1;

If our Monte Carlo was not cylindrically symmetric and we needed to uniformly launch over the x-y plane of our circular spot, then after choosing the r of launch we must choose the portion of r that is assigned to x and the portion assigned to y. The programming statements would be:

rnd = RandomGen(1,0,NULL); /* creates a random number */
r = a*sqrt(rnd);
rnd = RandomGen(1,0,NULL); /* creates a second random number */
x = r*cos(rnd*2*pi);
y = r*sin(rnd*2*pi);
z = 0;
ux = 0;
uy = 0;
uz = 1;

The result is that photons will be launched uniformly over a circular spot of radius a at an irradiance of b.