## Simulating a day’s sky

A nice sky is very important for most games. However, creating a good-looking one is not that easy. Researchers of the University of Utah have found a good approximation for the sky during a day (Paper). In this article I am going to explain the fundamentals of this algorithm. The complete project can be downloaded at the end of the article.

Most important for the calculation of the sky color is the position of the sun. Positions – or better directions – in the sky are expressed as polar coordinates (azimuth, elevation). Azimuth is the horizontal angle: 0 is south, Pi/2 is west, Pi is north and -Pi/2 is east. The second coordinate, elevation, specifies the height above the horizon. If elevation is 0, the according direction is directly at the horizon; if elevation is Pi/2, the direction is above the observer. The opposite angle is the zenital distance (short: zenith) and can be calculated by `zenith = Pi/2 - elevation`. The following figure shows the relation of these angles:

## Calculating solar position

The position of the sun is determined by the geographical position of the observer and the date. The date influences the sun’s declination which describes the solar position in relation to the earth. This value can be calculated using the julian date. Another factor, called the solar time, is needed to calculate the sun’s position. A value of 0 specifies noon and +/- Pi specifies midnight. Given those values, the solar azimuth and zenital distance can be calculated as follows (assuming that the observer is at a standard meridian):

```const double latitude = 50 * Math.PI / 180;
//Calculate the julian day
var Y = Date.Year;
var M = Date.Month;
var firstDayOfMonth = new DateTime(Y, M, 1);
if(M == 1 || M == 2)
{
Y -= 1;
M += 12;
}
var D = (Date - firstDayOfMonth).Add(TimeSpan.FromDays(1)).TotalDays;
var julian = 365.25 * (Y + 4716) + 30.6001 * (M + 1) + D - 1524.5;

var solarTime = Math.PI/12*(Date.TimeOfDay.TotalHours-12
+ 0.17 * Math.Sin(4 * Math.PI * (julian - 80) / 373)
- 0.129 * Math.Sin(2 * Math.PI * (julian - 8) / 355));

var declination = 0.4093 * Math.Sin(2 * Math.PI * (julian - 81) / 368);
var zenith = Math.ACos(Math.Sin(latitude)*Math.Sin(declination) + Math.Cos(latitude)*Math.Cos(declination)*Math.Cos(solarTime));
var azimuth = solarTime;
```

## Creating geometry for the sky

Now we have the solar position, it’s time to draw the sky. In order to draw something, we need some geometry. Most applications use a sky box or a sky dome. However, there is no reason why we should use that much polygons. Instead we will use the simplest volumetric object, the three dimensional simplex, which is a tetrahedron that surrounds the observer. For this object is that simple, it can easily be created in the vertex shader. There is no need to send geometry to the GPU:

```PS_IN VertexShader(uint id : SV_VertexID)
{
//create a tetrahedron
PS_IN output;
output.world_pos = float4(
(float)((id &lt;&lt; 1) &amp; 2) - 1,
(float)((id +  1) &amp; 2) - 1,
(float)( id       &amp; 2) - 1,
1);

output.position = mul(output.world_pos, view_proj);

return output;
}
```

We can draw the tetrahedron by issuing a draw call for a triangle strip with six vertices setting the input layout to null. This will draw the following tetrahedron:

One thing we should have in mind is other geometry. The sky should be behind everything else and centered around the observer. If we draw the sky as the last object in the scene, we can utilize the z-buffer to draw only to spots where there is nothing drawn yet. We can use the following states for that:

```DepthStencilState DSS
{
DepthEnable	= TRUE;
DepthFunc = EQUAL;
StencilEnable = FALSE;
};

RasterizerState RS
{
DepthClipEnable = FALSE;
};
```

The DepthStencilState will ensure that only pixels whith a z-buffer value of 1 will be overwritten and the RasterizerState will ensure that the sky tetrahedron will not be clipped due to the projection matrix’ znear and zfar plane.

Centering the tetrahedron around the observer is pretty easy. When we send the view matrix to the GPU, we ensure that the last column (the translation part) is set to 0:

```view.M41= 0;
view.M42 = 0;
view.M43 = 0;
viewProj.SetMatrix(view * proj);
```

## Calculating sky color

Now we have the geometry and solar position, it is time to actually draw the sky. To determine the sky color we use Perez’ model. This model can be used to calculate the sky color in CIE Yxy color space. Each color channel is calculated separately and needs some parameters: the absolute zenital values and constants A through E. The latter can be chosen or calculated with the air’s turbidity T. Here is how to calculate the values:

```private void CalculateZenitalAbsolutes()
{
var Yz = (4.0453 * Turbidity - 4.9710) * Math.Tan((4.0 / 9 - Turbidity / 120f) * (Math.PI - 2 * SolarZenith)) - 0.2155 * Turbidity + 2.4192;
var Y0 = (4.0453 * Turbidity - 4.9710) * Math.Tan((4.0 / 9 - Turbidity / 120f) * (Math.PI)) - 0.2155 * Turbidity + 2.4192; ;
this.Yz = (float)(Yz / Y0);

var z3 = (float)Math.Pow(SolarZenith, 3);
var z2 = SolarZenith * SolarZenith;
var z = SolarZenith;
var T_vec = new Vector3(Turbidity * Turbidity, Turbidity, 1);

var x = new Vector3(
0.00166f * z3 - 0.00375f * z2 + 0.00209f * z,
-0.02903f * z3 + 0.06377f * z2 - 0.03202f * z + 0.00394f,
0.11693f * z3 - 0.21196f * z2 + 0.06052f * z + 0.25886f);
var xz = Vector3.Dot(T_vec, x);
this.xz = xz;

var y = new Vector3(
0.00275f * z3 - 0.00610f * z2 + 0.00317f * z,
-0.04214f * z3 + 0.08970f * z2 - 0.04153f * z + 0.00516f,
0.15346f * z3 - 0.26756f * z2 + 0.06670f * z + 0.26688f);
var yz = Vector3.Dot(T_vec, y);
this.yz = yz;
}

private void CalculateCoefficents()
{
var coeffsY = new Coeff();
coeffsY.A = 0.1787f * Turbidity - 1.4630f;
coeffsY.B = -0.3554f * Turbidity + 0.4275f;
coeffsY.C = -0.0227f * Turbidity + 5.3251f;
coeffsY.D = 0.1206f * Turbidity - 2.5771f;
coeffsY.E = -0.0670f * Turbidity + 0.3703f;
this.coeffsY = coeffsY;

var coeffsx = new Coeff();
coeffsx.A = -0.0193f * Turbidity - 0.2592f;
coeffsx.B = -0.0665f * Turbidity + 0.0008f;
coeffsx.C = -0.0004f * Turbidity + 0.2125f;
coeffsx.D = -0.0641f * Turbidity - 0.8989f;
coeffsx.E = -0.0033f * Turbidity + 0.0452f;
this.coeffsx = coeffsx;

var coeffsy = new Coeff();
coeffsy.A = -0.0167f * Turbidity - 0.2608f;
coeffsy.B = -0.0950f * Turbidity + 0.0092f;
coeffsy.C = -0.0079f * Turbidity + 0.2102f;
coeffsy.D = -0.0441f * Turbidity - 1.6537f;
coeffsy.E = -0.0109f * Turbidity + 0.0529f;
this.coeffsy = coeffsy;
}
```

Most important for the calculation of the zenital Y value is the normalization with Y0. This is the absolute value if the sun was at zenith, which is the highest reachable value. By dividing by that value we ensure that Yz is in the valid range [0, 1].

Now we have all parameters and can switch over to the pixel shader which will calculate the actual sky color. It uses Perez’ function as defined in the paper:

```float Perez(float zenith, float gamma, Coeff coeffs)
{
return	(1 + coeffs.A*exp(coeffs.B/cos(zenith)))*
(1 + coeffs.C*exp(coeffs.D*gamma)+coeffs.E*pow(cos(gamma),2));
}

float4 RGB(float Y, float x, float y)
{
float X = x/y*Y;
float Z = (1-x-y)/y*Y;
float4 rgb;
rgb.a = 1;
rgb.r =  3.2406f * X - 1.5372f * Y - 0.4986f * Z;
rgb.g = -0.9689f * X + 1.8758f * Y + 0.0415f * Z;
rgb.b =  0.0557f * X - 0.2040f * Y + 1.0570f * Z;
return rgb;
}

float Gamma(float zenith, float azimuth)
{
return acos(sin(solar_zenith)*sin(zenith)*cos(azimuth-solar_azimuth)+cos(solar_zenith)*cos(zenith));
}

float4 PixelShader(PS_IN p, out float depth : SV_DEPTH) : SV_TARGET
{
depth = 1;

float azimuth = atan2(p.world_pos.x, p.world_pos.z);
float zenith = acos(p.world_pos.y / length(p.world_pos.xyz));

float gamma = Gamma(zenith, azimuth);
zenith = min(zenith, 3.1415926f/2.0f);
float Yp = Yz * Perez(zenith, gamma, coeffsY) / Perez(0, solar_zenith, coeffsY);
float xp = xz * Perez(zenith, gamma, coeffsx) / Perez(0, solar_zenith, coeffsx);
float yp = yz * Perez(zenith, gamma, coeffsy) / Perez(0, solar_zenith, coeffsy);

return RGB(Yp, xp, yp);
}
```

The function RGB converts colors from CIE Yxy space to RGB space. The function Gamma calculates the angle between the sun and an observed direction. Azimuth and zenith of an observed direction can be easily computed from a pixel’s world position with the given formulas. We explicitly set a depth value of 1 to ensure that the sky is always in the background.

## Results

The resulting program contains some more features than described above. It renders a very simple terrain with a solid color and directional lighting based on the solar position. The described algorithm is not appropriate for night skies. Therefore, some values are clamped (e.g. solar zenith) and blended (e.g. light color) to give a nice visual. The light color during day is calculated with the sky color formula, providing the solar position as the observed direction. The camera can be controlled with the mouse (including buttons for movement).

Please be aware that the sample is not focussed on terrain or the night sky.

Executables (1 MB)

Visual Studio 2012 solution (C#, SlimDX, 14 KB)

1. #1 by racarate on June 4, 2014 - 10:23 pm

this looks slick, cant wait to try it… i find it hard to implement straight from the papers, so i appreciate layman writeups like this

2. #2 by nicoschertler on June 4, 2014 - 11:25 pm

Thanks for the comment. I’m always glad to shed some light in the dark š

3. #3 by Carlos Oliva on July 21, 2014 - 2:56 am

Hi Nico, I got through your post from your question on gamedev.stackexchange. I’m trying to achieve more or less the same thing but in a 2D environment. I already made my implementation of the Preetham/Perez model but for some reason I’m getting some wacky colors. If you don’t mind taking a look, I published all of the details on a gamedev.net forum entry here: http://www.gamedev.net/topic/658937-preetham-sky-for-2d-wicked-colors-sample-images-and-video/

4. #4 by dekurian on April 7, 2016 - 6:27 am

Looks pretty good. Really would like to try it, but I’m not able to download the Sample-Code, since the Links seem to be dead. Is there any chance that you might fix this? Thanks !

5. #5 by nicoschertler on April 7, 2016 - 7:26 am