Weather is among one of my hobbies and as such, I love to watch windy and various satellite images (Such as en.sat24.com ) as well as enjoying the rain and storms outside and eventually taking images. Unfortunately I don’t have the time to focus entirely on weather forecasts and mathematics models, as it would consume all of my free time – So far,I can only say, that there are some awesome books about this topic – most notably the Atmospheric Science. Some of those standard modern forecasts however also lets you see the sunset-sunrise timing during each day of the year along with a note of whether the day is going to be longer or shorter than the previous day. We are all quite used to it, aren’t we? I have realized that over the year, this “Day contraction and expansion” actually changes and varies between +- 5 minutes per day.
I have decided to verify that myself by simply simulating the movement of the Earth around the sun along with the rotation of the Earth. As strange as it may sound, it is however not that hard at least for a very good approximation with several simplifications such as that the earth trajectory is not an ellipse (The eccentricity is 0.0167, but since its almost circular I have omitted elliptical trajectory), the Earth is a perfect sphere with a radius of 6378 km (Not an ellipsoid or some advanced GNSS Earth model) and rotates 360° each 24 hours, which is also not quite true.
The next simplification I have done is to shrink the sun within a single discrete point in space. If you love watching sunsets as I do, you already know that it takes some time for the sun to come up over 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 more commonly known as discretization – in our case discretization in time. Obviously, if your intent is to obtain precision in minutes, the discretization should be smaller. A 1-second resolution leads to very good results, but requires quite heavy computation: evaluation of whether the sun is visible 60x60x24x365 times (~ 31 M. 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-plane = 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-coordinate 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 (Mostly involving rotation in 3D space). For each of our simulation points, we rotate 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 shining.
When this single information is evaluated in all of the discretization 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 Earth Coordinates chosen. At this point, I would like to just point out, that in fact, the longitude coordinate is irrelevant and 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 have had some problems during computing though, which I would like to mention. At first, I was using the Earth’s semi major Earth Orbital Axis to be 149.60×109 meters. This sounds reasonable (As its the real distance), but not as long as we are computing in single precision format and calculating the sqrt(x2 + y2 + z2). 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 the computation time is ~ 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 GTX1060. The simulation time went down from approximately 300 seconds to just under 1 second, which is a reasonable time. Both the CUDA codes and Matlab Simulation files should be freely downloadable below. The CUDA Project requires CUDA10.2 and is intended for VS2019. 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 – custom daylight detection algorithm.
The Project is freely available ☀️ ➡️ HERE ⬅️