Calculating sunrise and sunset times

ebvalaim 2015-12-22 0 Comments

Today is the winter solstice - the shortest day in the year and the longest night. It is not a well-known fact, though, that although the days will only be longer now, the Sun will still rise a bit later. It is caused by the shape of the Earth's orbit, which is not perfectly round, but elliptical. The Earth will reach its closest point to the Sun in a little while, and it also moves faster than usual because of that. This in turn delays the solar noon by a few seconds each day, causing the times of sunrise and sunset to be later and later.

When is the latest sunrise and the earliest sunset, then, if not on the day of the solstice? I could probably check somewhere on the internet, but why should I, if I can calculate it myself with my computer ;)

I chose Haskell for the task - mostly because I still don't grok it, and I think it is a very interesting language that changes the way one thinks. I decided then to exercise it a bit.

First step - finding appropriate equations

Although many of the equations can be derived from first principles, some of them come from the observational evidence, so using the internet is a must. For example, everybody knows that the vernal equinox comes on the 21st of March (which, by the way, isn't true since about 2000 - now it's 20th of March, it will come back to 21st around 2100), but it is less obvious how far from the Sun the Earth is then and how fast it moves. This knowledge would suffice for calculations to a precision of a few minutes, but since we want such details as the "movement" of the solar noon, it has to be done better.

Finding the data isn't hard - Google gives the link to the article on Wikipedia at once. Almost everything is written there step by step, but of course it is nice to know what one actually calculates ;) What's especially important are the notions of right ascension (RA) and declination.

A long time ago people already had the idea that a system of describing the positions of stars and planets on the sky would be useful. The problem is that the sky rotates continuously (which, as we know, is caused by the Earth's rotation). Thus, an idea was born to introduce a celestial coordinate system independent from this movement, standing still with respect to the stars. It is quite easy, actually.

First, we introduce an analog of the geographical latitude, called declination. The declination equal to +90° corresponds to the celestial northern pole, i.e. the point at which the axis of Earth's rotation is aiming at the North Pole. Declination -90° is the celestial southern pole, 0° is the projection of the Earth's equator on the celestial sphere. This takes care of one coordinate.

The second one is right ascension, or RA for short. Just as we have an arbitrarily chosen zero meridian in Greenwich on Earth, the RA equal 0 corresponds to an arbitrary point, the so-called First Point of Aries. It is the point where the ecliptic (the plane of the Earth's orbit, or - as seen from Earth - the path the Sun traces on the sky during the year) crosses the celestial equator. The Sun crosses this point at the moment of the vernal equinox. On the opposite side of the sky we have the First Point of Libra, the second crossing of the ecliptic with the equator.

For a change, RA is measured not in degrees, but in hours, minutes and seconds and it takes the values from 0h to 24h. It is actually reasonable - the RA means the time that passes until a given point on the equator will come to the same place on the sky where the First Point of Aries is at that moment.

As I already mentioned, the movement of the Sun along the ecliptic is a reflection of the orbital movement of the Earth. Knowing the celestial coordinates (RA and declination) of the Sun, the geographical coordinates of the place of observation and which celestial meridian passes through the highest point in the sky at the moment, we can calculate when the Sun will rise, and when it will set.

This is exactly the kind of data we can get from Wikipedia. After calculating the so-called Julian Date from the current date, we can easily get the ecliptical coordinates of the Sun ($\lambda$ and $\beta$ on the Wiki), which we then combine into RA ($\alpha$) and declination ($\delta$) using the Earth's obliquity angle ($\epsilon$). We get the Earth's distance from the Sun for free in a way ;)

What will be really useful in calculations of sunrises and sunsets is the solar elevation angle (which is 0 at sunrise/sunset). Wikipedia links us to the articles about solar zenith angle and solar azimuth angle for that. Simply put, these angles define a spherical coordinate system in which the horizon is the equator, the zenith is the north pole, the nadir - the south pole, and the north points in the direction of the zero meridian.

Of course all of that is impossible to calculate if we don't know the Earth's orientation in space. Notions such as the hour angle and the sidereal time appear here. The sidereal time is just the RA of the meridian which is highest in the sky at the moment - if it is the meridian of the First Point of Aries, we have sidereal time 0. There is the Greenwich Sidereal Time, which can be calculated from the date and time, and local sidereal time, which is just Greenwich Sidereal Time shifted by the value dependent on the geographical longitude of the point of observation. When we add RA of the given object to the mix (in our case - the Sun), we get the hour angle for that object.

The hour angle tells us how far the meridian passing through our object is from the highest point in the sky. Knowing it and the declination of our object, we can calculate the azimuth and elevation (which only takes a few rotations on the celestial sphere). The complete equations are written on the Wiki (although when I found them, there was a little mistake in them caused by a mix in terminology - I'll take credit here and say that I detected and corrected it ;)).

When we are able to calculate the solar elevation above the horizon, we are only a step away from sunrises and sunsets - we only have to find the moments at which the elevations is 0.

Second step - coding

So, now that we know the theory, we can start coding. The complete program is here.

Initially, I started by creating my own datatypes for storing dates and times. Good programmers know though that you shouldn't roll your own crypto, nor your own date handling code, so I decided to look for some ready solutions and I found Data.Time.

The first step - the calculation of the Julian Date - is almost taken care of for us, except that the module calculates the integer part of the Modified Julian Date. The fractional part, which is just the time since midnight divided by 24 hours, is easily added, as is the correction to the Julian Date, or, what's more interesting, to J2000 epoch - which is just subtracting a constant.

The following functions calculate the parameters of the Earth's orbital movement and the ecliptical longitude (the solar ecliptical latitude is always 0, which makes things a bit easier):

All functions take the time calculated in the J2000 epoch as the parameter.

Now we can calculate the RA and declination of the Sun:

The last thing we can get without knowing the location of the observer is the Greenwich Sidereal Time:

Now we need the location, so we could use a type to store it:

Now we can calculate the local sidereal time and the Sun's hour angle:

And from here it is a straight path to the zenith angle (and, which is almost the same, the elevation) and the azimuth:

Now, how do we check when the elevation is 0?

At first I used a naive algorithm: for a given date, choose the noon and midnight (hours 0 and 12), and then check if at one of those moments the Sun is below the horizon, and above the horizon at the other one. If not, we have a polar day/night. If yes, divide the interval in two and take the two points, at which the Sun is again once below, and once above the horizon. Repeat until the time difference falls below a given threshold and assume that one of the moments is the moment we are looking for.

It works nicely... as long as the location isn't too far from Greenwich and there is no polar day/night. In other cases, not so good.

Above all, the longitude affects the time of the solar noon quite a bit. Since we calculate everything in UTC, it doesn't matter near Greenwich, but go far enough and the actual noon can change so that hours 0 and 12 are both during the day or the night. And so we have a false polar day/night.

This was easy to correct (just shift the time by the longitude / 15), but that's not all. The second problem is less obvious and has to do with the elliptical shape of the Earth's orbit. You see, the solar noon isn't always at 12 local time; it can go back and forth by up to about 15 minutes. This can cause the 12 o'clock to fall for example during the night near the polar night, even though the Sun was rising that day. The program recognized that as the polar night, though.

Fortunately, the working solution isn't much more complicated. It is sufficient to find the solar noon first, i.e. the moment of the maximum solar elevation. This can be done with a similar algorithm - we start from local 11:30 and 12:30 (because the noon doesn't change by more than ~15 minutes, we don't have to check the entire day and risk finding the noon of another day). Then we calculate the middle point and substitute it for the one in the pair that has the lower elevation. By repeating that, we find the solar noon.

Having the solar noon, we can choose moments 12 hours before and after it, and thus go back to the initial algorithm, but without the risk of missing a day or a night. This way we can find the sunrise and sunset every time.

The calculation of the length of the day is now only the matter of subtracting the sunrise from the sunset - no need to go into details ;)

And here is the complete algorithm

Third step - an application and the answer to the initial question

All that's left is to create a simple application that would allow us to use the functions above. I used optparse-applicative for that - this library allows for easy creation of command line parsers. Using the application, we can now answer the question of when the Sun rises the latest, and when it sets the earliest:

For Warsaw, Poland, as it turns out, the latest sunrise is on 30th of December, and the earliest sunset - on 13/14th of December :)