Every smartphone nowadays has a weather forecast and some more advanced forecasts lets you see the sunrise and sunset time along with the information whether the day is shorter or longer than the day before and also the difference of the time. I was wondering for a while now how this is calculated and whether its even correct. Of course that from summer to winter, the days are getting shorter and longer from winter to summer. That is likely something no one will argue about, but how is it that sometime it says +5 minutes,sometimes +4 minutes and sometime there is no difference at all?

Imagine a real scenario. The Earth rotating 360° each 24 hours (This is one of the simplifications as days are not exactly 3600*24 seconds) at its axial tilt of approximately 23.48° Deg (At least for Earth  ). If there was no axis tilt and you were standing on the Equator, you would definitely measure a total daylight time of 12 Hours and nighttime of additional 12 hours. To be honest,if there will be no axis tilt and the Earth would move through its trajectory around the sun during the year, you would still find the same results given that the trajectory is a circle and not an ellipse. But thanks to “Wikipedia” ,we know that all the planets do move on an elliptical trajectory, but most of them are nearly circular. The Earth’s eccentricity is approximately 0.0167,which leads to additional simplification, that the trajectory is just a circle. Whats not as well visible on the first sight is that the earth is not a sphere, but an elipsoid. There are some models especially for GNSS systems,that do describe the earth more precisely, but we will as well simplify our model by stating that the Earth is a sphere with a radius of 6378 km.

The next simplification we have to make is to discuss the amount of daylight during sunset/sunrise,because we know that the Sun is not really a single point in space,but is in fact quite huge. As a result, the sun slowly appears/disappears behind the horizon. My best guess  of the time required is approximately 1 minute,but this will vary based on the longitude as sunsets/sunrises are fastest on the equator. There is still something worth considering however. That is the precision of the simulation,which shrinks to one “two” considerations. We know,that our model will not be extra precise,but estimated errors should not be high and some of them will cancel each other or just shift up/down the final chart. The talk is about the simulation time at first,It may not be completely necessary to evaluate the model for each second during the year,because we know,that model is not exact on a 1 second scale. I personally would say a resolution of 1 minute is well sufficient,which still means evaluation the earth movement and rotation around the sun 60*24*365 times. 

Now if we define our location it terms latitude/longitude and change coordinate system to cartesian xyz, we may define the Earth’s trajectory to be in z=0,which is not true for the solar system,but since the Sun is a Symmetric sphere,we can ignore inclination angles and whatever and just say that rotation plane is just where z is zero. Next we can define the sun to be the center of the universe, centered at coordinates [0,0,0]. What follows is basically some applied math from the high school and linear algebra from the university. For each of our simulation points,we rotate the our position by angle appropriate for the daytime and by Earth’s axis tilt. When this is done, compute the intersection of a line defined with 2 points (Our position and Position of the sun) with a Sphere (The Earth). Since our location is on the earth’s surface the solution will likely provide 2 results of intersection,where one will be likely our position and the other arbitrary. What needs to be decided is whether  our position is closer to the sun that the other intersection point. If so,we can say that the sun is visible and thus evaluate that for the given point in time,the sun is visible.

When this single information is evaluated from all 60*24*365 points,we can apply some custom algorithm to the data to extract amount of time when the sun is visible and finally plot the Daylight chart for x varying from 1 to 365 and showing the total amount of daylight for the position chosen. At this point,I would like to just point out,that in fact,the longitude coordinate irrelevant may be chosen arbitrary,since its just a rotation in time and represents a small shift (up to 1/365) to the final chart.

As you can see, the results match the expected daylight difference (-5,5) minutes even though we have made some simplifications. I had some problems at first, which i would like to mention. At first, I was using the Earth’s semi Major Earth Orbital Axis to be 149.60e9 meters. This sounds reasonable, but not as long as we are talking in single precision format and calculating the sqrt(x^2 + y^2 + z^2) here, we are already at the half of the floating point precision limits and the usage of double precision format becomes mandatory. At second, I have found that even using double precision format doesn’t necessary provide error-free results and the whole model needs to be scaled down a bit. What is more is that I was not using the recommended resolution of 1 minute discretization in time, but 60x higher resolution. This leads to computation time of 5 minutes on the latest i7-8700K CPU.

Since I don’t like long computation times, I have recoded the simulation to run in nVidia CUDA on my GeForce GTX960. The simulation time went down from approximately 300 seconds to 1 second, which is a reasonable time. Both the CUDA codes and Matlab Simulation files should be freely downloadable. The CUDA Project requires CUDA10 and is intended for VS2017. The single/double precision computation could be coded smarter, but who cares? At last, I would like to point out, that estimating the Daylight of each day is not as easy as it sounds. For sure we know that the year has 365 days, but lets say that somebody uses location at longitude of 85°. What is going to happen is that you would likely see a polar day and a polar night. In that case, one cannot just iterate through the edges visible/not visible. In that case, its likely better to just divide the total simulation time into 365 pieces and sum them up

If you find any errors, or problems, please write below in the comments or feel free to contact me.
The code is freely downloadable and modifiable for non-commercial purposes.