ramsdell@linus.UUCP (John D. Ramsdell) (08/26/85)
;;; -*- Mode: LISP; Base: 10; Package: COMMON-LISP-GLOBAL; Syntax: Common-lisp -*-
; Functions that compute the Julian date
; and the position of the sun and moon.
; Use print-julian-date and print-astro-time.
(provide 'astro-time)
(in-package 'astro-time)
(export '(print-julian-date print-astro-time))
;J2000.0 = January 1, 2000; Noon GMT.
(defconstant J2000.0 (encode-universal-time 0 0 12 1 1 2000 0))
(defconstant julian-date-at-j2000.0 2451545.0d0)
(defconstant time-units-per-day (* 60 60 24))
(defconstant days-per-time-unit (/ (float time-units-per-day)))
(defconstant days-per-julian-century 36525.0)
(defun time-units-from-j2000.0 (ut) ; Should be declared inline.
(float (- ut j2000.0)))
(defun ut->jd (ut) ; => Julian date as a dfloat.
(+ julian-date-at-j2000.0
(float (* days-per-time-unit
(time-units-from-j2000.0 ut))
0.0d0)))
(defun print-julian-date (&optional (ut (get-universal-time))
(stream *standard-output*))
(format stream "JD~2,7,,$" (ut->jd ut)))
; Radians and degrees.
(defconstant -pi (- pi))
(defconstant 2pi (* 2.0 pi))
(defconstant radians-per-degree (/ pi 180.0))
; Reduce an angle so that -pi < angle <= pi.
(defun normalize-angle (angle)
(if (> angle pi)
(loop while (> angle pi) do (setq angle (- angle 2pi)))
(loop while (<= angle -pi) do (setq angle (+ angle 2pi))))
angle)
; Astronomical almanac
; All formulas from
; The Astronomical Almanac for the Year 1984,
; US Naval Observatory and Royal Greenwich Observatory,
; US Goverment Printing Office, Washington DC, 1984.
; Position of the sun. (Page C24).
; Convert degrees to radians at compile time.
(defconstant sun0 (* radians-per-degree 280.460)) ; for mean longitude
(defconstant sun1 (* radians-per-degree 0.9856474 days-per-time-unit))
(defconstant sun2 (* radians-per-degree 357.528)) ; for mean anomaly
(defconstant sun3 (* radians-per-degree 0.9856003 days-per-time-unit))
(defconstant sun4 (* radians-per-degree 1.915)) ; for eliptic longitude
(defconstant sun5 (* radians-per-degree 0.020))
(defun sun-position (ut) ; Input: universal time.
"Returns the ecliptic longitude and the mean anomaly of the sun."
(let* ((time (time-units-from-j2000.0 ut))
(mean-longitude-of-sun (normalize-angle (+ sun0 (* sun1 time))))
(mean-anomaly (normalize-angle (+ sun2 (* sun3 time))))
(eliptic-longitude
(normalize-angle
(+ mean-longitude-of-sun
(* sun4 (sin mean-anomaly))
(* sun5 (sin (* 2.0 mean-anomaly)))))))
(values eliptic-longitude mean-anomaly)))
; Position of the moon (Page D46).
(defconstant julian-centuries-per-time-unit
(/ days-per-time-unit days-per-julian-century))
(defconstant moon0 (* radians-per-degree 218.32))
(defconstant moon1 (* radians-per-degree 481267.883 julian-centuries-per-time-unit))
(defconstant moon2a (* radians-per-degree 6.29))
(defconstant moon2b (* radians-per-degree 134.9))
(defconstant moon2c (* radians-per-degree 477198.85 julian-centuries-per-time-unit))
(defconstant moon3a (* radians-per-degree -1.27))
(defconstant moon3b (* radians-per-degree 259.2))
(defconstant moon3c (* radians-per-degree -413335.38 julian-centuries-per-time-unit))
(defconstant moon4a (* radians-per-degree 0.66))
(defconstant moon4b (* radians-per-degree 235.7))
(defconstant moon4c (* radians-per-degree 890534.23 julian-centuries-per-time-unit))
(defconstant moon5a (* radians-per-degree 0.21))
(defconstant moon5b (* radians-per-degree 269.9))
(defconstant moon5c (* radians-per-degree 954397.70 julian-centuries-per-time-unit))
(defconstant moon6a (* radians-per-degree -0.19))
(defconstant moon6b (* radians-per-degree 357.5))
(defconstant moon6c (* radians-per-degree 035999.05 julian-centuries-per-time-unit))
(defconstant moon7a (* radians-per-degree -0.11))
(defconstant moon7b (* radians-per-degree 186.6))
(defconstant moon7c (* radians-per-degree 966404.05 julian-centuries-per-time-unit))
(defun moon-position (ut) ; Input: universal time.
"Returns the eliptic longitude of the moon."
(let ((time (time-units-from-j2000.0 ut)))
(normalize-angle
(+ moon0
(* moon1 time)
(* moon2a (sin (+ moon2b (* moon2c time))))
(* moon3a (sin (+ moon3b (* moon3c time))))
(* moon4a (sin (+ moon4b (* moon4c time))))
(* moon5a (sin (+ moon5b (* moon5c time))))
(* moon6a (sin (+ moon6b (* moon6c time))))
(* moon7a (sin (+ moon7b (* moon7c time))))))))
(defun moon-and-sun-position (ut)
" => moon and sun position, mean-anomaly in radians."
(multiple-value-bind (sun-eliptic-longitude mean-anomaly)
(sun-position ut)
(let ((relative-moon-eliptic-longitude
(normalize-angle (- (moon-position ut)
sun-eliptic-longitude))))
(values relative-moon-eliptic-longitude
sun-eliptic-longitude
mean-anomaly))))
; Prints the Julian date and the position of the sun and moon.
(defun print-astro-time (&optional (ut (get-universal-time))
(stream *standard-output*))
(multiple-value-bind (relative-moon-eliptic-longitude sun-eliptic-longitude mean-anomaly)
(moon-and-sun-position ut)
(format stream
"~&Date JD~1,7,,$,
Moon phase ~1,2,6,$% (0.0% at new moon),
Sun phase ~1,2,6,$% (0.0% at spring equinox),
Sun anomaly ~1,2,6,$% (0.0% at perigee)."
(ut->jd ut)
(* (/ relative-moon-eliptic-longitude pi) 100.0)
(* (/ sun-eliptic-longitude pi) 100.0)
(* (/ mean-anomaly pi) 100.0))))