#!/local/bin/siod -v0

;; The Phase of the Moon.
;;
;; This program is in the public domain.  It is derived largely from
;; calendar algorithms published in _Numerical Recipes in Pascal_,
;; and was written for SIOD Scheme <ftp://ftp.std.com/pub/gjc/siod.tgz>.
;;
;; The file "julian.scm" converts dates between civil and Julian formats,
;; and can be found at <http://twp.rootsweb.com/~twp/julian.scm>.
;;
;; Tim Pierce <twp@rootsweb.com>

(load "julian.scm" nil t)

;; (define PI 3.14159265358979323846)
(define PI *pi*)
(define (dsin ang) (sin (/ (* ang PI) 180)))
(define epoch 2415020)  ; 1899 Dec 31

(define (flmoon n nph)  ; return (jd . frac)
  ;; From _Numerical Recipes_.  `flmoon' returns the date (returned as
  ;; a Julian Day, a real value) on which a given phase of the moon
  ;; will fall.  NPH indicates the desired moon phase (0 = new, 1 =
  ;; first quarter, 2 = full, 3 = last quarter) and N the total number
  ;; of moon cycles since the epoch (January 1, 1900).
  (if (or (> nph 3) (< nph 0))
      (error "illegal nph")
      (let* ((c (+ n (/ nph 4)))
             (t (/ c 1236.85))
             (t2 (* t t))
             (as (+ 359.2242 (* c 29.105356)))
             (am (+ 306.0253
                    (* c 385.816918)
                    (* t2 0.010730)))
             (jd (+ epoch  ; 1899 Dec 31
                    (* 28 n)
                    (* 7 nph)))
             (xtra (+ 0.75933
                      (* c 1.53058868)
                      (* t2
                         (- 1.178e-4 (* 1.55e-7 t)))
                      (if (or (= nph 0) (= nph 2))
                          (+ (* (dsin as) (- 0.1734 (* t 3.93e-4)))
                             (* (dsin am) -0.4068))
                          (+ (* (dsin as) (- 0.1721 (* t 4.0e-4)))
                             (* (dsin am) -0.6280))))))
        (let ((i (if (> xtra 0)
                     (trunc xtra)
                     (trunc (- xtra 1)))))
          (+ jd i (- xtra i))))))

(define (latest-phase jd)
  ;; Find the most recent moon phase to Julian Day JD.
  ;; Return a total number of moon *phases*; the number
  ;; of full moon cycles may be obtained from this by dividing by 4
  ;; and taking the quotient.
  ;; Start by guesstimating the number of moon phases since the epoch.
  (let* ((elapsed-days  (- jd epoch))
         (elapsed-years (/ elapsed-days 365.2422))
         (elapsed-cycles (* elapsed-years 12.3685))
         (full-cycles (floor elapsed-cycles))
         (total-phases (* full-cycles 4)))

    ;; Hunt down to find a phase earlier than the target date.
    (let ((try-prev-phase nil)
	  (first-phase nil))
      (set! try-prev-phase
	    (lambda (nph)
	      (if (< (flmoon (quotient nph 4) (remainder nph 4)) jd)
		  nph
		  (try-prev-phase (- nph 1)))))
	 (set! first-phase (try-prev-phase total-phases))
	 
	 ;; Now hunt up to find the *closest* early phase.  The result
	 ;; is returned to the calling function.
	 (let ((try-next-phase nil))
	   (set! try-next-phase
		 (lambda (nph)
		   (if (>= (flmoon (quotient nph 4) (remainder nph 4)) jd)
		       (- nph 1)
		       (try-next-phase (+ nph 1)))))
	      (try-next-phase first-phase)))))

(define (phase-of-moon jd . verbosity)
  ;; Count the number of moon phases from the epoch.
  (let* ((latest (latest-phase jd))
         (full-cycles (quotient latest 4))
         (last-phase (remainder latest 4))
         (verbose? (and (not (null? verbosity))
                        (car verbosity))))

    ;; if debugging, print out dates for all moon cycles
;;    (if verbose?
;;        (for-each (lambda (n)
;;                    (prin1 (phase-name n))
;;                    (prin1 ":\t")
;;                    (prin1 (unix-ctime (julian->UTC (flmoon full-cycles n)))))
;;                  '(0 1 2 3)))

    (puts "The Moon is ")
    (let* ((next-full-cycles (quotient (+ latest 1) 4))
           (next-phase (remainder (+ latest 1) 4))
           (time-from-last-phase (- jd (flmoon full-cycles last-phase)))
           (time-to-next-phase (- (flmoon next-full-cycles next-phase) jd))
           (time-since-new (- jd (flmoon full-cycles 0)))
           ;(time-since-full (- jd (flmoon full-cycles 2)))
           (days-old (floor (+ .5 time-since-new))))
           ;(days-old (floor (+ .5 (if (<= next-phase 2)
	   ;			      time-since-new
	   ;			      time-since-full)))))
      (puts
       (cond ((< time-from-last-phase .5) (phase-name last-phase))
             ((< time-to-next-phase   .5) (phase-name next-phase))
             ('else
              (string-append (cond
                               ((= last-phase 0)  "Waxing Crescent, ")
                               ((= last-phase 1)  "Waxing Gibbous, ")
                               ((= last-phase 2)  "Waning Gibbous, ")
                               ('else "Waning Crescent, "))
                             (number->word days-old)
                             (if (> days-old 1)
                               " days"
                               " day")
                             " old"))))
      (puts ".\n"))))

(define (number->word n)
  (let ((list-words '("zero" "one" "two" "three" "four" "five" "six" "seven"
		      "eight" "nine" "ten" "eleven" "twelve" "thirteen"
		      "fourteen" "fifteen")))
    (if (<= n 15)
	(lref-default list-words n)
	(number->string n))))

(define (phase-name phase)
  (cond
    ((= phase 0) "New")
    ((= phase 1) "First Quarter")
    ((= phase 2) "Full")
    ((= phase 3) "Last Quarter")
    ('else (error "illegal phase specifier"))))

;;;; Program begins here.

(let ((now (gmtime (unix-time))))
  (phase-of-moon
   (civil->julian (+ 1900 (cdr (assq 'year now)))
                  (+ 1    (cdr (assq 'mon now)))
                          (cdr (assq 'mday now))
                          (cdr (assq 'hour now))
                          (cdr (assq 'min now))
                          (cdr (assq 'sec now)))
   nil))