[ create a new paste ] login | about

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

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


Output:
1
2
-1
114


Create a new paste based on this one


Comments: