(define (solar lat-north long-west year month day gmt-offset)
(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)))
(define (mod360 x) (if (< x 0) (+ x 360) (if (<= 360 x) (- x 360) x)))
(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)))))
(let* (
; julian date
(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))))
; julian cycle
(n (round (- d 2451545.0009 (/ long-west 360))))
; approximate solar noon
(j (+ 2451545.0009 (/ long-west 360) n))
; solar mean anomaly
(m (mod360 (+ 357.5291 (* 0.98560028 (- j 2451545)))))
; equation of the center
(c (+ (* 1.9148 (sin-d m))
(* 0.0200 (sin-d (* 2 m)))
(* 0.0003 (sin-d (* 3 m)))))
; ecliptic longitude
(e (mod360 (+ m 102.9372 c 180)))
; solar transit
(t (+ j (* 0.0053 (sin-d m)) (- (* 0.0069 (sin-d (* 2 e))))))
; declination of the sun
(s (asin-d (* (sin-d e) (sin-d 23.45))))
; hour angle
(h (acos-d (/ (- (sin-d -0.83) (* (sin-d lat-north) (sin-d s)))
(* (cos-d lat-north) (cos-d s)))))
; sunrise and sunset
(set (+ 2451545.0009
(+ (/ (+ h long-west) 360)
n
(* 0.0053 (sin-d m))
(- (* 0.0069 (sin-d (* 2 e)))))))
(rise (- t (- set t))))
; return sunrise and sunset
(values (time rise) (time set))))
(call-with-values
(lambda () (solar 38.623944 90.187235 2009 11 24 6))
(lambda (rise set)
(display "Sunrise: ") (display rise) (newline)
(display "Sunset: ") (display set) (newline)))