Writing a Day/Night Simulation

Something I’ve been itching to implement in my game is a day/night cycle. Over the last week, I did exactly that! I could have done it in under an hour by just adding a DirectionalLight2D and scripting it to change its direction over a fixed time period representing the movement of the sun. I went way further. I did an unreasonably accurate simulation of both the sun and moon, accounting for orbital and latitudinal effects. It has summer and winter solstices! It has 24-hour days and nights near the poles during the equinoxes! The sun moves north and south in accordance with the time of year! It even accurately simulates the phases of the moon!

Ok, “simulation” is not the right word. Simulation implies time steps and numerical integration. What I did was find the exact equations where I could plug in all the parameters and get out the positions of the sun and moon. So today I want to go over how I did this.

I could have initially worked directly from the equations on wikipedia for the Solar zenith angle and the Solar azimuth angle, but I find that simply using equations doesn’t give a great intuitive understanding of what they mean. I actually initially wanted to do this so that I could gain a better intuitive understanding of orbital mechanics from the surface of the Earth. I noticed some time back that the sun appears to rise from the northeast during the summer, even far north of the tropics, which feels extremely counterintuitive. I had a logical understanding to why (it surely has to do with the earth being tilted), but the intuition was not there.

The method that ended up working best was to make use of the subsolar point as described on the above wikipedia pages. The concept of the subsolar point is self-explanatory - it’s the point on the earth where the sun is directly overhead - so once I saw that there is a derivation starting from that, I backed away to experiment with doing it myself. My research had given me before this point knowledge of the concept of the sun’s declination and the hour angle, so I would work from those three concepts.

Our goal is to find the direction vector for the sun at a given latitude and time of day.

First off, we define mathematically where the subsolar point is. For any given latitude theta and longitude phi with a defined radius r, the cartesian position is

$$\hat{r} = \left( r \cos{\theta} \cos{\phi}, r \cos{\theta} \sin{\phi}, r \sin{\theta} \right)$$

We have here the equator on the xy-plane and the north pole along +z. We are using a right handed coordinate system for convenience, with +x pointing towards us and +y pointing to the right. Finally, we are also going to be working on the unit sphere so our direction vector is normalized.

We begin with the subsolar point. Orient our sun direction vector a longitude of 0. The subsolar point will be at latitude equal to the declination delta.

$$\hat{r}_S = \left( \cos{\delta}, 0, \sin{\delta} \right)$$

As noted before, the north pole is along +z. That makes the axis of rotation for the planet be the z-axis. This actually makes things rather convenient for us. Let t be the number of hours past noon for our target location. We can then define an “hour angle”

$$h = \frac{t}{24} 2\pi = \frac{\pi t}{12}$$

This hour angle is effectively the longitude of our target coordinate relative to the subsolar point. Hence we can write our target coordinate as

$$\hat{r}_T = \left( \cos{\theta} \cos{h}, \cos{\theta} \sin{h}, \sin{\theta} \right)$$

To be clear on coordinates here - positive longitude is moving east, negative longitude is moving west. Relative to the sun which appears to move west across the sky, a higher hour angle would be further east. Hence a positive hour angle would correspond with a positive longitude.

Our next step is to rotate our target coordinate into the xz plane (0 longitude) around the axis of rotation and see where that puts the subsolar point. The relevant rotation matrix is

$$R_z = \left[ \begin{matrix} \cos{h} & \sin{h} & 0 \\ -\sin{h} & \cos{h} & 0 \\ 0 & 0 & 1 \end{matrix} \right]$$

Giving our solar vector at latitude 0 (the equator) and our correct longitude as

$$\hat{r}_S = \left( \cos{\delta} \cos{h}, -\cos{\delta} \sin{h}, \sin{\delta} \right)$$

With the new location of our target coordinates in this reference frame as

$$\hat{r}_T = \left( \cos{\theta}, 0, \sin{\theta} \right)$$

Just one last step. We want to rotate the target point to be along the +x axis and see where that puts the subsolar point. This will give us a coordinate system where +x is radial out from the planet, and the yz plane represents north, east, south, and west, with +z pointing north and +y pointing east. Our rotation matrix is

$$R_y = \left[ \begin{matrix} \cos{\theta} & 0 & \sin{\theta} \\ 0 & 1 & 0 \\ -\sin{\theta} & 0 & \cos{\theta} \end{matrix} \right]$$

Our target coordinate goes, as expected, to

$$\hat{r}_T = \left( 1, 0, 0 \right)$$

And our subsolar point goes to

$$\hat{r}_S = \left( \cos{\theta} \cos{\delta} \cos{h} + \sin{\theta} \sin{\delta}, -\cos{\delta} \sin{h}, \cos{\theta} \sin{\delta} - \sin{\theta} \cos{\delta} \cos{h} \right)$$

Now that we have these equations, we need to define our solar zenith and azimuth angles in this reference frame. This is thankfully fairly simple - the zenith angle is our angle away from directly overhead, so a vector entirely in the +x direction would have a zenith of 0. Meanwhile, the azimuth is our angle clockwise from north. With +z pointing north and +y pointing east, the solar position equation using zenith and azimuth is

$$\hat{r}_S = \left( \cos{\Theta}, \sin{\Theta} \sin{\Phi}, \sin{\Theta} \cos{\Phi} \right)$$

Equating the two representations, we have

$$\cos{\Theta} = \cos{\theta} \cos{\delta} \cos{h} + \sin{\theta} \sin{\delta}$$

$$\sin{\Theta} \sin{\Phi} = -\cos{\delta} \sin{h}$$

$$\sin{\Theta} \cos{\Phi} = \cos{\theta} \sin{\delta} - \sin{\theta} \cos{\delta} \cos{h}$$

The last thing we need to have fully determined our equations is to define the solar declination. This is one that thankfully does make at least some intuitive sense. The declination can be defined as the zenith of the sun at local solar noon at the equator. If the earth’s tilt were zero, the declination would also always be zero. The earth’s tilt is not zero! This gives us seasons as the sun appears to “drift” north and south over the course of a year. During the equinoxes, the declination is zero as those are the points where the sun “crosses” the equator. At the solstices the sun is at its most extreme latitudes equal to the earth’s tilt. The expression for this is fairly simple.

$$\sin{\delta} = \sin{\Phi_T} \sin{EL}$$

Here phi is the earth’s tilt while EL is the “ecliptic longitude” - effectively where in its orbit the earth is. Because the earth’s orbit is eccentric (non-circular), this ecliptic longitude does not map directly to the time of year. There are small deviations which are most extreme, coincidentally, around the solstices.

With the equations in hand, next is implementing this in Godot! It ends up being pretty simple thanks to DirectionalLight2D. Simply set the height of the light to the x component of our solar position vector (and turn the light off when the x component is negative as that indicates the sun has set), use the other two components to calculate the rotation of the light (keeping in mind that Godot’s DirectionalLight2D’s rotation directly maps to the azimuth!), and set variables appropriately.

We simplify the ecliptic longitude by assuming no orbital eccentricity so it can be calculated directly as

$$EL = \frac{t_{days}}{T_{year}} 2\pi$$

The axial tilt of our planet can be set freely, and with the ecliptic longitude we have the declination. Note that we constrain the declination to be

$$-\frac{\pi}{2} \le \delta \le \frac{\pi}{2}$$

as it is effectively the latitude of the sun, and latitudes have the same range.

The latitude is freely set by the location you’re running your day/night cycle at, and the hour angle is naturally set to the passage of the clock. Longitudes can be accounted for by adding a “time zone” offset to the hour angle.

The last thing I added was a moon simulation. It is… nearly identical to the sun simulation. Instead of calculating the ecliptic longitude on the length of a year, I used the length of a month. There was a small corrective factor I needed to add where I converted from solar time to sidereal time, but this was more to force the moon to “drift” over the year for stylistic reasons. The correction I used was

$$t' = \frac{T_{year}-1}{T_{year}} t$$

but again this is honestly just a stylistic thing. I equivalently could have just changed the month length so it wasn’t an even divisor of the year.

I also needed to calculate the relative brightness of the moon at different phases. I was able to intuitively guess an equation of

$$B = \frac{1 - \cos{\theta}}{2}$$

where theta is the angle between sun and moon from the surface of the planet with the cosine calculated by taking the dot product of the solar and lunar direction vectors. My logic for this was - when the sun and moon are aligned, the brightness should be zero because we’re in a new moon. When the sun and moon are on opposite sides of the planet, the brightness should be 1 because we’re in a full moon. The variation should likely follow a sine or cosine curve as everything in orbital mechanics is circles and spheres, therefore we’re using trigonometry. Since we’re hitting maxima and minima at angles of 0 and pi, it’s going to be cosines specifically. The equation naturally falls out from this, and it turned out to be the correct one.

Beyond that, there were a few other small stylistic things I added. When the sun’s position is low enough in the sky (just looking at the raw x coordinate), I do a lerp on the sun’s color to simulate sunset/rayleigh scattering becoming more intense. I also adjust the brightness of the sun and moon as they eclipse the horizon so they naturally fade to nothing until they are completely blocked by the planet.

It ends up using around 12 trig calls every time I update the sun and moon positions which might be a little on the expensive side. If it ever becomes a performance issue, I can always reduce how frequently the positions are updated and possibly just lerp or slerp between fully calculated positions on the non-update frames.

It all honestly came out looking really good! A little timelapse of my testing environment is attached.