Modeling Diffuse Reflection (or How to Sample Cosine Distribution)

Molecules impacting surfaces may reflect either specularly or diffusely, depending on surface properties. The diffuse reflection is by far more common and hence it’s important that it’s implemented in the code correctly. In a diffuse interaction, the molecule actually temporarily attaches to the target. The time spent on the surface is known as residence time, and it will be tiny for surfaces at temperatures above the molecule condensation temperature. In this process, the molecule completely “forgets” its original direction and leaves the surface in a new random direction. But there is one important caveat: flux of outgoing molecules is proportional to the cosine from the surface normal. This is known as Lambert’s cosine law. Sampling from this cosine distribution is the topic of this article.

Incorrect Way

Before I show you the correct way to model diffuse reflection, I wanted to talk a bit about the incorrect way. The reason I am mentioning this is that on first sight, this incorrect way actually seems acceptable. In fact, it is what was being used in my plume simulation code Draco for a while. Oops! But first, a bit of geometry.

Computationally, we define a planar surface element (be it a triangle, a quad, or so on) with three local unit vectors: a surface normal, and two tangents. The actual orientation of the tangents is not important since we are sampling random orientation. The important characteristic is that \(\mathbf{n}=\mathbf{t}_{1}\times\mathbf{t}_{2}\). The obvious algorithm for sampling diffuse reflections then may seem to be to select a random multiplier for each of the three unit vectors, add the resulting components together, and convert to a unit vector. Since we are sampling only in the positive half-space of the surface normal (assuming normal pointing outward), the algorithm is as follows (remember, this is incorrect!):

#pick three random numbers
a = -1 + 2*random()
b = -1 + 2*random()
c = random()
 
#multiply by corresponding directions
v1 = [a*n for n in tang1]
v2 = [b*n for n in tang2]
v3 = [c*n for n in norm]
 
#add up to get velocity, vel=v1+v2+v3
vel = []
for i in range(0,3):
    vel.append(v1[i]+v2[i]+v3[i])
 
#make a unit vector
vel_mag = math.sqrt(vel[0]*vel[0]+vel[1]*vel[1]+vel[2]*vel[2])
for i in range(0,3):
   vel[i] /= vel_mag

The full code is listed in random3.py, available at the link below. Now, if we do this for a large number of particles, and for each compute the cosine from the surface normal, we will obtain an angular distribution as shown below. As we can see, the resulting distribution is not right. Majority of particles bounce off at angles closer to 50 degrees, instead of the expected zero.

Sampling from Cosine Distribution

Instead, in order to model diffuse reflections correctly, we need to perform two steps:

1. Sample component along the normal, \(\mathbf{v}_|\), from the cosine distribution
2. Rotate the perpendicular component \(\mathbf{v}_\perp\) about a random \(0:2\pi\) angle

This is given by the following set of equations:
$$\begin{align}
\sin \theta &= \sqrt{\mathcal{R}_1}\\
\cos \theta &= \sqrt{1-\sin^2\theta}\\
\psi &= 2\pi\mathcal{R}_2\\
\mathbf{v}_1 &= \cos\theta \mathbf{n}\\
\mathbf{v}_2 &= \sin\theta \cos\psi\mathbf{t}_1\\
\mathbf{v}_3 &= \sin\theta \sin\psi\mathbf{t}_2\\
\mathbf{v} &= \mathbf{v}_1+\mathbf{v}_2+\mathbf{v}_3
\end{align}
$$
where \(\mathcal{R}\) are two random numbers in the range of [0:1). The computation of cosine comes from the identity that \(\cos^2\theta + \sin^2\theta = 1\).

Numerically, this is written as follows:

sin_theta = math.sqrt(random())   
cos_theta = math.sqrt(1-sin_theta*sin_theta)
 
#random in plane angle
psi = random()*2*math.pi;
 
#three vector components
a = sin_theta*math.cos(psi)
b = sin_theta*math.sin(psi)
c = cos_theta
 
#multiply by corresponding directions
v1 = [a*n for n in tang1]
v2 = [b*n for n in tang2]
v3 = [c*n for n in norm]
 
#add up to get velocity, vel=v1+v2+v3
vel = []
for i in range(0,3):
    vel.append(v1[i]+v2[i]+v3[i])

If we again plot a histogram of the angular distributions, we will obtain a trace as shown below in Figure 2. A much better agreement is seen!

Why sample from sine?

If you look at the code closely, you will notice that we are actually picking a random number that corresponds to the sine of theta, not the cosine. Is this a mistake?

No, it’s not. The reason the sampling works this way is that in order to sample from an arbitrary distribution, you actually need to sample from the cumulative distribution function (CDF) of the probability density function. This is known as the inverse transform method, and you will find more info in the references [1], [2], and [3], noted below. But in short, the CDF tells you the probability of finding some value (x) in the range \(-\infty:x\). In other words, it is the sum of the probabilities up to the value point. Numerically it is written as
$$F_X(x) = \int_{-\infty}^x f_x(\xi) d\xi$$
The importance of CDF is that it ranges from 0 at the left boundary to 1 at the right boundary, which is also what we get from the random number generator. In the case of a periodic function like cosine, we have
$$\begin{align}
F_{cos}(\theta) &= \int_{0}^\theta \cos(\xi) d\xi\\
&=\sin(\theta) – \sin(0)\\
&=\sin(\theta)
\end{align}
$$

Source Code

You can download the two example Python scripts below. Also, I recently started adding the blog examples to Github. You may be interested in cloning the following repo: github.com/particleincell/PICCBlog.
random3.py (illustrates the incorrect way)
cosine.py (correct sampling for diffuse reflection)

Refrences

[1] Generating Random Variates
[2] Lecture 3: Inverse Transform Method
[2] Sampling via Inversion of the cdf