1
E(r)

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

The irradiance E(r) [W/cm²] 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

2
p(r)

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:

3
F(r)

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

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

4
r₁

Equate F(r₁) = rnd and solve for r₁.

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

5
Launch

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.