; autumn equinox
(define (string-split sep str)
(define (f cs xs) (cons (list->string (reverse cs)) xs))
(let loop ((ss (string->list str)) (cs '()) (xs '()))
(cond ((null? ss) (reverse (if (null? cs) xs (f cs xs))))
((char=? (car ss) sep) (loop (cdr ss) '() (f cs xs)))
(else (loop (cdr ss) (cons (car ss) cs) xs)))))
(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))))
(define (julian year month day hour minute second)
(let* ((a (quotient (- 14 month) 12))
(y (+ year 4800 (- a)))
(m (+ month (* 12 a) -3)))
(+ day (/ hour 24) (/ minute 24 60) (/ second 24 60 60)
(quotient (+ (* 153 m) 2) 5)
(* 365 y)
(quotient y 4)
(- (quotient y 100))
(quotient y 400)
(- 32045))))
(define (gregorian julian)
(let* ((j (floor julian))
(hms (* (- julian j) 24 60 60))
(hour (floor (/ hms 60 60)))
(hms (- hms (* hour 60 60)))
(minute (floor (/ hms 60)))
(second (- hms (* minute 60)))
(j (+ j 32044))
(g (quotient j 146097))
(dg (modulo j 146097))
(c (quotient (* (+ (quotient dg 36524) 1) 3) 4))
(dc (- dg (* c 36524)))
(b (quotient dc 1461))
(db (modulo dc 1461))
(a (quotient (* (+ (quotient db 365) 1) 3) 4))
(da (- db (* a 365)))
(y (+ (* g 400) (* c 100) (* b 4) a))
(m (- (quotient (+ (* da 5) 308) 153) 2))
(d (+ da (- (quotient (* (+ m 4) 153) 5)) 122))
(year (+ y (- 4800) (quotient (+ m 2) 12)))
(month (+ (modulo (+ m 2) 12) 1))
(day (+ d 1)))
(values year month day hour minute second)))
(define (duration j-start j-finish)
(if (< j-finish j-start) (duration j-finish j-start)
(let* ((diff (- j-finish j-start))
(days (floor diff))
(hms (* (- diff days) 24 60 60))
(hours (floor (/ hms 60 60)))
(hms (- hms (* hours 60 60)))
(minutes (floor (/ hms 60)))
(seconds (- hms (* minutes 60))))
(values days hours minutes seconds))))
(call-with-values
(lambda () (solar 0 0 2012 9 21 0))
(lambda (set rise)
(let* ((start-hour (string->number (car (string-split #\: rise))))
(finish-hour (string->number (car (string-split #\: set))))
(start-minute (string->number (cadr (string-split #\: rise))))
(finish-minute (string->number (cadr (string-split #\: set)))))
(call-with-values
(lambda () (duration (julian 2012 9 21 start-hour start-minute 0)
(julian 2012 9 21 finish-hour finish-minute 0)))
(lambda (d h m s)
(display d) (newline)
(display h) (newline)
(display m) (newline)
(display s) (newline))))))