Math behind a world sunlight map


World sunlight map fractionMy neighbour has a map of the world on the wall. You can see it from the street in front of his house. It has a backlight but that only illuminates half of the map. The transition from day to night is shaped like a sine wave most of the time. It actually is a physical world sunlight map. Of course, you can simulate this with a computer too. There even is an instance using Google maps.

As many roads lead to Rome multiple ways are possible to this simulation. One could model the sun, earth, maybe more and start ray tracing. This approach would include solar eclipses but is quite heavy by means of the load on the processor. Because of the number of calculations involved in ray tracing is quite high. The way I choose to describe fully in this article is one close to it. Using vectors pointing from a sphere (earth) to a point (sun) I map a Mercator projected map of the world on the sphere. The challenges included are the yearly orbit of earth around the sun and it’s 23.5° tilted 24 hour spin.

About a decade ago I made a similar program in highschool. I didn’t know vectors and trigonometry then as I do now. Back then I used to have the map of the world where for each longitude (the ones parallel to the equator) I calculated sunset and sunrise times for that longitude. This way I knew when to start painting night over the daylight map I had. This result was pretty accurate and an algorithm to determine those solar times is not hard to implement. This time I used another more advanced approach to the problem.

Space

First let us define a mathematical 3D space, with three axis: x (pointing from left to right), y (pointing up) and z (pointing into your monitor). Pointing directions are given to help your imagination. Axis point from small (negative infinite) to big (positive infinite), they intersect at O_{x,y,z} = (0, 0, 0).

Earth

Define a sphere on point O, with a radius of 1. This will be the earth. The funny thing is that for every point\vec{p}(x,y,z coordinate) on the earth, because the earth is centered around O, automatically is the same a the vector pointing orthoganally away from the earth starting on that point\vec{p}. This comes in handy when we do calculations later on.

\vec{p} = (x, y, z)

Sun

We don’t need the sun for this method. We only need the direction to it from the earth. So define a vector S which points in the direction of the sun from O. To do this we must use sine and cosine. On time t = 0 the direction of the sun will be (1, 0, 0), which is a vector pointing to the right. On t = 12 (hour) that should be pointing right (-1, 0, 0). The same goes for 6 o’clock am and pm, being respectively (0, 1, 0) and (0, -1, 0). This gives us insight that the direction \vec{d} to the sun can be given as (t should be scaled to a value from 0 to 1, so we need to divide t by 24):

\vec{d} =(\cos{\frac{2\pi t}{24}}, 0, \sin{\frac{2\pi t}{24}})

With this approach we actually make the sun rotate around the earth.

Tilted earth

Earth tilted axis

Tilted orbit of the spinning earth around the sun. The method discussed here turns it around as it models the sun orbitting the earth.

The above formula still is not complete. The axis around which the earth is spinning is not orthogonal to the plane defined by the orbit of the earth around the sun. It is tilted 23.5°. This is how we get the four seasons on earth of course. To model this a nice trick can be applied to d. We know that on or closely around June 21th the sun is right above the 23.5° latitude (Tropic of Cancer), the summer solstice. From here it takes a sine with a period of a year to go to -23.5° and back again. June 21th is the 172th day of the year. The standard cosine function has the form:

y = a \cdot \cos{(bx + c)} + d

We can simply fill this out. The a must be 23.5, b is \frac{2\pi}{365}, c is -172 and d equals 0. If you define y the observable tilt (because the real tilt will always be 23.5) and x is the day of the year:

observableTilt = 23.5 \cos{(\frac{2\pi}{365} ({dayOfYear} - 172))}

Now this observable tilt is the angle there is between what is the real direction to the sun and the already calculated direction d. To calculate the vector (y component, pointing up or down) which should be added to \vec{d} to get the real direction you can use tangens.

\vec{tiltCorrection} = \vec{d} \cdot \tan{ 2\pi \frac{observableTilt}{360}}

Then the real direction vector \vec{rd} is the addition of \vec{tiltCorrection} to \vec{d}:

\vec{rd} = \vec{tiltCorrection} + \vec{d}

Calcultating illumination and mapping

Using the dot product (or scalar product or inner product) you can calculate the angle between two vectors. Or, if the vectors are both normalized (length = 1) it is the projection of the one vector onto the other. So if we calculate the dot product for every normal vector from the surface of the earth (\vec{p}) with \vec{rd} we should get some meaningful results.

illumination = \vec{rd} \cdot \vec{p}

If the illumination < 0 the sun is in another direction so \vec{p} is pointing to the night side of the earth. If illumination \geq 0 that piece of earth recieves solar light.

We have a 2D map (Mercator projection) we need to map on the 3D sphere. We need to do this to determine which we need to calculate dot products for. To avoid confusion, lets call x and y from the 2D map u and v. We can iterate for each u and each v for each u and map it to a \vec{p} in 3D. To do this we need to determine phi (\phi, logitude: v)and theta (\theta, latitude: u). In the iteration (maxV is the height of the map and maxU is the width of the map):

\phi = 2\pi (\frac{v}{maxV})

\theta = 2\pi (\frac{u}{maxU})

Now u and v map to \vec{p} as follows:

\left\{\begin{matrix} \vec{p}_x = \sin{\phi} \cos{\theta} \\ \vec{p}_y = \cos{\theta} \\ \vec{p}_z = \sin{\phi} sin{\theta} \end{matrix}\right.

Now we can for each coordinate (u, v) on the map if it should be displayed as day or night.

Fun with shading

Because you know the angle the sunlight makes with the earth’s surface you can make add some shading by making it increasingly dark on the edges of the day. It the dot product of \vec{p} and \vec{rd} is smaller than 0.1 (illumination < 0.1) the point is in dusk or dawn and you could mix day and night as a gradient to make the difference between them look more fluently. Another joke is the reflection of the sun: if the product is greater than, let’s say, 0.95 (illumination > 0.95) the sun is approximately orthogonally above that coordinate an you could make it more white to have it look like the sun is reflecting in the map.

Demo and discussion

A demontration of this method is available at the edesign example site. There are a number of assumptions done and restrictions set while creating this demo.

Live sunlight map generated using the algorithm discussed in this article. Click the image to go to the example page where you can override current time.

  • I assumed the solar rays are parallel. The distance from the sun to the earth is huge, but it is actually wrong to assume rays are parallel. They are only by approximation so it is good enough for this simulation.
  • Because of this assumption, the night ’starts’ exactly at the half of the world in the shade while actually this is different. The sun is far bigger than the earth and therefor some rays should be able to reach the shadow half of the earth (near the limit depicted in the first image). By approximation I eliminated this too.
  • Timing is not very accurate. As well as for the 24h clock as for the solstices rough estimates are used.
  • A perfect sphere is taken as world. No fattening is modelled (polar radius should be smaller than equatorial radius).
  • The ‘fun with shading’ should implement something to make land mass not reflect sunlight.
VN:F [1.3.1_645]
Rating: 9.5/10 (14 votes cast)
  • Share/Save/Bookmark
  1. #1 by Uira Cocha - June 5th, 2009 at 22:18

    Beautiful!

    Thanks for your article and for making public the source code of its implementation in PHP

  2. #2 by WorldSunRise - November 9th, 2009 at 22:42

    Tried your script, but it doesn’t generate map picture. Not sure, but it might be because of “header already set” problem. t.i. do you have any tip? What web server you use?

  3. #3 by Admin - January 23rd, 2010 at 12:47

    These “headers already set” error is usually caused by some echo/print output or any content outside of your tags. You should locate the line causing the output of your script. Btw, I use Apache 2 httpd.

(will not be published)
  1. No trackbacks yet.