; solar compass
(define (julian year month day)
(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))))
; trig functions in degrees
(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)))
; yr, mon, day, hr, min -- normal reckoning of time (24-hour clock)
; tz -- hours east (positive) or west (negative) of greenwich
; dst? -- #t if daylight saving time, otherwise #f
; lat -- degrees north (positive) or south (negative) of equator
; long -- degrees east (positive) or west (negative) of greenwich
(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))))
(call-with-values
(lambda () (solar-compass 2012 2 7 7 0 -6 #f 38.6 -90.2))
(lambda (elev azi) (display elev) (newline) (display azi) (newline)))