; phases of the moon

(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))))

(define (days year month day)
  (define (normalize x)
    (let ((x (- x (floor x))))
      (if (negative? x) (+ x 1) x)))
  (let ((new 2451550.1) (moon 29.530588853)
        (j (julian year month day)))
    (* (normalize (/ (- j new) moon)) moon)))

(define (phase year month day)
  (let ((d (days year month day)))
    (cond ((< d  1.84566) "New")
          ((< d  5.53699) "Waxing crescent")
          ((< d  9.22831) "First quarter")
          ((< d 12.91963) "Waxing gibbous")
          ((< d 16.61096) "Full")
          ((< d 20.30228) "Waning gibbous")
          ((< d 23.99361) "Last quarter")
          ((< d 27.68493) "Waning crescent")
          (else           "New"))))

(display (phase 2000 1 6)) (newline)
(display (days 2010 1 22)) (newline)
(display (phase 2010 1 22)) (newline)

New
7.106982227906547
First quarter