[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Nov 22:
(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)))


Output:
1
2
Sunrise: 6:54
Sunset: 16:44


Create a new paste based on this one


Comments: