Simple Atmosphere
By: Gijs Bellaard
On earth the white light from the sun is scattered/absorbed by the atmosphere, creating the beautiful colors in the sky we are all familiar with. In this article I will derive a simple approximate mathematical model of how the atmosphere looks, which can be used to render the sky for computer graphics purposes.
Let \(u \in \bbR^3\) be a normalized vector indicating which way is up in our world. We will model the atmosphere as a completely homogeneous volume of gas occupying: \[ \Omega = \{ p \in \bbR^3 \mid \Ang{p, u} \lt 0 \} \] where \(\Ang{\cdot,\cdot}\) is the standard dot product on \(\bbR^3\). In other words, \(\Omega\) is the halfspace in the opposite direction from the up direction. Pick any point \(p \in \Omega\) in this halfspace and suppose we're looking from this point \(p\) in the normalized direction \(d \in \bbR^3\). We want to determine how much light \(L(p,d)\) enters our eyes in this configuration.
Let \(s\) be the normalized vector pointing to where the sun is in the sky. We will model the intensity of the sun in the sky \(S(d)\) as a perfect disc having an angular radius of \(\theta\) in radians. \[ S(d) = \begin{cases} 1 & \text{if } \arccos \Ang{d,s} \leq \theta \\ 0 & \text{otherwise} \end{cases} \] where \(\arccos \Ang{d,s} \) calculates the angle between the sun vector \(s\) and the direction vector \(d\). As an indication, the angular radius of the sun is about \( 0.0045 \) radians or equivalently \(0.26\) degrees, pretty small!
In real life every frequency of light is scattered/absorbed differently in the atmosphere: this is precisely how the white light of the sun can turn into different colors in the sky. To simplify the situation we will only consider one frequency of light and look at how this frequency scatters/is absorbed in our atmosphere. Once we know how one frequency behaves we can make an approximation for the complete visible spectrum of light by only considering three frequencies our eyes are most sensitive to: one red, one green, and one blue frequency.
Consider one ray of light from the sun with one specific frequency. Once this ray of light is in the atmosphere it has a chance to bump into a particle. When this happens the ray of light either reflects in a random direction, i.e. scatters, and continues its course, or is absorbed by the particle. Suppose that the probability that the ray of light travels a distance of \(l \gt 0\) without being interrupted by particles is \(0 \lt p \lt 1\). This chance is called the transmittance. We can deduce, from the assumed homogeneity of the atmosphere, that the transmittance for a distance of \(2l\) is \(p^2\). More generally, we have that the transmittance over a distance of \(al\) is \(p^a\) for every \(a \gt 0\). Let \(T(l)\) denote the function that gives the transmittance for every distance \(l\). We can deduce from the previous argument that the transmittance must be of the form: \[T(l) = e^{-\sigma l}\] where \(\sigma\) is called either the scattering or absorption coefficient, depending on what kind of phenomenon is being discussed. Indeed, we have that: \[ T(a\cdot l) = e^{-\sigma a l} = \Par{e^{-\sigma l}}^a = T(l)^a \] Stated in a different way, transmittance has the following property: \[ T(l_1 + l_2) = T(l_1) T(l_2) \] For all \(l_1, l_2 \geq 0\). This property, that is satisfied by all exponential functions, will be of use later.
There are an infinite amount of ways light from the sun can reach our eyes: either directly without being scattered, or indirectly because it scattered once and then made it to our eyes, or it scattered twice, or thrice, etc... Simulating every possible way is impossible, so we will make the approximation of only considering direct light, denoted by \(D(p,d)\), and light that has scattered once, denoted by \(I(p,d)\). We denote this with an \(I\) because this kind of scattering is known as inscattering: it scatters in towards our direction. We thus set the total received light \(L(p,d)\) equal to: \[ L(p,d) = D(p,d) + I(p,d) \]
For the direct light the computation is quite simple. We only need to know through how much atmosphere the light from the sun needs to travel before reaching our eye. Let \(l(p,d)\) denote this length. If we are looking down into the atmosphere, which is when \(\Ang{u,d} \leq 0\), the light will travel through an infinite amount of atmosphere. Otherwise, the light will travel a distance of \(-\Ang{u,p} / \Ang{u,d}\). \[ l(p,d) = \begin{cases} -\Ang{u,p} / \Ang{u,d} & \text{if } \Ang{u,d} \gt 0 \\ \infty & \text{otherwise} \end{cases} \] The intuition here is that \( -\Ang{u,p} \) captures how deep we are within the atmosphere, and that \(\Ang{u,d}\) captures how fast we're moving up out of the atmosphere. One interesting property of \(l(p,d)\) that we will need later is that: \[ l(p_1+p_2,d) = l(p_1,d )+l(p_2,d) \] or more specifically: \[ l(a\cdot p,d) = a\cdot l(p,d) \] for every \(a \in \bbR\). Anyway, we can now calculate that the direct light \(D(p,d)\) is equal to the sun intensity in the direction of \(d\) multiplied by the transmittance determined using \(l(p,d)\): \[ D(p,d) = S(d) T(l(p,d)) \]
For the inscattering the formula is a bit more involved. This stems from the fact there are still an infinite amount of ways light can scatter in our direction. To capture this behavior we need to integrate over all these possibilities. Without going into detail, the following formula holds: \[ I(p, d) = \int_0^{l(p, d)} \sigma S(s) T(l(p+\tau d, s)) T(\tau) d\tau \] Using the properties of \(T\) and \(l\) we can transform this equation into: \[ = S(s) T(l(p, s)) \int_0^{l(p, d)} \sigma T(f \tau) d\tau \] Where we have used the shorthand \(f = l(d,s)+1\). Using the definition of transmittance we get that the integral term is equal to: \[ \int_0^{l(p, d)} \sigma T(\tau f) d\tau = \int_0^{l(p, d)} \sigma e^{-\sigma f \tau} d\tau \] This is an easy integral to solve: \[ = \frac{1}{f}\Par{1-e^{-\sigma f l(p,d)}} = \frac{1}{f}\Par{1-T(fl(p,d))} \] So, all in all, the inscattering term is: \[ I(p,d) = S(s)T(l(p,s)) \frac{1}{f} \Par{1-T(fl(p,d))} \]
By putting everything together, we have that: \[ L(p,d) = S(d) T(l(p,d)) + S(s)T(l(p,s)) \frac{1}{f} \Par{1-T(fl(p,d))} \] This formula, together with the scattering coefficients \(\sigma=(0.1, 0.3, 0.7)\) for the red, green and blue channel, and a position \(p\) that is \(0.3\) units within the atmosphere, gives the following result. You might have to reset the shader for it to work properly.