# Modeling Diffuse Reflection (or How to Sample Cosine Distribution) Posted on April 2nd, 2015
Previous Article :: Next Article

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*vel+vel*vel+vel*vel)
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. Figure 1. Randomly sampling three directions will not result in the correct cosine distribution. Black line is the expected distribution, the red line is what we got.

## 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 , , and , 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

Subscribe to the newsletter and follow us on Twitter. Send us an email if you have any questions.

### 9 comments to “Modeling Diffuse Reflection (or How to Sample Cosine Distribution)”

1. Jonas Collberg
March 13, 2016 at 12:02 pm

What’s the reasoning behind the sqrt in:

sin_theta = math.sqrt(random())

2. Jonas Collberg
March 13, 2016 at 12:17 pm

I guess if x is uniform in [0, 1] then also sqrt(x) is uniform in [0, 1], but it only seems to be an unnecessary computation?

3. Jonas Collberg
March 17, 2016 at 3:33 pm

Lubos, thanks for the link. I found this (and reading up on the inverse transform method) to be helpful in understanding the derivations.

Another question though: In your Python code example, you calculate the area of each “bin” using the formula:
 A = 2 * pi * ((1-cos(t2)) - (1-cos(t1))) 

Could you explain this a bit more? Is this taking Lambert’s cosine factor into account, or is it just plain geometry? Any clarification would help, thanks.

4. David Miller
April 14, 2016 at 7:59 pm

Thank you for this post. I was pulling my hair out assuming particles left a surface in completely random directions.

What sort of surfaces reflect specularly? Also, how does that affect the desorption of particles?

• April 15, 2016 at 9:08 am

Not sure – I generally use diffuse reflection for all surface impacts.

5. David Miller
April 16, 2016 at 8:08 pm

Thanks Lubos;

I know most surfaces reflect light specularly. Do you know of any that reflect gas particles specularly? Perhaps it depends on surface preparation?

• April 19, 2016 at 7:50 am

I spoke to some colleagues about this and it seems that in general molecules reflect diffusely. However, there may be some cases of high speed molecules impacting cryogenic surfaces that result in specular-like reflections: https://www.princeton.edu/~gscoles/theses/sean/chapter2.html

6. Pierre
July 29, 2016 at 10:14 am

Hi,

I’m doing a simulation of particles in a box in 2D, and instead of reflection for the collisions on the wall I would like to use the cosine law.

so for example when a particle hit the upper wall for the reflection case :

v[j]=-1*(v[j])#Vy reflected

but now if I want to use cosine law:

sin_theta = math.sqrt(random.random())
cos_theta = math.sqrt(1-sin_theta*sin_theta)

#random in plane angle
psi = random.random()*2*math.pi;

v[j]=math.sqrt(v[j]*v[j]+v[j]*v[j])*sin_theta #Vx
v[j]=math.sqrt(v[j]*v[j]+v[j]*v[j])*cos_theta #Vy

With that I thought I would end up with a distribution equivalent to the cosine distribution with phi=0, i.e a quarter of a circle. But I ends up with a probability of reflection around 50 degrees.

Do you have any idea why ?

Anyway thank you for this very cool website.