How can I compute sunrise/sunset times?

Question

I like to keep track of sunrise and sunset times. For the past couple of years I've been doing this with a small program written with a popular library for my favorite programming language. The last two months I've been keeping track of these times more regularly than usual, and I happened to notice that on the day of the equinox the sunrise time jumped eight minutes as compared to the day before! I knew this was impossible and compared with NOAA, finding out that my rise and set times had been off for several days and in fact seemed to be off by about a minute for most of the year.

At this point, I'd like to just implement the calculations myself. What algorithms or formulas are available to do this computation?

Answer 1

You may consider reading Wikipedia's article on sunrise equations. The lead paragraph gives the equation:

cos(ω_o) = -tan(φ) * tan(δ)

where:

• ω_o is the hour angle in degrees at either sunrise (when negative value is taken) or sunset (when positive value is taken) in degree (°)
• φ is the latitude of the observer on the Earth in degrees
• δ is the sun declination in degrees

Answer 2

If you want to match NOAA, you'll have to consult Jean Meeus' Astronomical Algorithms (mostly Chapter 15). And it's complicated! Martin Beckett is correct, you have to define the sunset. Typically this is the apparent rise or set of upper limb of the sun, which makes your "standard" altitide -5/6 degrees (not zero). And you can't calculate the sunrise or sunset directly with NOAA's accuracy. You'll have to create a governing set of equations for apparent right ascension and declination for the day in question and then interpolate the apparent right ascension and declination over time to find the exact rising and setting times at the standard altitude.

Hope this helps. I spent about a month digesting AA and re-writing all our solar code when I came across the same thing, and it still took over a year to sort out some of the corner cases where my code broke. So it will take some time to figure out. I'm not aware of any public code examples of this algorithm and don't have any to share at this moment, but I'm happy to help you through some headaches if I can.

Answer 3

Definitions of what constitutes sunrise/sunset can vary. For example, in ephem "rising and setting are defined as the moments when the upper limb of the body touches the horizon (that is, when the body’s alt plus radius equals zero)" [PyEphem Quick Reference].

#!/usr/bin/env python
import datetime
import ephem # to install, run `pip install pyephem`

o = ephem.Observer()
o.lat, o.long, o.date = '34:3', '-118:15', datetime.datetime.utcnow()
sun = ephem.Sun(o)
print "Los Angeles"
print "sunrise:", o.next_rising(sun), "UTC"
print "sunset:",o.next_setting(sun), "UTC"

Output

Los Angeles, CA
sunrise: 2010/3/30 13:42:43 UTC
sunset: 2010/3/30 02:11:50 UTC

If it is open-source library then you could fix it instead of creating a new one with new bugs.

Answer 4

For accuracy in the 5min range you have to consider 'which' sunset.
Do you want the time the bottom of the sun touches the horizon or the time the top of the sun passes below the horizon?

It takes 2mins for the sun to cross the horizon.
Below the 1min level you also need to take into account atmospheric refraction.

Answer 5

Have a look at this book:
"Practical Astronomy with your Calculator (Paperback)" by Peter Duffett-Smith.
It is quite old but still in print... link

Answer 6

In Ruby I wrote this for equation of time.

include Math

# degrees to radians = PI/180
to_r = PI/180.0

#radians to degrees = 180/PI
to_d = 180.0/PI

puts "Day, Declination, EofT"

# test a celestial year worth of values.
for jday in 1..366

  et = -7.633 * sin(jday * (2 * PI)/365.24) + 9.65 * sin((jday - 78) * 180/92 * to_r)
  a_sin = sin(23.433 * to_r) * sin((2 * PI/366) * (jday - 81))
  declination = asin(a_sin) * to_d
  puts "#{jday}, #{declination}, #{et}"

end

Then for the above equation:

# center disk and refraction factor have been considered.
cos_omega = sin(-0.83 * to_r) - tan(latitude * to_r) * tan(declination * to_r)
semi_diurnal_arc = acos(cos_omega)

There is a whole website devoted to this at http://www.analemma.com/

• The key to these calculations are good libraries like the one for python above.
• Also usage of Date and Time classes. I would look on rubyforge for something like ephem.

Answer 7

I did see some interesting things on this NOAA page under the Technical Definitions and Computational Details heading, but I'm sure you've read that already.

The answer to the SO question "Position of the sun given time of day, and lat/long" and the above may actually be all you need.

As a side note (it doesn't answer your question directly), is there a reason you can't pull the NOAA data and use it as a lookup table instead of calculating it? Storage tends to be relatively cheap these days.

Answer 8

You might want to look at an earlier entry on calculating the position of the sun. Specifically, the Solpos program that I pointed to has support for sunrise/sunset.