The direction cosines of any scattered ray are determined using a Monte Carlo technique. Two numbers (RN1 and RN2) are drawn from a random number generator designed to produce values between 0 and 1 with uniform probability. These numbers then define two direction cosines for the scattered ray. The third direction cosine is determined from the constraint that the magnitude of the scattered ray vector is 1.

The question is then how the direction cosines are related to the random numbers. This relationship is determined from integration of the scattering distribution. Suppose that the scattering model is given in terms of a probability, e.g. say that P(q,f) describes the probability of the ray scattering into a range of angles dq (centered about q) and df (centered about f):



where q is the polar angle with respect to the surface normal n and f is the azimuthal angle, as shown above. The necessary integral, performed in polar coordinates, is then:

I = ∫∫ P(q,f) sinq dq df

This indefinite integral (which is still a function of q and f) is then normalized by the Total Integrated Scatter (TIS). The TIS is calculated using the exact same formula given above, but it is a definite integral, i.e. the integration is performed over the full limits of q and f (i.e. the integration goes from 0 to p/2 in q and from 0 to 2p in f). The normalized integral (I_N = I/TIS) is then a function of q and f which – when taking into account the constant of integration that arises from the indefinite integration P(q,f) – is constrained to be between 0 and 1.

For the case when P(q,f) is separable in q and f – i.e. P(q,f) = F(q)*G(f) – I_N can be separated into two terms:

I_N(q,f) = I_q(q)*I_f(f)

Values for the scattered angles are determined by setting each of the two random numbers equal to I_q and I_f, respectively, and then inverting the equations to solve for q and f:

RN1 = I_q(q); q = I_q^-1(RN1)
RN2 = I_f(f); f = I_f^-1(RN2)

If P(q,f) is not separable, then an intermediate step is necessary, but the basic idea is the same. I_q(q,f) is first calculated by integrating over f:

I_q(q,f) = ∫ P(q,f) df

This indefinite integral is normalized to the definite integral (= I_q(q)) calculated by performing the same integration over the full limits of f (0 to 2p). I and TIS are then calculated by performing the indefinite and definite integration of I_q(q), respectively:

I = ∫ I_q(q) sinq dq
I_N = I/TIS

The value for the polar angle is determined by setting the first random number equal to I_N and inverting the equation to solve for q:

RN1 = I_N(q); q = I_N^-1(RN1)

The value for the azimuthal angle is then determined by setting the second random number equal to I_q(q,f) for the known value of q, and inverting to solve for f:

RN2 = I_q(q,f); f = I_q^-1(I_N^-1(RN1),RN2)

For this case when P(q,f) is not separable, the order of integration is in fact not important. We could have chosen to instead integrate P(q,f) over q first, and then over f. The technique for relating the random numbers RN1 and RN2 to q and f would have been identical:

RN1 = I_N(f); f = I_N^-1(RN1)
RN2 = I_f(q,f); q = I_f^-1(I_N^-1(RN1),RN2)

Integration may therefore be performed in whichever order is most convenient for the given probability function.