Programming Praxis

Sunrise and sunset are processes, rather than events, and picking one moment as the exact moment of sunrise or sunset is difficult. Astronomical twilight occurs when the sun begins to obscure faint stars and other celestial objects. Nautical twilight occurs when the horizon becomes visible. Civil twilight occurs when the sky is illuminated. Sunrise and sunset occur when the highest portion of the sun’s disk is just barely visible above the horizon. “True” sunrise and sunset occur when the center of the solar disk touches the geometric horizon. It is odd but true that sunrise and sunset occur when the sun is still technically below the geometric horizon because the earth’s atmosphere causes the light to bend, as the atmosphere refracts light differently than the vacuum of space. Our calculations below, which are based on the Wikipedia article “Sunrise equation,” calculate the true geometric time of sunrise and sunset, and thus may differ by a minute or two from the sunrise and sunset times given in the local newspaper, which typically account for refraction; for simplicity, we ignore those times and places where the sun will not rise or set during the day.

The simple formula for sunrise and sunset is cos(ω) = −tan(φ) × tan(δ), where ω is the hour angle of the sun at sunrise (the negative cosine) or sunset (the positive cosine), φ is the latitude of the observer on the earth, and δ is the declination of the sun, all angles in degrees. Since the earth rotates at 15° per hour, ω/15 gives the time of sunrise as the number of hours before the local noon (the time when the sun transits the observer’s meridian) and the time of sunset as the number of hours after the local noon.

We will write a single function solar that takes six arguments and returns two values, the times of sunrise and sunset. The six arguments are the latitude (north positive, south negative) and longitude (west positive, east negative) of the observer, the four-digit year, month and day of the observation, and the number of hours (west positive, east negative) of the local time zone from Greenwich mean time. In outline, the function looks like this:

(define (solar lat-north long-west year month day gmt-offset)
  ... local definitions ...
  (let ( ... preparatory calculations ...)
    (values (time rise) (time set))))


We will follow Wikipedia‘s calculations, giving only the briefest of explanations of the astronomy involved, since astronomy is a very large topic and this is a very small exercise. D is the julian date, the number of days since January 1, 4713 BC; this is the same calculation we made in two previous exercises:

(d (let* ((a (quotient (- 14 month) 12))
          (y (+ year 4800 (- a)))
          (m (+ month (* 12 a) -3)))
     (+ day
        (quotient (+ (* 153 m) 2) 5)
        (* 365 y)
        (quotient y 4)
        (- (quotient y 100))
        (quotient y 400)
        (- 32045))))

N is the julian cycle, which adjusts from the Greenwich noon of the julian date to local noon at the observer’s location:

(n (round (- d 2451545.0009 (/ long-west 360))))

J is the approximate solar noon at the observer’s location; we’ll refine this calculation momentarily:

(j (+ 2451545.0009 (/ long-west 360) n))

M is the solar mean anomaly, which tells how far along the earth is on its annual orbit around the sun:

(m (mod360 (+ 357.5291 (* 0.98560028 (- j 2451545)))))

C is the equation of the center, which adjusts for the difference between a perfect circle for the orbit of the earth around the sun and the actual elliptical path of the orbit:

(c (+ (* 1.9148 (sin-d m))
      (* 0.0200 (sin-d (* 2 m)))
      (* 0.0003 (sin-d (* 3 m)))))

E is the ecliptic longitude, which specifies the position of the earth on the celestial sphere, alternately specifies the path the sun appears to follow across the sky over the course of a year:

(e (mod360 (+ m 102.9372 c 180)))

T is the solar transit, which is the exact time of solar noon when the sun transits the celestial meridian:

(t (+ j (* 0.0053 (sin-d m)) (- (* 0.0069 (sin-d (* 2 e))))))

S is the declination the sun, alternately the tilt of the earth’s axis:

(d (asin-d (* (sin-d e) (sin-d 23.45))))

H is the hour angle, the main equation given above that calculates the time of sunrise and sunset, adjusted for the apparent size of the solar disc, approximately 50 arc-minutes:

(h (acos-d (/ (- (sin-d -0.83) (* (sin-d lat-north) (sin-d d)))
              (* (cos-d lat-north) (cos-d d)))))

Finally, we calculate the exact julian times of sunrise and sunset, in Greenwich mean time, at the observer’s location:

(set (+ 2451545.0009
        (+ (/ (+ h long-west) 360)
           n
           (* 0.0053 (sin-d m))
           (- (* 0.0069 (sin-d (* 2 e)))))))
(rise (- t (- set t))))

The conversion to local time is given below.

The formulas we have used are all stated in degrees, but Scheme provides all trigonometric functions in radians. The following conversions are part of the local definitions:

(define pi 3.141592654)

(define (d->r d) (* d pi 1/180))
(define (r->d r) (* r 180 (/ pi)))

(define (sin-d x) (sin (d->r x)))
(define (cos-d x) (cos (d->r x)))
(define (tan-d x) (tan (d->r x)))
(define (asin-d x) (r->d (asin x)))
(define (acos-d x) (r->d (acos x)))
(define (atan-d x) (r->d (atan x)))

Another local definition forces an angle to be in the range 0° ≤ θ < 360°:

(define (mod360 x)
  (if (< x 0) (+ x 360)
    (if (<= 360 x) (- x 360) x)))

The final local definition converts a fractional day into an hour-and-minute time, considering the local gmt-offset:

(define (time j)
  (let* ((j1 (+ (- j (floor j)) (/ gmt-offset 24)))
         (j2 (round (* j1 24 60)))
         (j3 (floor (/ j2 60)))
         (j4 (- j2 (* j3 60))))
    (string-append
      (number->string (inexact->exact j3)) ":"
      (number->string (inexact->exact j4)))))

The entire function is assembled at http://programmingpraxis.codepad.org/6EHh1miI. Here’s an example:

> (solar  38.623944 90.187235  2009 11 24  6)
6:54
16:44

It is a fascinating exercise to calculate the various twilight times. If you are interested, you may wish to consult the Sunrise/Sunset Algorithm from the Almanac for Computers from the Nautical Almanac Office of the United States Naval Observatory, available at http://williams.best.vwh.net/sunrise_sunset_algorithm.htm.

Pages: 1 2