Solar Compass
February 7, 2012
We begin by specifying trigonometric functions that work in degrees, rather than the radians used by Scheme; it’s easier that way than to constantly convert back and forth:
(define pi (* 4 (atan 1)))
(define (degrees->radians d) (* d pi (/ 180)))
(define (radians->degrees r) (* r 180 (/ pi)))
(define (sine d) (sin (degrees->radians d)))
(define (cosine d) (cos (degrees->radians d)))
(define (arc-sine x) (radians->degrees (asin x)))
(define (arc-cosine x) (radians->degrees (acos x)))
Then the calculations are tedious but not tricky:
(define (declination yr mon day)
(let ((d (- (julian yr mon day)
(julian yr 1 1) -1))
(len (- (julian yr 12 31)
(julian yr 1 1) -1)))
(* 23.45 (sine (* 360 (/ len) (- d 81))))))
(define (equation-of-time yr mon day)
(let* ((d (- (julian yr mon day)
(julian yr 1 1) -1))
(len (- (julian yr 12 31)
(julian yr 1 1) -1))
(b (* 360 (/ len) (- d 81))))
(- (* 9.87 (sine (+ b b)))
(* 7.53 (cosine b))
(* 1.5 (sine b)))))
(define (time-correction yr mon day tz dst? long)
(+ (* 4 (- long (* 15 (+ tz (if dst? 1 0)))))
(equation-of-time yr mon day)))
(define (local-solar-time yr mon day hr min tz dst? long)
(+ hr (/ min 60) (if dst? 1 0)
(/ (time-correction yr mon day tz dst? long) 60)))
(define (hr-angle yr mon day hr min tz dst? long)
(* 15 (- (local-solar-time yr mon day hr min
tz dst? long) 12)))
(define (elevation yr mon day hr min tz dst? lat long)
(let* ((d (declination yr mon day))
(hra (hr-angle yr mon day hr min tz dst? long)))
(arc-sine (+ (* (sine d) (sine lat))
(* (cosine d) (cosine lat) (cosine hra))))))
(define (azimuth yr mon day hr min tz dst? lat long)
(let* ((d (declination yr mon day))
(hra (hr-angle yr mon day hr min tz dst? long))
(alpha (+ 90 (- lat) d))
(azi (arc-cosine (/ (- (* (sine d) (cosine lat))
(* (cosine d) (sine lat) (cosine hra)))
(cosine alpha)))))
(if (< hra 0) azi (- 360 azi))))
(define (solar-compass yr mon day hr min tz dst? lat long)
(define (approx x) (inexact->exact (round x)))
(values (approx (elevation yr mon day hr min tz dst? lat long))
(approx (azimuth yr mon day hr min tz dst? lat long))))
We round the results because of the inherent inaccuracies in the calculation. Here is the function in action:
> (solar-compass 2012 2 7 7 0 -6 #f 38.6 -90.2)
-1
114
We used julian from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/1LPgblwn, or check your calculations at http://www.srrb.noaa.gov/highlights/sunrise/azel.html. The web pages at http://www.pveducation.org/pvcdrom/properties-of-sunlight/azimuth-angle were the primary resource for this exercise.
Does anybody know how to develop Android applications?
Programming Praxis 
There is a low-tech way to build a solar compass. Create a level table. Put a pointer above it so it casts a shadow on the table. Mark the shadow’s position. Wait until the shadow has moved, then mark its new position. The ray from the first to second position points due east.
This has some error, as the sun isn’t moving* _exactly_ east-to-west. The error increases
with latitude and at times closest to the equinoxes. Here at 45° north latitude, on the equinox, the error is about 0.2°.
*Pedantically, the sun isn’t moving, the earth is rotating.
Burt’s compass does much the same thing. The latitude scale builds a sundial. Then the third scale marks the solar hour angle, from which you can calculate any compass point. And you don’t have to wait an hour.