A Journey Through the Atmosphere

Oct 10, 2017

—

in Maths, Science, Tutorial, Uncategorized

This post describes how to model the density of the atmosphere at different altitude. This is a critical step, since the atmospheric density is one of the parameters necessary to correctly calculate the Rayleigh scattering.

You can find all the post in this series here:

Part 1. Volumetric Atmospheric Scattering
Part 2. The Theory Behind Atmospheric Scattering
Part 3. The Mathematics of Rayleigh Scattering
Part 4. A Journey Through the Atmosphere
Part 5. A Shader for the Atmospheric Sphere
Part 6. Intersecting The Atmosphere
Part 7. Atmospheric Scattering Shader
🔒 Part 8. An Introduction to Mie Theory

You can download the Unity package for this tutorial at the bottom of the page.

Atmospheric Density Ratio

Something that we have not addressed yet is the role of the atmospheric density ratio $\rho$ . From a logical point of view, it makes sense that the strength of the scattering is proportional to the atmospheric density. More molecules per square metre means more chances for photons to be scattered. The challenge is that the composition of the atmosphere is very complex, consisting of several layers with different pressures, densities and temperatures. Luckily enough, most of the Rayleigh scattering happens in the first 60Km of the atmosphere. Within the troposphere, the temperature decreases linearly and the pressure decreases exponentially.

The diagram below shows the relationship between density and altitude in the lower atmosphere.

The value of $\rho \left(h\right)$ represents the atmospheric measure at an altitude of $h$ metres, normalised so that it starts at zero. In many scientific papers, $\rho$ is also called the density ratio because it can be also defined as:

$\rho\left(h\right) = \frac{density\left(h\right)}{density\left(0\right)}$

Dividing the actual density by $density\left(0\right)$ causes $\rho\left(h\right)$ to be $1$ at sea level. As highlighted before, however, computing $density\left(h\right)$ is far from being trivial. We can approximate it with an exponential curve; some of you might in fact have recognised the density in the lower atmosphere follows an exponential decay.

If we want to approximate the density ratio with an exponential curve, we can do it like this:

$\rho\left(h\right) = exp\left\{-\frac{h}{H}\right\}$

where $H_0$ is a scaling factor called scale height. For the Rayleigh scattering in the lower atmosphere of Earth, it is often assumed $H=8500$ metres (diagram below). For the Mie scattering, it is often around $1200$ metres.

The value used for $H$ does not give the best approximation for $\rho\left(h\right)$ . This, however, is not really an issue. Most of the quantities presented in this tutorial are subjected to severe approximations. For best looking results, it will be much more effective to tweak the available parameters to match reference images.

Exponential Decay

In the previous parts of this tutorial, we have derived an equation that shows how to account for the out-scattering that a ray of light is subjected to after interacting with a single particle. The quantity used to model this phenomenon was called the scattering coefficient $\beta$ . We have introduced the coefficients $\beta$ to take that into account.

In the case of the Rayleigh scattering, we have also provided a closed form to calculate the amount of light that is subjected to atmospheric scattering per single interaction:

$\beta \left(\lambda, h \right )=\frac{8\pi^3 \left(n^2-1 \right )^2}{3}\frac{\rho\left(h\right)}{N}\frac{1}{\lambda^4}$

When evaluated at sea level, which means using $h=0$ , it produces the following results:

$\beta\left(680nm\right) = 0.00000519673$

$\beta\left(550nm\right) = 0.0000121427$

$\beta\left(440nm\right) = 0.0000296453$

Where $680$ , $550$ and $440$ are the wavelength which loosely maps to red, green and blue.

What is the meaning of those numbers? They represent the ratio of light that is lost by a single interaction with a particle. If we assume a ray of light has initial intensity $I_0$ and traverses a segment of the atmosphere with (a generic) scattering coefficient $\beta$ , then the amount of light that is not lost to scattering is:

$I_1=\underset{\text{initial energy}}{\underbrace{I_0}} - \underset{\text{energy lost}}{\underbrace{I_0 \beta}}=I_0 \left(1-\beta\right)$

While this holds for a single collision, we are interested in understanding how much energy is scattered over a certain distance. This means that, at each point, the remaining light is subjected to this process.

When light travels through a uniform medium with scattering coefficient $\beta$ , how can we calculate how much of it survives after travelling for a certain distance?

For those of you who have studied Calculus, this should sound familiar. Whenever a multiplicative process like $\left(1-\beta\right)$ is repeated over a continuous segment, the Euler’s number makes its grand appearance. The amount of light that survives scattering after travelling for $x$ metres is:

$I = I_0 exp \left\{-\beta x \right\}$

Once again, we encounter an exponential function. This is not in any way related to the exponential function used to describe the density ratio $\rho$ . The reason why both phenomena are described by exponential functions is that they are both subjected to exponential decay. Besides that, they are no other connections between them.

❓ Where does exp come from?

If you are unfamiliar with Calculus, you might struggle to understand how a process as simple as $I_1 = I_0 \left(1-\beta\right)$ translated into $I = I_0 exp \left{-\beta x \right}$ .

We can see the original expression as a first approximation of our solution. If we want to get a closer approximation, we can take into account two scattering events, halving the scattering coefficient for each one:

$I_1 = I_0 \left(1-\frac{\beta}{2}\right)$

$I_2 = \boxed{I_1} \left(1-\frac{\beta}{2}\right)$

Which leads to:

$I_2 = \boxed{I_0 \left(1-\frac{\beta}{2}\right)} \left(1-\frac{\beta}{2}\right)=$

$=I_0 \left(1-\frac{\beta}{2}\right)^2$

This new expression for $I_2$ indicates how much energy survives two collisions. What about three? Or four? Or five? We can generalise this with the following expression:

$I=\lim_{n\rightarrow \infty }I_0 \left(1-\frac{\beta}{n}\right)^n$

where $\lim_{n\rightarrow \infty }$ is the mathematical construct that allows repeating a process infinitely many times. This is necessary since $\infty$ is technically not a number, hence it does not make sense in this context to write something like \frac{\beta}{\infty}.

That expression is the very definition of the exponential function:

$\lim_{n\rightarrow \infty } \left(1-\frac{\beta}{n}\right)^n = e^{-\beta}=exp\left{-\beta\right}$

which describes the amount of energy that survives a multiplicative process over a continuous interval.

Uniform Transmittance

In the second part of this tutorial we have introduces the concept of transmittance $T$ , as the ratio of light that survives the scattering process after travelling through the atmosphere. We now have all the elements we need to finally derive an equation that describes it.

Let’s look at the diagram below, and see how we can calculate the transmittance factor for the segment $\overline{CP}$ . It’s easy to see that light rays reaching $C$ travel through empty space. Hence, they are not subjected to scattering. The amount of light at $C$ is, therefore, the sun intensity $I_S$ . During its journey to $P$ , some of the light is scattered away from the path; hence, the amount of light reaching $P$ , $I_P$ , will be lower than $I_S$ .

The amount of light scattered depends on the distance travelled. The longer the journey, the strongest the attenuation will be. According to the law of exponential decay, the amount of light at $I_P$ can be calculated as follow:

$I_P = I_S \exp{\left\{-\beta \overline{CP}\right\}}$

where $\overline{CP}$ is the length of the segment from $C$ to $P$ , and $\exp{\left\{x\right\}}$ is the exponential function $e^{x}$ .

Atmospheric Transmittance

We have based our equation on the assumption that the chance of being deflected (the scattering coefficient $\beta$ ) is the same for each point along $\overline{CP}$ . Sadly, that is not the case.

The scattering coefficient strongly depends on the atmospheric density. More air molecules per cubic metre mean higher chances of impact. The density of a planet’s atmosphere is not uniform, but changes depending on the altitude. This also means that we cannot calculate the out-scattering over $\overline{CP}$ in a single step. To overcome this issue, we need to calculate the out-scattering at each point, using its own scattering coefficient.

To understand how this work, let’s start with an approximation. The segment $\overline{CP}$ is split in two, $\overline{CQ}$ and $\overline{QP}$ .

We calculate first the amount of light from $C$ that reaches $Q$ :

$I_Q = I_S \exp{\left\{-\beta{\left(\lambda, h_0\right)} \overline{CQ} \right\}}$

Then, we use the same approach to calculate the amount of light that reaches $P$ from $Q$ :

$I_P = \boxed{I_Q} \exp{\left\{-\beta{\left(\lambda, h_1\right)} \overline{QP} \right\}}$

If we subsitute $I_Q$ in the second equation and simplify, we get:

$I_P = \boxed{I_S \exp{\left\{-\beta{\left(\lambda, h_0\right)} \overline{CQ} \right\}}} \exp{\left\{-\beta{\left(\lambda, h_1\right) \overline{QP} \right\}}=$

$=I_S \exp{\left\{-\beta{\left(\lambda, h_0\right)} \overline{CQ} -\beta{\left(\lambda, h_1\right) \overline{QP} \right\}}$

If both $\overline{CQ}$ and $\overline{QP}$ have same length $ds$ , we can further simplify our expression:

$I_P=I_S \exp{\left\{-\boxed{\left({\beta{\left(\lambda, h_0\right)} +\beta{\left(\lambda, h_1\right)\left) ds} }} \right\}}$

In the case of two segments of equal length with different scattering coefficients, the out-scattering can be calculated by summing up the scattering coefficient of the individual segments, multiplied by the segment lengths.

We can repeat this process with an arbitrary number of segments, getting closer and closer to the actual value. This leads to the following equation:

$I_P = I_S\exp\left\{-\boxed{ \sum_{Q \in \overline{CP}}{\beta\left( \lambda, h_Q \right)}\, ds}\right \}$

where $h_Q$ is the altitude of the point $Q$ .

The approach of splitting a line into multiple segments just like we have done is called numerical integration.

If we assume that the initial amount of light received is equal to $1$ , we obtain the equation for the atmospheric transmittance over an arbitrary segment:

$T\left(\overline{CP}\right) =\exp\left\{-\sum_{Q \in \overline{CP}}{\beta\left( \lambda, h_Q \right)}\, ds\right \}$

We can further expand this expression by replacing the generic $\beta$ with the actual value used by the Rayleigh scattering, $\beta$ :

$T\left(\overline{CP}\right) =\exp\left\{-\sum_{Q \in \overline{CP}}{\boxed{\frac{8\pi^3 \left(n^2-1 \right )^2}{3}\frac{\rho\left(h_Q\right)}{N}\frac{1}{\lambda^4}}}\, ds\right \}$

Many factors of $\beta$ are constant and can be factored out of the summation:

$T\left(\overline{CP}\right) =\exp\left{\underset{\beta\left(\lambda\right)}{\underset{\text{constant}}{\underbrace{\frac{8\pi^3 \left(n^2-1 \right )^2}{3}\frac{1}{N}\frac{1}{\lambda^4}}}}\overset{\text{optical depth}\,D\left(\overline{CP}\right)}{\overbrace{\sum_{Q \in \overline{CP}}{\rho\left(h_Q\right)}\, ds}}\right }$

The quantity expressed by the summation is referred to as optical depth $D\left(\overline{CP}\right)$ , and is what we will actually calculate in the shader. The remaining part is a multiplicative coefficient which can be calculated only once, and corresponds to the scattering coefficient at sea level. In the final shader, we will calculate only the optical depth, and provide the scattering coefficients at sea level $\beta$ as an input.

To sum it up:

$T\left(\overline{CP}\right) =\exp\left\{- \beta\left(\lambda\right)D\left(\overline{CP}\right)\right\}$

If you are interested in this topic, I would also suggest reading Carl Davidson post on atmospheric scattering, where he used an improved version of this iterative approach.

Coming Next…

This post explained how to model Earth’s atmosphere. With the next post, we will start writing the shader code necessary to simulate atmospheric scattering.

You can find all the post in this series here:

Part 1. Volumetric Atmospheric Scattering
Part 2. The Theory Behind Atmospheric Scattering
Part 3. The Mathematics of Rayleigh Scattering
Part 4. A Journey Through the Atmosphere
Part 5. A Shader for the Atmospheric Sphere
Part 6. Intersecting The Atmosphere
Part 7. Atmospheric Scattering Shader
🔒 Part 8. An Introduction to Mie Theory

Download

Become a Patron!
You can download all the assets necessary to reproduce the volumetric atmospheric scattering presented in this tutorial.

Feature	Standard	Premium
Volumetric Shader	✅	✅
Clouds Support	❌	✅
Day & Night Cycle	❌	✅
Earth 8K	✅	✅
Mars 8K	❌	✅
Venus 8K	❌	✅
Neptune 2K	❌	✅
Postprocessing	❌	✅
Download	Standard	Premium

atmospheric scattering Mie Rayleigh volumetric rendering

💖 Support this blog

This website exists thanks to the support of patrons on Patreon. If you think my work has helped or inspired you, please consider joining!

👉

Become a member

👈

Comments

8 responses to “A Journey Through the Atmosphere”

王旭

March 31, 2023

I dont understand the “exp{x}”,what does this meaning?

Reply
cc

March 17, 2023

thanks，but i want to know more detail about “exp come from”,what should i do?

Reply
1. Alan Zucconi
  
  March 20, 2023
  
  If you click on it, it will expand a section explaining how exp is derived!
  
  Reply
WeiChen

October 31, 2022

In the tip Where does exp come from, I dont understand why in the limit expression ，the bottom number is n not 2? The process is I0*(1-β/2)*(1-β/2)*(1-β/2)*(1-β/2)….so the result is I0*(1-β/2)^n，isnt it? I am unfamiliar with Calculus,so can you say more about this ? Or you can give me some link about Calculus that help me to understand this. Thank you very much!

Reply
1. Alan Zucconi
  
  November 1, 2022
  
  I think I understand what you mean!
  
  To get CLOSER to the real value, we first do an approximation using TWO steps.
  So we take the value, we multiply it (1-β/2) and then what remains by (1-β/2).
  Here the /2 is added because we repeat the process twice, so we half the contribution of each multiplication by two to compensate.
  
  Then, we can do the same to get an even more precise value, by repeating the process THREE time. But this time we divide β by 3. So it is (1-β/3)*(1-β/3)*(1-β/3).
  Here we repeat the process three time, so we divide each by 3 to compensate.
  
  Then we repeat 4 times, 5 time, …etc.
  If we use the limit, we can effectively apply (1-β/n) exactly n times, with n that goes to infinity.
  
  Reply
Toni

October 27, 2017

Hi Alan!

I just supported you on Patreon, but I can’t seem to figure out how to access to the tutorial, where should I put the password? Could you please help me on this?

Thanks,
Toni

Reply
1. Alan Zucconi
  
  October 27, 2017
  
  Hi Toni!
  First of all, thank you!
  I have recently switched to the Patreon API, but it seems they are not working.
  Other people have been unable to access the tutorial. 🙁
  I am switching back to the previous system.
  If you are a $5+ patron, you should be able to access the entire series from here: https://www.patreon.com/posts/14814109 !
  Let me know if you’re experiencing any other issue, and thank you for your patience!
  
  Reply
  1. Toni
    
    October 27, 2017
    
    Hi Alan!
    
    Thank you for the quick support, now I can access to the tutorials! 🙂
    
    Toni
    
    Reply