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 ( and on the Wiki), which we then combine into RA () and declination () using the Earth's obliquity angle (). 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):

1 2 3 4 5 6 7 8 9 10 11 |
meanLongitude :: Double -> Double meanLongitude n = 280.46 + 0.9856474*n meanAnomaly :: Double -> Double meanAnomaly n = 357.528 + 0.9856003*n -- ecliptic longitude eclipticLong :: Double -> Double eclipticLong n = let g = deg2rad $ meanAnomaly n in meanLongitude n + 1.915 * sin g + 0.02 * sin (2*g) |

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

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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
-- Earth's obliquity obliquity n = 23.439 - 0.0000004 * n -- Right Ascension of the Sun sunRA :: Double -> Double sunRA n | ra < 0 = ra * 12 / pi + 24 | otherwise = ra * 12 / pi where eps = deg2rad $ obliquity n l = deg2rad $ eclipticLong n x = cos eps * sin l y = cos l ra = atan2 x y -- Declination of the Sun sunDec :: Double -> Double sunDec n = let eps = deg2rad $ obliquity n l = deg2rad $ eclipticLong n in rad2deg $ asin (sin eps * sin l) |

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

1 2 3 |
-- Greenwich Mean Sidereal Time gmst :: Double -> Double gmst n = 18.697374558 + 24.06570982441908 * n |

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

1 |
data Location = Loc { lat :: Double, lon :: Double } |

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

1 2 3 4 5 6 7 |
-- local sidereal time lst :: Location -> Double -> Double lst l n = gmst n + lon l / 15 -- Sun's Local Hour Angle sunLHA :: Location -> Double -> Double sunLHA l n = (lst l n - sunRA n) * 15 |

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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 |
-- zenith angle sunZenith :: Location -> Double -> Double sunZenith l n = let d = deg2rad $ sunDec n h = deg2rad $ sunLHA l n lat' = deg2rad $ lat l in rad2deg $ acos (sin lat' * sin d + cos lat' * cos d * cos h) -- elevation sunElev :: Location -> Double -> Double sunElev l n = 90 - sunZenith l n -- azimuth sunAzim :: Location -> Double -> Double sunAzim l n | az < 0 = az + 360 | otherwise = az where lat' = deg2rad $ lat l dec = deg2rad $ sunDec n h = deg2rad $ sunLHA l n sinAz = - sin h * cos dec cosAz = (sin dec * cos lat') - (cos h * cos dec * sin lat') az = rad2deg $ atan2 sinAz cosAz |

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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 |
-- Sunrise, sunset diffThreshold :: (Fractional a) => a diffThreshold = 0.2 timeDiff :: (Num a) => Location -> a timeDiff l = fromIntegral (floor (lon l * 240)) findNoon :: Location -> Day -> UTCTime findNoon l d = findNoonBin l time1 time2 where time1 = addUTCTime (41400 - timeDiff l) (UTCTime d 0) time2 = addUTCTime (45000 - timeDiff l) (UTCTime d 0) findNoonBin l t1 t2 | diffUTCTime t2 t1 < diffThreshold = t1 | elev1 < elev2 = findNoonBin l timeMid t2 | otherwise = findNoonBin l t1 timeMid where elev1 = sunElev l (j2000 t1) elev2 = sunElev l (j2000 t2) timeMid = addUTCTime (diffUTCTime t2 t1 / 2) t1 elevMid = sunElev l (j2000 timeMid) sunBinSearch :: Location -> UTCTime -> UTCTime -> Maybe UTCTime sunBinSearch l time1 time2 | elev1 * elev2 > 0 = Nothing | diffUTCTime time2 time1 < diffThreshold = Just time1 | elev1 * elevMid <= 0 = sunBinSearch l time1 timeMid | elevMid * elev2 <= 0 = sunBinSearch l timeMid time2 | otherwise = Nothing where elev1 = sunElev l (j2000 time1) elev2 = sunElev l (j2000 time2) timeMid = addUTCTime (diffUTCTime time2 time1 / 2) time1 elevMid = sunElev l (j2000 timeMid) sunrise :: Location -> Day -> Maybe UTCTime sunrise l d = sunBinSearch l midnight noon where noon = findNoon l d midnight = addUTCTime (-43200) noon sunset :: Location -> Day -> Maybe UTCTime sunset l d = sunBinSearch l noon midnight where noon = findNoon l d midnight = addUTCTime 43200 noon dayLength :: Location -> Day -> Maybe DiffTime dayLength l d = maybeDiff time1 time2 where maybeDiff Nothing _ = Nothing maybeDiff _ Nothing = Nothing maybeDiff (Just t1) (Just t2) = Just $ secondsToDiffTime $ floor $ diffUTCTime t1 t2 time1 = sunset l d time2 = sunrise l d |

## 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:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 |
$ day-length-hs -s 2015-12-10 -e 2016-01-10 Location: Lat 52.2; Lon 20.9 Time zone offset: +0100 2015-12-10: rise 07:40:48, set 15:17:00, length = 7h 36min 12s 2015-12-11: rise 07:41:53, set 15:16:50, length = 7h 34min 56s 2015-12-12: rise 07:42:56, set 15:16:44, length = 7h 33min 47s 2015-12-13: rise 07:43:56, set 15:16:41, length = 7h 32min 44s 2015-12-14: rise 07:44:54, set 15:16:41, length = 7h 31min 47s 2015-12-15: rise 07:45:48, set 15:16:46, length = 7h 30min 58s 2015-12-16: rise 07:46:39, set 15:16:54, length = 7h 30min 15s 2015-12-17: rise 07:47:27, set 15:17:05, length = 7h 29min 38s 2015-12-18: rise 07:48:11, set 15:17:21, length = 7h 29min 9s 2015-12-19: rise 07:48:53, set 15:17:39, length = 7h 28min 46s 2015-12-20: rise 07:49:31, set 15:18:02, length = 7h 28min 31s 2015-12-21: rise 07:50:05, set 15:18:28, length = 7h 28min 22s 2015-12-22: rise 07:50:37, set 15:18:57, length = 7h 28min 20s 2015-12-23: rise 07:51:04, set 15:19:30, length = 7h 28min 25s 2015-12-24: rise 07:51:29, set 15:20:07, length = 7h 28min 37s 2015-12-25: rise 07:51:50, set 15:20:47, length = 7h 28min 56s 2015-12-26: rise 07:52:07, set 15:21:30, length = 7h 29min 22s 2015-12-27: rise 07:52:21, set 15:22:16, length = 7h 29min 55s 2015-12-28: rise 07:52:31, set 15:23:06, length = 7h 30min 35s 2015-12-29: rise 07:52:37, set 15:23:59, length = 7h 31min 21s 2015-12-30: rise 07:52:41, set 15:24:55, length = 7h 32min 14s 2015-12-31: rise 07:52:40, set 15:25:55, length = 7h 33min 14s 2016-01-01: rise 07:52:36, set 15:26:57, length = 7h 34min 21s 2016-01-02: rise 07:52:28, set 15:28:02, length = 7h 35min 33s 2016-01-03: rise 07:52:17, set 15:29:10, length = 7h 36min 53s 2016-01-04: rise 07:52:03, set 15:30:21, length = 7h 38min 18s 2016-01-05: rise 07:51:44, set 15:31:35, length = 7h 39min 50s 2016-01-06: rise 07:51:23, set 15:32:51, length = 7h 41min 28s 2016-01-07: rise 07:50:58, set 15:34:10, length = 7h 43min 12s 2016-01-08: rise 07:50:29, set 15:35:31, length = 7h 45min 1s 2016-01-09: rise 07:49:57, set 15:36:54, length = 7h 46min 56s 2016-01-10: rise 07:49:22, set 15:38:20, length = 7h 48min 57s |

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