[ create a new paste ] login | about

Project: programmingpraxis
Link: http://programmingpraxis.codepad.org/lb9mRWuW    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Sep 20:
; 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))))))


Output:
1
2
3
4
0
11
53
0


Create a new paste based on this one


Comments: