Ever been stuck in an office where you couldn't see outside? Or wanted to know if the sun was up somewhere else on the globe? You might find this little page useful.

While working at AFGWC on a radiation package in the '80s for a climate / weather forecast model, I wanted something that could give me the angle of the sun from a grid point on the planet. This algorithm is the result.

It's not terribly fast, but the original was extremely accurate for the 50 years' worth of astronomical data I could find. The requirement wasn't for accuracy, but speed. Oh, well. Over the decades I lost the code and any notes I'd made, and this is the recreation from memory.

The dot-product of two normalized vectors in cartesian space is the cosine of the angle of separation (see the box below). The apparent latitude of the sun doesn't change much over the course of a day so I use an approximation. I simply did the dot-product and solved for the unknown: the longitude of the sun. This is all observer-centric (no offset from the center of the planet - though I think the original used a translation of coordinates).

The two constants I calculated for the two waves for the Equation Of Time are from memory. The only part I vaguely remembered was the math.

Derivation of the formula for θ_{s}:

s⃗ * o⃗ = cos( δ )

s⃗ = ( x_{s},y_{s},z_{s} ) o⃗ = ( x_{o},y_{o},z_{o} )

x = r * cos( φ ) * cos( θ )

y = r * cos( φ ) * sin( θ )

z = r * sin( φ )

normalize: r ≡ 1

cos( δ ) = x_{s} * x_{o} +y_{s} * y_{o} +z_{s} * z_{o}

cos( δ ) = a * cos( θ_{s} ) + b * sin( θ_{s} ) + c

a = cos( φ_{o} ) * cos( θ_{o} ) * cos( φ_{s} )

b = cos( φ_{o} ) * sin( θ_{o} ) * cos( φ_{s} ) which drops if θ_{o} ≡ 0

c = sin( φ_{o} ) * sin( φ_{s} )

d = a * cos( θ_{s} ) + b * sin( θ_{s} )

using the double angle identity:

d = e * sin( α + θ_{s} )

d = e * cos( α ) * cos( θ_{s} ) - e * sin( α ) * sin( θ_{s} )

a = e * cos( α )

b = -e * sin( α )

tan( α ) = sin( α ) / cos( α )

tan( α ) = ( -b/e ) / ( a/e ) = -b/a

α = tan^{-1}( -b/a )

a^{2} = e^{2} * cos^{2}( α )

b^{2} = e^{2} * sin^{2}( α )

a^{2} + b^{2} = e^{2} * ( cos^{2}( α ) +sin^{2}( α ) )

a^{2} + b^{2} = e^{2}

e = √( a^{2} + b^{2} )

d = e * cos( α + θ_{s} )

d/e = cos( α + θ_{s} )

α = cos^{-1}( d/e )

θ_{s} = cos^{-1}( d/e )

To see what's actually implemented... best to look at the code.

s⃗ * o⃗ = cos( δ )

s⃗ = ( x

x = r * cos( φ ) * cos( θ )

y = r * cos( φ ) * sin( θ )

z = r * sin( φ )

normalize: r ≡ 1

cos( δ ) = x

cos( δ ) = a * cos( θ

a = cos( φ

b = cos( φ

c = sin( φ

d = a * cos( θ

using the double angle identity:

d = e * sin( α + θ

d = e * cos( α ) * cos( θ

a = e * cos( α )

b = -e * sin( α )

tan( α ) = sin( α ) / cos( α )

tan( α ) = ( -b/e ) / ( a/e ) = -b/a

α = tan

a

b

a

a

e = √( a

d = e * cos( α + θ

d/e = cos( α + θ

α = cos

θ

To see what's actually implemented... best to look at the code.

Equation of Time approximation (eccentricity +obliquity):

eot( date ) = effectEccentricity +effectObliquity

effectEccentricity = EOT_DEVIATION_PERIGEE *sin( 2 *π *( date -APSIS_NEAR_PERIHELION) /YEAR_LENGTH_ANOMALISTIC )

effectObliquity = EOT_DEVIATION_SOLSTICE *sin( 4 *π *( date -SOLSTICE_EVENT) /YEAR_LENGTH_TROPICAL )

EOT_DEVIATION_PERIGEE ≊ 480s (from memory)

EOT_DEVIATION_SOLSTICE ≊ 600s (from memory)

eot( date ) = effectEccentricity +effectObliquity

effectEccentricity = EOT_DEVIATION_PERIGEE *sin( 2 *π *( date -APSIS_NEAR_PERIHELION) /YEAR_LENGTH_ANOMALISTIC )

effectObliquity = EOT_DEVIATION_SOLSTICE *sin( 4 *π *( date -SOLSTICE_EVENT) /YEAR_LENGTH_TROPICAL )

EOT_DEVIATION_PERIGEE ≊ 480s (from memory)

EOT_DEVIATION_SOLSTICE ≊ 600s (from memory)

Derivation of the double angle identity:

e^{± i * α} = cos( α ) ± i * sin( α ) (by Euler)

e^{± i * β} = cos( β ) ± i * sin( β )

e^{i * α} * e^{i * β} = e^{i * ( α + β )}

e^{i * ( α + β )} = cos( α + β ) + i * sin( α + β )

∴

cos( α + β ) + i * sin( α + β ) =

( cos( α ) + i * sin( α )) * ( cos( β ) + i * sin( β ))

or

cos( α + β ) = cos( α ) * cos( β ) - sin( α ) * sin( β )

sin( α + β ) = cos( α ) * sin( β ) + sin( α ) * cos( β )

e

e

e

e

∴

cos( α + β ) + i * sin( α + β ) =

( cos( α ) + i * sin( α )) * ( cos( β ) + i * sin( β ))

or

cos( α + β ) = cos( α ) * cos( β ) - sin( α ) * sin( β )

sin( α + β ) = cos( α ) * sin( β ) + sin( α ) * cos( β )

The yellow dot is the position where the sun is directly overhead (or within a half of a degree). The blue and yellow circles are where the sun is at 30° and 60° from the vertical (take note if your skin is sensitive to harsh sunlight).

The three stripes show where it's civil, nautical, or astronomical twilight.

The ECMAScript code starts and creates however many worker threads it needs to have one worker per core on the machine. It then divides the problem up into vertical pixel lines and calculates the angle at each pixel.

The page should be updated every minute, unless the checkbox is unselected.

φ_{s} approximation:

φ_{s} ≊ 23.44 * cos(( date -SOLSTICE_EVENT ) *2π / YEAR_LENGTH_TROPICAL )

(yes, this could be greatly improved)

φ

(yes, this could be greatly improved)

latitude | ||

longitude | ||

Noon |

Event | GMT | |

Begins | Ends | |

Daylight | ||

Twilight - Civil | ||

Twilight - Nautical | ||

Twilight - Astronomical |

To get the times for sunrise, sunset, and the twilights, enter the latitude and longitude into the input fields. The code will display these values and the time of noon (when the sun is over your longitude).

All times are given in GMT. This makes it easy (or less confusing) when inspecting other latitudes and longitudes.

I haven't checked the accuracy of any of this at present, and my constants for my version of the Equation of Time are only approximations. And since I only work on Linux, and use Firefox, I suspect this might not work on other environments.