Introduction
If you are living in the northern hemisphere then you surely enjoy that now, at the beginning of February, the days have become noticeably longer again. The shortest day was just a few weeks ago (December 21) but the sun already sets nearly an hour later than on this dark day!
I guess most people would think that the shortest day of the year is also the day of latest sunrise and earliest sunset. That sounds most natural. I also did so until I looked at rise and set times in more detail. And to my big surprise, the earliest sunset happens on December 12 whereas the latest sunrise happens on December 30!
The reason for this is simple, but not obvious: the days in December are 30 seconds longer than 24 hours! This phenomenon is called 'Equation Of Time' and there are many articles about it in the world wide web, but unfortunately none of them helped me to really understand the "why" and so I wrote this 3D Sun-Earth simulation which calculates rise and set times depending on various parameters.
Background
At first glance, it might be surprising that a day in December - or any other day - does not exactly last 24 hours. The length of a day should be determined by the speed of rotation of the earth around itself and that shouldn't change throughout the year!
In fact, the speed of rotation of the earth around itself is very stable and it takes about 23 hours and 56 minutes for the earth to complete a full rotation. But the earth has to rotate a little bit more until a day is finished, because it also travels around the sun about one degree during that time. You can see this effect in the following picture from Wikipedia:
In fact, it is this extra rotation which causes the problem! Depending on the earth's location and orientation, it sometimes takes more than 4 minutes to finish a solar day and sometimes less.
When our ancestors decided to split the time between two consecutive noons into 24 hours, they didn't have atomic clocks or things like that. I guess they even didn't care about minutes or seconds. If they had to arrange a meeting, it probably was something like "let's meet at the tenth hour" or so. Happy times then!
What If the Days in December were Exactly 24 Hours?
Well, the world would be as we'd expect: at the beginning of December, the sun would rise a minute later than the day before and set a minute earlier. On December 12, the difference would have been reduced to half a minute and on December 21, the shortest day, there would be no difference anymore between the rise and set times.
From that day on, the differences grow again, but in the opposite direction. So at the end of December, the sun would rise half a minute earlier and set half a minute later than the day before.
What's the Effect of the Additional 30 Seconds?
First of all, it means that the sun reaches a certain position at the sky 30 seconds later than the day before. So if the sun happened to be in the exact south direction at 12:00:00 on a certain day, it would reach the same position at 12:00:30 on the next day. And the same is true for sun rise and sun set. Everything is shifted by the same amount of 30 seconds.
If we apply this shift to the differences that we saw in the previous section, we come to this:
At the beginning of December, the sun rises 1.5 minutes later than the day before but sets 0.5 minutes earlier.
On December 12, it rises one minute later, but does not set earlier.
On December 21, it rises half a minute later and sets half a minute later (so no net change).
At the end of December, it rises no earlier but sets one minute later.
So Let's Celebrate on December 13!
Because this is the first day when the sun sets later than the day before! Even if the difference is just 2 seconds, it is definitely the beginning of a new, brighter era.
But What If Every Day is Longer than 24 Hours?
Well, that would be bad! That would mean that sooner or later the sun's position would be completely out of sync with our clocks. So when the sun is in the south position at noon, our clocks would tell something about afternoon, evening or even night! Since we do not experience such a strange thing, there must be something which compensates for the longer days in December.
And in fact, in spring and in autumn, the days are shorter than 24 hours, whereas in summer they are slightly longer again!
In the end, if we look at the phenomenon for longer than a year, we will find that the time between two consecutive noons varies periodically about an average value. And this value is exactly 24 hours!
Now isn't it magic that something in the universe cares about our clocks and their 24 hours?
For sure it's not, because it's not the earth which cares about that time, but we defined 24 hours to be the average time span between two noons!
So Why Isn't Each Day Exactly 24 Hours?
There are two reasons:
- The earth's orbit around the sun is not a perfect circle but an ellipse
- The earth's rotation axis is not perpendicular to the plane of its orbit
According to Kepler's laws of planetary motion, the orbit of a planet is an ellipse with the sun at one focus. That means that there are times when the planet is nearer to the sun than at other times. And that also means that the planet moves faster around the sun in these times.
If you look at the picture at the beginning of this article, it gets clear that if a planet moves to a position beyond position 2 after one rotation around itself, it needs to do more extra rotation until the sun is at noon position again.
The earth reaches its nearest position to the sun on January 3 (despite the fact that it's hotter in summer than in winter). If the earth's rotation axis were perpendicular to its orbit, that would cause the days to be 8 seconds longer than 24 hours at that time. But in summer, when the earth is further away from the sun, the days would be 8 seconds shorter, and in between there are two times in spring and in autumn where each day takes exactly 24 hours.
That's for the easy part. The earth's speed around the sun varies periodically with a maximum in winter and a minimum in summer and that makes the days longer and shorter.
The not-so-easy part is the effect of the earth's so called obliquity or axial tilt, the fact that the earth's rotation axis is not perpendicular to its orbit around the sun. There is an angle of about 23 degrees which most prominently causes our seasons but also causes our noon-to-noon intervals to vary with two maxima and two minima per year.
If we neglected the eccentricity and assumed that the earth's orbit is a perfect circle, so its speed around the sun were the same throughout the year, then the days would be 20 seconds longer than 24 hours in winter and in summer, and they would be 20 seconds shorter in spring and in autumn depending on this obliquity.
All in all, these two effects sum up in something we can see here, again from Wikipedia:
The second contribution to the Equation of Time (the effect of obliquity) is really hard to explain in words and that's where the simulation comes in.
Simulation of the Sun-Earth System
Simulations are great! They help to understand things that are either too big, too small, too fast or too slow for our human dimensions. In case of our problem, the thing is too big (the distance between the sun and the earth is about 150 million kilometers) and too slow (we have to wait 24 hours to measure the difference between two rise times of the sun).
But when doing a simulation, we can choose our model to fit in our world (let's say the screen of our computer) and we can choose how fast time flies (e.g. 24 hours in ten seconds). We can play around with parameters that have an influence on what we are studying and we can repeat the same simulation over and over again until we get it.
But the simulation will only give us true answers, if the models we use are correct. So it's crucial to start with the right model. It doesn't help if the model is too complex but for sure it's also bad if it's too simple. Something in between sounds right. But where exactly is that?
What Do We Want to Achieve?
We want our simulation to show the times for sun rise, noon and sun set for every day of the year with an accuracy of at least one second. The following parameters need to be adjustable: the obliquity of the earth's rotation axis, the eccentricity of the earth's orbit and the latitude of the observer on earth.
We do not care if the abolute values match real values (there are web sites which can be used to calculate sun rise and sun set times, for e.g. London or any other place on earth). But it's important that the difference between two values is like the real difference.
The Earth Model
Although some people still think that the Earth is a flat disc, let's assume it's a nearly perfect ball which travels around the Sun on a nearly circular orbit. So let's start with a perfect ball for the Earth and a perfect circle for the orbit. We just have to make sure that there is a small eccentricity in the orbit, which ends up in a variable speed throughout the year.
The Earth is spinning around itself in a very constant way. It takes 23 hours, 56 minutes and 4.1 seconds to complete a rotation, so the angular velocity (in radians per second) is this:
double wEarth = MathUtils.PIx2 / MathUtils.ToSeconds(0, 23, 56, 4.1);
We can also specify a mean angular velocity for the Earth orbiting the Sun: it takes 365 days, 6 hours, 9 minutes and 9.54 seconds to complete a revolution. In other words:
double wSun = MathUtils.PIx2 / MathUtils.ToSeconds(365, 6, 9, 9.54);
Note that wSun does not mean the angular velocity of the Sun spinning around itself, but means the mean angular velocity of the Earth orbiting the Sun! If we add the following declarations...
public double Latitude = 51;
public double Obliquity = 23.44;
public double Eccentricity = 0.0167;
double oneDay = MathUtils.ToSeconds(1, 0, 0, 0);
...we can write a method which calculates the Earth's position and rotation based on a time value (in seconds) like this:
public void Update()
{
double angle = MathUtils.ToDegrees(wEarth * time);
EarthRotation = new Quaternion(Math3D.UnitZ, angle);
AxialTilt = new Quaternion(Math3D.UnitY, -Obliquity);
angle = wSun * time;
EarthAngle = angle - 2 * Eccentricity * Math.Sin(angle - 15 * wSun * oneDay);
EarthPosition = new Point3D(3 * Math.Cos(EarthAngle), 3 * Math.Sin(EarthAngle), 0);
}
The first quaternion (EarthRotation) shows how much the Earth turned around itself based on the time, whereas the second one (AxialTilt) just shows the constant obliquity.
The calculation of the Earth's angle on its way around the Sun is a two-step process: first the angle is calculated based on the mean angular velocity wSun. Then a correction to that angle is made. This correction takes into account that the orbit is an ellipse with a perihelion at January 3 and an aphelion half a year later. This approximation is valid for small eccentricities and is taken from here.
Now that EarthAngle is known, we can calculate the Earth's position in cartesian coordinates. We are free to choose a coordinate system that makes our calculations easy, so we choose the xy plane to be plane of the ecliptic, i.e. the plane of the Earth's orbit. The Sun will be located at the origin.
The factor of 3 in the calculation of the Earth position is somewhat arbitrary. It has no effect on rise or set times. I've chosen this value to make the 3D model of the moving Earth (which is a sphere of radius 1) look OK.
About the Coordinate Axes
If we choose a time value of 0 in the above Update() method, all resulting angles will be 0 and the Earth will be located at (3, 0, 0), so somewhere at the positive x axis. The Earth's rotation axis is always tilted by -23 degrees around the y axis. So at the time of 0, the north pole is tilted exactly towards the origin, i.e. the Sun. This is exactly the summer condition.
So the x axis of our global coordinate system is the winter-summer axis, whereas the y axis is the autumn-spring axis.
Setting the Observer Location
Sun rise and sun set times depend decisively on the latitude of the observer. If you've ever been in the tropics near the equator, you surely noticed that the sun rises and sets nearly always at the same time: 6 o'clock in the morning for sun rise and 6 o'clock in the evening for sun set. Every day. Throughout the year!
In contrast, people living in the polar regions have non-stopping daylight in summer (the sun never sets) and non-stopping darkness in winter (the sun never rises).
We will start the simulation with an observer at a latitude of 51 degrees, which could be London. We don't really care about the longitude of our observer. For sure, it should be constant and we choose it in a way that when time is zero (so the Earth is in the exact summer position), the local time of our observer is midnight.
This is how the situation looks like for time = 0:
The Earth is located somewhere at the positive red x axis, its rotation axis is tilted towards the Sun and the observer (the blue spot at the parallel of 51°) is facing the night.
And Yet It Moves!
When the simulation is started, a loop is entered in a background thread which increases the time value by 1/100 second (this value is adjustable) and then calls the above Update() method. If a second of the simulated time has elapsed, the code checks four conditions:
- Did the Sun rise above the horizon?
- Did the Sun reach noon position?
- Did the Sun fall below the horizon?
- Did the Sun reach midnight position?
If one of these conditions is true, the observed time value is stored for comparison with the corresponding value of the next day.
While the background thread is updating the core state values of the earth as fast as it can, the UI thread picks up these values within a 30 millisecond timer method and updates the position and rotation properties of the corresponding 3D model objects.
The 30 millisecond UI update should be fast enough to ensure a smooth UI experience, but it also should be slow enough to give our computer enough resources to perform the steps which are required to check for the above four conditions.
The Math
The answer to each of the above questions cannot be given in one shot. Since our simulation is running at certain time steps, we will not be able to identify the exact time when for example the sun hits the horizon. We will only be able to identify that the sun was above the horizon at step n and that it is below the horizon at step n + 1. But if we know those two time values, we are able to calculate the exact hit time by doing simple linear interpolation.
The first thing we need to know is the observer position in global coordinates at a certain time. When we have this, we can draw an imaginary line to the Sun's position (which is always the origin of our coordinate system) and decide, if this is below or above the horizon or if it's left or right of the so called meridian (the south direction which marks the noon position).
The Observer Position in Global Coordinates
The calculation of the observer position is very easy if we take the same approach that 3D libraries do: use transformation matrices! If we knew the transformation matrix of the Earth, we can transform points from the Earth's coordinate system into the global coordinate system simply by calling earthMatrix.Transform().
Building up this matrix again is very easy, because the above Update() method provides everything we need. Before doing so in a formal way, let's play with an imaginary construction kit. Our Earth is a small ball with radius 1 which is currently located at the center of our global coordinate system. We need to apply three transformations so that this ball represents the Earth at a certain time:
- Rotate the ball about the z axis to show the correct spin state
- Rotate the ball about the y axis to show the correct axial tilt
- Translate the ball to the correct xy coordinates
In code, it looks like this:
Matrix3D earthMatrix = new Matrix3D();
earthMatrix.Rotate(EarthRotation);
earthMatrix.Rotate(AxialTilt);
earthMatrix.Translate((Vector3D)EarthPosition);
Remember that the Update() method calculates the required quantities EarthRotation, AxialTilt and EarthPosition. With this matrix, we can transform points from the Earth system to the global system by calling earthMatrix.Transform(Point3D pt).
So the only thing that's left to do is to calculate the observer position in the Earth's coordinate system. We can do this with another transformation matrix. Let's come back to our imaginary construction kit again. What would be the transformations that we need to apply to an observer who is currently located at the center of the Earth?
- Translate the observer location to the surface of the Earth
- Rotate the observer location to the required latitude
Since our model of the Earth has a radius of 1, the code looks like this:
Matrix3D locationMatrix = new Matrix3D();
locationMatrix.Translate(new Vector3D(1, 0, 0));
locationMatrix.Rotate(new Quaternion(Math3D.UnitY, -Latitude));
If we combine this matrix with the transformation matrix of the Earth (by simply multiplying them) we get a matrix that transforms points from the observer coordinate system to the global coordinate system. Since our observer is sitting at the center of its own system, i.e. at the point (0,0,0), its global coordinates are:
Point3D location = (locationMatrix * earthMatrix).Transform(new Point3D(0, 0, 0));
So now that we know the observer location in global coordinates, we can perform the next steps.
Calculate the Sun's Position with Respect to the Horizon
To answer the question if the Sun is above or below the horizon, look at this picture:
The Earth position is marked with E, the observer is located at L. Z is a point in the sky which is straight above the head of the observer. It's the so called zenith. It's direction is perpendicular to the horizontal plane which is shown as a brown line denoted with H.
The zenith direction immediately helps us to identify if the sun is visible. If the angle between zenith direction and Sun direction is smaller than 90 degrees, the Sun is above the horizon and visible. Otherwise, if the angle is greater than 90 degrees, the sun is below the horizon and invisible.
So we just need to know the angle between the vector to the Sun and the vector to the zenith. This can be done manually by using the dot product of the two vectors, but the .NET libraries already come with some handy methods. But first, we need to calculate the two vectors!
The zenith direction is given by L - E and the Sun's direction is simply -E. The latter sounds wrong because in our model the direction to the Sun would be S - L, or simply -L because the Sun is located in the origin, but our model is too simple in this case.
The relation between Earth radius and distance to the Sun is very wrong in this model and if we calculate the direction to the Sun by S - L, we would do something wrong. In practice, we can assume that the Sun is infinitely far away (behind the origin) and that all rays from the Sun hit the Earth in parallel, which ends up in -E being the direction to the sun.
So the code for calculating the angle between Sun and zenith is this:
Vector3D dirSun = -(Vector3D)EarthPosition;
Vector3D zenith = location - EarthPosition;
double angle = Vector3D.AngleBetween(dirSun, zenith);
Identifying Sun Rise and Sun Set Times
Watching the angle between the Sun and the zenith during the simulation allows for identification of sun rise and sun set:
- If the angle was greater than 90° at the previous step and smaller than 90° at the current step, a sun rise did happen
- If the angle was smaller than 90° at the previous step and greater than 90° at the current step, a sun set did happen
In both cases, we want to find out when exactly the angle was 90°. This is a simple linear interpolation problem and since the time step between two of the above checks is exactly one second, the code to calculate the corrected time value can be simplified to this:
double corrTime = time - (angle - 90) / (angle - oldAngle);
where time is the actual time, angle is the actual angle and oldAngle is the angle from the step before.
So What About Noon?
To identity the time of noon, we need to be sure what noon means. More specifically what solar noon means.
In my simple words, it is solar noon when the sun crosses the exact south direction on its daily walk from its rising point somewhere in the east to its set point somewhere in the west. At least on our northern hemisphere.
For every observer on Earth, there is exactly one south direction and the sun is either to the left of this direction (before noon) or to the right of this direction (after noon). We can visualize this direction by an imaginary half-disc that we draw into the sky which starts at the north pole, goes through the observer's zenith and ends at the south pole.
OK, that might sound confusing, so take a look at this:
Hope that this half-disc gets clear now. The tiny green spot of the Earth's surface is our observer at 51°. Note that the half-disc moves around the Earth in the same way that the observer does (which is counter-clockwise when looking down to the north pole). We take our south direction with us while we spin around the polar axis! Also note that I did set the axial tilt in the above picture to zero just for a clearer illustration.
And this is my definition of noon: on noon the half-disc crosses the sun, or in other words, the sun crosses our half-disc. So it is noon, when the rotating half-disc hits the Sun's center.
The name of this half-disc is meridian which comes from the Latin word meridies which just means midday. And as a side note, this is the meaning of AM and PM in the Anglo-American countries: before noon, the Sun is left of the meridian or ante meridiem (AM) in Latin, and after noon, it is right of the meridian or post meridiem (PM).
You can see the meridian in action here:
Identifying the Time of Noon
Now that we are sure about the definition of noon, we are able to calculate the exact time of noon if we observe the following:
- The Sun is right of the meridian at the current step
- The Sun was left of the meridian at the previous step
Now what do left and right mean?
With a little help of Analytic Geometry, we will find that both left and right are talking about the distance of an object to a certain plane. Left means negative distance, right means positive distance. And the plane we are talking about is the meridian. So how can we calculate the apparent distance of the Sun to the meridian?
Again it is the vector's dot product which helps here. The first vector for sure is the direction to the Sun (which we already have) and the second vector is the normal vector of the meridian, i.e. a vector which is perpendicular to the meridian. We can see its direction in the second last picture: it is H - L.
If the Sun is at noon position, then H - L is perpendicular to the Sun's direction and the dot product is zero. In any other case, the dot product is either a positive number or a negative number.
To calculate the normal vector, we'll go back to the untransformed Earth. Remember that the untransformed Earth is just a sphere with radius 1 at the origin of the coordinate system with the polar axis showing into the z direction. The observer was located at (1, 0, 0) so the meridian in that case is just the xz plane. And hence the normal vector of this plane is just the y direction, i.e. the vector (0, 1, 0).
So for the calculation of the normal vector at a specific time value, we just have to find out what happened to the Earth's y direction after it has been transformed. But that's an easy thing to do because we already have the Earth's transformation matrix. If we know the transformed y direction, we have to subtract it from the Earth's current position and that's the desired normal vector.
So the code to calculate the distance of the Sun to the meridian is this:
Point3D earthUnitY = earthMatrix.Transform(new Point3D(0, 1, 0));
Vector3D normalOfMeridian = EarthPosition - earthUnitY;
double distance = Vector3D.DotProduct(normalOfMeridian, dirSun);
It doesn't matter if dirSun is a unit vector, because we will just watch out for a change of sign.
It turns out that the above distance also helps us to identify the time of midnight:
- If the distance was negative at the previous step and it is positive at the current step, solar noon did happen
- If the distance was positive at the previous step and it is negative at the current step, midnight did happen
Again, if one of the above happens, we can calculate the exact time of the event by linear interpolation:
double corrTime = time - distance / (distance - oldDistance);
with time being the current time, distance being the current distance and oldDistance being the distance at the previous step.
The Application
If you start the application and press the Start button, the Earth will move around the Sun in early December. If you wait for about a minute, it will be in the January position and a table to the right will show rise, noon and set times plus their differences in seconds for a certain day:
Pressing the Start button again stops the simulation. You can choose the starting day and month with the controls to the left of this button. The Noon button lets the simulation start at the specified time but it will stop at the next noon.
You can change the camera position at any time by using the combobox to the right of the Demo button. Some of these position are called "fixed" which means that you cannot modify the position interactively. "Free" positions allow you to control the camera's position and orientation with the mouse. The following positions are available:
Fix: Overview 1
The camera is located at a certain point in 3D space (0.5, -1, 12) and is pointed to the origin with an up vector that has a positive z value. In this view the Sun stands still and the Earth moves around the Sun counter-clockwise.
Fix: Earth Eclip. Pole
The camera is located straight above the Earth at a distance of 10 in the global z direction (which is the direction of the ecliptic pole). It is pointed downwards the z direction with the up vector in the y direction. Now the Earth stands still but spins around itself and the Sun moves around the Earth.
Fix: Earth North Pole
Now the camera is located straight above the north pole at a distance of 10 and is pointed down the pole axis with an up vector of the global y direction. Like before, the Earth stands still but spins and the Sun moves around the Earth.
Fix: Earth Location
This is an interesting view from straight above the observer down to the Earth. Consider the camera was mounted at a giant skyscraper at the observer's location and pointed down to the observer with an up vector showing to the north pole. In this view the Earths stands completely still, so even does not spin around itself, because the camera is spinning around in the same way. But now the Sun apparently moves around the Earth once a day!
Fix: Earth Antipodal
This is also a very interesting view from the position of the so called Counter-Earth, which is a hypothetical planet located always at the other side of the Sun from the Earth. I remember a film called Journey to the Far Side of the Sun which picks up this hypothesis. Well, that's a long time ago!
In this view, both Earth and Sun stand still (with the Earth spinning) but now we can see the remarkable apparent spinning of the Earth's polar axis as the Earth orbits the Sun. For sure, it is an essential characteristic of the Earth's axis that it is fixed with respect to the global coordinate system, but looking at the Earth from this position reveals much about what's going on.
Free: From Location
This is the view that we all know best. The camera is at the observer's location and pointing to the south direction. Now we can see the Sun moving across the sky as we see it in real life. The camera is not fixed, so you can turn around with left mouse dragging. You can also zoom in and out with the mouse wheel.
Free: Overview 2 and 3
Like in Overview 1, the camera is located at fixed positions in space. But this time, you can change the heading with left mouse dragging and the position with Ctrl + arrow keys (and Ctrl PgUp PgDn).
More Controls
In the middle of the control area, you'll find controls to change the eccentricity, obliquity, observer's latitude and the simulation speed. A value of 10 will set the time step to 1/100 second. This step can be increased/decreased by a factor of 2 with the spin control. You can even reverse time by pressing the Inv button:
The next set of controls is used to add or remove additional 3D models which might help to understand what's going on: the global coordinate system, the ecliptic, the observer's location, the horizon, the meridian and the shadow border:
Building Up the 3D Scene
The 3D scene is created with my small WFTools3D library, which you can find here on Codeproject and on Github. Depending on the chosen settings, it contains more or less 3D models but it will always show the Sun and the Earth.
Both Sun and Earth are of type WFTools3D.Sphere. The Sun has a radius of 0.2 and is always located at the origin, whereas the Earth has a radius of 1 and its position and rotation state get updated when the rendering timer is called:
void Update()
{
earth.LockUpdates(true);
earth.Position = simulator.EarthPosition;
earth.Rotation1 = simulator.EarthRotation;
earth.Rotation2 = simulator.AxialTilt;
earth.LockUpdates(false);
...
}
Both Rotation1 and Rotation2 are applied to the model's transformation matrix before the translation which corresponds to Position is applied. There is another property, Rotation3, which gets applied after the translation, but it's not used here.
Method LockUpdates(true) prevents the transformation matrix to be recalculated whenever a single 3D property is set. After all properties are set, LockUpdates(false) will update the matrix.
Many of the additional 3D models, e.g. observer location, pole axis, horizon or shadow border, are children of the Earth model. So they do not have to be updated in the timer method. They will follow the Earth while spinning around itself and on its way around the Sun!
The Demo Mode
The Demo mode is some kind of course which tries to explain the Equation of Time by text and with the help of the simulator. It starts with the rather simple explanation of the seasons on Earth and ends with the more complex influence of the axial tilt to the length of a solar day.
The parameters for the simulation, i.e. start time, observer latitude, obliquity and eccentricity are adjusted automatically to highlight certain situations. Once you've gone through the 27 pages, your understanding of the phenomenon should at least be improved. At least that's what I hope!
Find Real Sun Rise and Sun Set Times
If you want to find real sun rise and sun set times for your home town, visit calsky.com. You can enter your location, a start date and a duration and you'll get the rise, set and transit times (which means solar noon time).
The End
I know it was again a long journey, but I hope you liked it! Happy coding!
History
- 2-Feb-2017: Initial upload
- 20-Mar-2017: Updated the first picture and both the download files
