This tutorial finally concludes our journey to simulate Rayleigh Scattering for a planet’s atmosphere.

The next (and final) part will show how to change the shader to also include an additional type of scattering, known as Mie 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 refer to the Atmospheric Scattering Cheatsheet for a complete reference of all the equations used.

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

#### Sampling the View Ray

Let’s recall the equation for the atmospheric scattering that we have recently derived:

The amount of light that we receive is equal to the amount of light emitted from the sun, , multiplied by the sum of the individual contributions from each point in the segment .

We could go one and implement that function directly in our shader. However, there are few optimisations that can be done. It was hinted, in a previous tutorial, that the expression could be simplified even further. The first step we can take is to decompose the scattering function into its two basic components:

The **phase function** and the **scattering coefficient at sea level** are constant with respects to the summation, since the angle and the wavelength do not depend on the sampled point. Hence, they can be factored out:

This new expression is mathematically equivalent to the previous one, but is more efficient to calculate since some of the most heavy parts have been taken out of the summation.

We are not ready to start implementing it. There are infinitely many points we should take into account. A reasonable approximation to is to split into several smaller segments of length , and to accumulate the contribution of each individual segment. In doing so, we assume that each segment is small enough to have a constant density. Generally speaking, that is not the case, but if the is small enough, we can still reach a reasonably good approximation.

The number of segments in is referred to as the **view samples**, since all segments lie on the view ray. In the shader, this will be the
_ViewSamples property. By having it as a property, it is accessible from the material inspector. This allows us to reduce the precision of the shader, in favour of its performance.

The following piece of code allows looping through all the segments in the atmosphere.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
// Numerical integration to calculate // the light contribution of each point P in AB float3 totalViewSamples = 0; float time = tA; float ds = (tB-tA) / (float)(_ViewSamples); for (int i = 0; i < _ViewSamples; i ++) { // Point position // (sampling in the middle of the view sample segment) float3 P = O + D * (time + ds * 0.5); // T(CP) * T(PA) * ρ(h) * ds totalViewSamples += viewSampling(P, ds); time += ds; } // I = I_S * β(λ) * γ(θ) * totalViewSamples float3 I = _SunIntensity * _ScatteringCoefficient * phase * totalViewSamples; |

The variable time is used to keep track of how far we are from the origin point , and is incremented by ds after each iteration.

#### Optical Depth PA

Each point along the view ray brings its own contribution to the final colour of the pixel that we are drawing. That contribution is, mathematically speaking, the quantity inside the summation:

Like we did in the previous paragraph, let’s try to simplify it. We can expand the expression above even further, by replacing with its actual definition:

The product of the transmittance over and becomes:

The **combined transmittance** is modelled as an exponential decay with which coefficient is the sum of the **optical depths** over the path travelled by the light ( and ), multiplied by the *scattering coefficient at sea level* ( with ).

The first quantity that we start calculating is the optical depth for the segment , which goes from the point of entry through the atmosphere, up to the point we are currently sampling in the for loop. Let’s recall the definition of optical depth:

If one had to implement this naively, we would create a function called opticalDepth that samples points between and in a loop. This is possible, but very inefficient. As a matter of fact, is the optical depth of the segment that we are already analysing the outermost for loop. We can save a lot of calculations if we calculate the optical depth the current segmenet centred at ( opticalDepthSegment), and keep accumulating it in the for loop ( opticalDepthPA).

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
// Accumulator for the optical depth float opticalDepthPA = 0; // Numerical integration to calculate // the light contribution of each point P in AB float time = tA; float ds = (tB-tA) / (float)(_ViewSamples); for (int i = 0; i < _ViewSamples; i ++) { // Point position // (sampling in the middle of the view sample segment) float3 P = O + D * (time + viewSampleSize*0.5); // Optical depth of current segment // ρ(h) * ds float height = distance(C, P) - _PlanetRadius; float opticalDepthSegment = exp(-height / _ScaleHeight) * ds; // Accumulates the optical depths // D(PA) opticalDepthPA += opticalDepthSegment; ... time += ds; } |

#### Light Sampling

If we look back at the expression for the light contribution of , we see that the only quantity needed is the optical depth of the segment ,

We will move the code that calculates the optical depth of the segment into a function called
lightSampling. The name comes from the **light ray**, which is the segment that starts at and is directed towards the sun. We have called the point where it exits the atmosphere is .

The lightSampling function, however, will not just calculate the optical depth of . So far we have only considered the contribution of the atmosphere, ignoring the role of the actual planet. Our equations do not account for the possibility that the light ray that goes from towards the sun might hit the planet. If this happens, all calculations done so far must be discarded, as no light would actually reach the camera.

In the diagram above, it’s easy to see that the light contribution of should be ignored, since no sun light is reaching . While looping through the points in between and , the lightSampling function will also check if the planet has been hit. This can be done by checking if the altitude of a point is negative.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 |
bool lightSampling ( float3 P, // Current point within the atmospheric sphere float3 S, // Direction towards the sun out float opticalDepthCA ) { float _; // don't care about this one float C; rayInstersect(P, S, _PlanetCentre, _AtmosphereRadius, _, C); // Samples on the segment PC float time = 0; float ds = distance(P, P + S * C) / (float)(_LightSamples); for (int i = 0; i < _LightSamples; i ++) { float3 Q = P + S * (time + lightSampleSize*0.5); float height = distance(_PlanetCentre, Q) - _PlanetRadius; // Inside the planet if (height < 0) return false; // Optical depth for the light ray opticalDepthCA += exp(-height / _RayScaleHeight) * ds; time += ds; } return true; } |

The function above first calculates the point using rayInstersect. It then divides the segment into _LightSamples segments of length ds. The calculation for the optical depth is the same used in the outermost loop.

The function returns false if the planet has been hit. We can use this to update the missing code the outermost loop, replacing the ....

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
// D(CP) float opticalDepthCP = 0; bool overground = lightSampling(P, S); if (overground) { // Combined transmittance // T(CP) * T(PA) = T(CPA) = exp{ -β(λ) [D(CP) + D(PA)]} float transmittance = exp ( -_ScatteringCoefficient * (opticalDepthCP + opticalDepthPA) ); // Light contribution // T(CPA) * ρ(h) * ds totalViewSamples += transmittance * opticalDepthSegment; } |

Now that we have taken into account all the elements, our shader is complete.

#### Coming Next…

This post (finally!) completes the volumetric shader that simulates atmospheric scattering. So far, we have only taken into account the contribution of the Rayleigh scattering. There are many optical phenomena that cannot be explained with Rayleigh scattering alone. The next post will introduce a second type of scattering, which is known as Mie 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 refer to the Atmospheric Scattering Cheatsheet for a complete reference of all the equations used.

#### Download

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 |

##### 📧 Stay updated

A new tutorial is released every week.

##### 💖 Support this blog

This websites exists thanks to the contribution of patrons on Patreon. If you think these posts have either helped or inspired you, please consider supporting this blog.

How does this shader look like when inside the atmosphere?

Hey!

This particular version doesn’t really work inside the atmosphere.

But it can be changed very easily. The first step is to have a material that does not cull back faces.

Hey Alan, wanted to thank you for the great tutorial. I was able to create my own atmospheric renderer thanks in large part to this series (you can see it in action here: http://davidson16807.github.io/tectonics.js/).

I wanted to point out something I found: you can dramatically improve performance by using a closed form approximation for the “lightSampling()” function in your tutorial. I wrote up a blog post about this, which I hope will make sense: http://davidson16807.github.io/tectonics.js//2019/03/24/fast-atmospheric-scattering.html

Thank you, that is really interesting!

I hope other people will be able to use your implementation!