;;; --------------------------------------------------------------------------- ;;; MoonPhase.lsp ;;; code compiled by YZ August 2026 ;;; --------------------------------------------------------------------------- ;;; PHASE OF THE MOON, CALCULATED AND DRAWN ;;; ;;; PURPOSE ;;; Works out the phase of the moon for today or any date you give, reports it, ;;; and draws the moon at that phase. ;;; ;;; The drawing is the novelty. The CALCULATION is the substance: it finds the ;;; times of the new moons either side of the date by summing periodic terms ;;; in the moon's and sun's mean anomalies, and interpolates between them. It ;;; is good to about two minutes. ;;; ;;; WHY THAT IS HARDER THAN IT LOOKS ;;; The lunar month is not constant. It averages 29.53 days but varies by ;;; several hours either side, because the moon's orbit is elliptical and the ;;; sun perturbs it. Dividing the days since a known new moon by 29.53 is ;;; wrong by up to most of a day - which is the difference between a crescent ;;; and a quarter. ;;; ;;; So the series below carries terms in: ;;; ;;; M5 the sun's mean anomaly ;;; M6 the moon's mean anomaly ;;; B6 the moon's argument of latitude ;;; ;;; with secular corrections in T, the centuries since 1899.5. The largest ;;; correction, the -0.4068 sin(M6) term, is worth nearly ten hours on its own. ;;; ;;; WHAT GETS DRAWN ;;; The disc as a circle, and the terminator - the line between light and dark ;;; - as an ellipse across it. The terminator's width is the cosine of the ;;; phase angle, which is what makes it a straight line at the quarters, a ;;; narrow ellipse near new and full, and correctly bowed either way between. ;;; ;;; The routine this came from drew the same shape by making a block of a ;;; semicircular arc and re-inserting it with an X scale, which left a block ;;; definition named after the time of day in every drawing it touched. ;;; ;;; MOONPHASE - report and draw the phase of the moon ;;; --------------------------------------------------------------------------- ;;; --------------------------------------------------------------------------- ;;; SESSION MEMORY ;;; --------------------------------------------------------------------------- (if (null *MoonPhase:Prefs*) (setq *MoonPhase:Prefs* (list (cons "TZONE" 0.0) ; hours to add to UT for local time (cons "SIZE" 1.0) ; last radius used ) ) ) (defun MoonPhase:Get ( key ) (cdr (assoc key *MoonPhase:Prefs*))) (defun MoonPhase:Put ( key val ) (setq *MoonPhase:Prefs* (cons (cons key val) (vl-remove-if '(lambda (p) (= (car p) key)) *MoonPhase:Prefs*))) val ) ;;; --------------------------------------------------------------------------- ;;; JULIAN DATE TO CALENDAR YEAR ;;; ;;; Only the year is needed - it seeds the lunation count - but the month and day ;;; fall out of the same reduction and are returned for the report. ;;; --------------------------------------------------------------------------- (defun MoonPhase:JYear ( jd / j y d m ) (setq j (- (fix jd) 1721119.0) y (fix (/ (1- (* 4 j)) 146097.0)) j (- (* j 4.0) 1.0 (* 146097.0 y)) d (fix (/ j 4.0)) j (fix (/ (+ (* 4.0 d) 3.0) 1461.0)) d (- (+ (* 4.0 d) 3.0) (* 1461.0 j)) d (fix (/ (+ d 4.0) 4.0)) m (fix (/ (- (* 5.0 d) 3) 153.0)) d (- (* 5.0 d) 3.0 (* 153.0 m)) d (fix (/ (+ d 5.0) 5.0)) y (+ (* 100.0 y) j)) ;; The reduction counts March as month 1, so anything from January or ;; February belongs to the following calendar year. (if (< m 10.0) (setq m (+ m 3)) (setq m (- m 9) y (1+ y))) (list y m d) ) ;;; --------------------------------------------------------------------------- ;;; PHASE ;;; ;;; Returns ( phase age-in-days lunation-length ) where phase runs 0 to 1: ;;; ;;; 0.00 new 0.50 full ;;; 0.25 first quarter 0.75 last quarter ;;; ;;; Every working variable is local. The routine this came from used TT's name ;;; for the centuries term - T, LISP's true constant - and left it overwritten ;;; for the rest of the session. ;;; --------------------------------------------------------------------------- (defun MoonPhase:Phase ( jd tzone / year rad k0 tt tt2 tt3 j0 f0 m0 m1 b1 k9 j f k m5 m6 b6 prev this found ) (setq year (car (MoonPhase:JYear jd)) rad (/ pi 180.0) ;; Lunations since 1900. 12.3685 is the mean number per year. k0 (fix (* (- year 1900.0) 12.3685)) ;; Centuries since 1899.5, and its powers - the secular terms below ;; are polynomials in it. tt (/ (- year 1899.5) 100.0) tt2 (* tt tt) tt3 (* tt2 tt) ;; Mean new moon: integer part, then the fractional day. j0 (+ 2415020.0 (* 29 k0)) f0 (+ (- (* 0.0001178 tt2) (* 1.55e-7 tt3)) 0.75933 (* 0.53058868 k0)) f0 (- f0 (* 8.370001e-4 tt) (* 0.000335 tt2)) ;; Sun's mean anomaly at that new moon. m0 (* k0 0.08084821133) m0 (+ (* 360.0 (- m0 (fix m0))) 359.2242) m0 (- m0 (* 0.0000333 tt2) (* 3.47e-6 tt3)) ;; Moon's mean anomaly. m1 (* k0 0.07171366128) m1 (+ (* 360.0 (- m1 (fix m1))) 306.0253 (* 0.0107306 tt2) (* 1.236e-5 tt3)) ;; Moon's argument of latitude. b1 (* k0 0.08519585128) b1 (+ (* 360.0 (- b1 (fix b1))) 21.2964) b1 (- b1 (* 0.0016528 tt2) (* 2.39e-6 tt3)) k9 0.0 prev 0.0 found nil) ;; Step forward two lunations at a time until the date is bracketed. The ;; half-steps are what let the loop land on quarters as well as new moons. (while (and (not found) (< k9 29)) (setq j (+ j0 (* 14.0 k9)) f (+ f0 (* 0.765294 k9)) k (/ k9 2.0) m5 (* (+ m0 (* k 29.10535608)) rad) m6 (* (+ m1 (* k 385.81691806)) rad) b6 (* (+ b1 (* k 390.67050646)) rad) ;; The periodic corrections. The first is by far the largest - ;; nearly ten hours at its extremes. f (+ f (* -0.4068 (sin m6)) (* (- 0.1734 (* 0.000393 tt)) (sin m5)) (* 0.0161 (sin (* 2 m6))) (* 0.0104 (sin (* 2 b6))) (* -0.0074 (sin (- m5 m6))) (* -0.0051 (sin (+ m5 m6))) (* 0.0021 (sin (* 2 m5))) (* 0.001 (sin (- (* 2 b6) m6)))) j (+ j (fix f)) f (- f (fix f)) f (+ f (/ tzone 24.0))) ;; Carry the fractional day into the integer day if it has run over. (if (>= f 1.0) (setq f (1- f) j (1+ j))) (if (< f 0.0) (setq f (1+ f) j (1- j))) (setq prev this this (+ j f)) ;; Bracketed: the date lies between these two phase events. (if (and prev (>= jd prev) (< jd this)) (setq found t) (setq k9 (+ k9 2)))) (if (and found prev (> (- this prev) 0.0)) (list (/ (- jd prev) (- this prev)) (- jd prev) (- this prev))) ) ;;; The conventional name for a phase value. (defun MoonPhase:Name ( ph ) (cond ((< ph 0.03) "New") ((< ph 0.22) "Waxing crescent") ((< ph 0.28) "First quarter") ((< ph 0.47) "Waxing gibbous") ((< ph 0.53) "Full") ((< ph 0.72) "Waning gibbous") ((< ph 0.78) "Last quarter") ((< ph 0.97) "Waning crescent") (t "New") ) ) ;;; --------------------------------------------------------------------------- ;;; MAIN COMMAND ;;; --------------------------------------------------------------------------- (defun c:MOONPHASE ( / *error* vars vals jd tzone res ph age lun ymd cen size lit ratio v ) (setq vars '("CMDECHO" "OSMODE" "BLIPMODE" "CLAYER") vals (mapcar 'getvar vars)) (defun MoonPhase:Restore ( ) (mapcar 'setvar vars vals) (while (= 8 (logand 8 (getvar 'undoctl))) (command "_.UNDO" "_End")) (vl-catch-all-apply '(lambda ( ) (*pop-error-mode*)) '()) (princ) ) (defun *error* ( msg ) (MoonPhase:Restore) (if (and msg (not (wcmatch (strcase msg t) "*break*,*cancel*,*exit*"))) (princ (strcat "\n** MOONPHASE error: " msg " **"))) (princ) ) (setvar "CMDECHO" 0) ;; AutoCAD 2015 and later refuse (command) inside an *error* handler ;; unless the routine says up front that it will use one. Restore does, ;; to close this undo group. (vl-catch-all-apply '(lambda ( ) (*push-error-using-command*)) '()) (command "_.UNDO" "_Begin") ;; Time zone, as hours to add to UT. The routine this came from had 8 hard ;; coded, which is right for exactly one part of the world. (initget 4) (setq v (getreal (strcat "\nHours to add to UT for local time <" (rtos (MoonPhase:Get "TZONE") 2 2) ">: "))) (if v (MoonPhase:Put "TZONE" v)) (setq tzone (MoonPhase:Get "TZONE")) (setq jd (getvar "DATE") ymd (MoonPhase:JYear jd) res (MoonPhase:Phase jd tzone)) (if (null res) (princ "\n** Could not bracket the date between two phase events. **") (progn (setq ph (car res) age (cadr res) lun (caddr res) ;; Illuminated fraction: 0 at new, 1 at full, 0.5 at quarters. lit (/ (- 1.0 (cos (* 2.0 pi ph))) 2.0)) (princ (strcat "\n" (MoonPhase:Name ph) " - " (rtos (* lit 100.0) 2 1) "% lit" "\n Age " (rtos (float (fix age)) 2 0) " days " (rtos (* 24.0 (- age (fix age))) 2 0) " hours" "\n This lunation " (rtos lun 2 4) " days" "\n Date " (itoa (fix (car ymd))) "-" (itoa (fix (cadr ymd))) "-" (itoa (fix (caddr ymd))))) (initget "Yes No") (if (/= "No" (getkword "\nDraw the moon [Yes/No] : ")) (progn (setq cen (getpoint "\nCentre of the moon: ")) (if cen (progn (initget 6) (setq v (getdist cen (strcat "\nRadius <" (rtos (MoonPhase:Get "SIZE") 2 3) ">: "))) (if v (MoonPhase:Put "SIZE" v)) (setq size (MoonPhase:Get "SIZE")) (setvar "OSMODE" 0) (setvar "BLIPMODE" 0) (if (tblsearch "layer" "MOON") (command "_.LAYER" "_ON" "MOON" "_THAW" "MOON" "_UNLOCK" "MOON" "") (command "_.LAYER" "_NEW" "MOON" "")) (setvar "CLAYER" "MOON") ;; The disc. (command "_.CIRCLE" cen size) ;; The terminator. Its half-width across the disc is ;; the cosine of the phase angle, so it collapses to a ;; straight line at the quarters. (setq ratio (abs (cos (* 2.0 pi ph)))) (if (< ratio 1e-6) (command "_.LINE" (polar cen (* pi 0.5) size) (polar cen (* pi 1.5) size) "") (entmake (list '(0 . "ELLIPSE") '(100 . "AcDbEntity") (cons 8 "MOON") '(100 . "AcDbEllipse") (cons 10 cen) ;; Major axis vertical, so the ellipse ;; narrows across the disc as the phase ;; approaches a quarter. (cons 11 (list 0.0 size 0.0)) '(210 0.0 0.0 1.0) (cons 40 ratio) (cons 41 0.0) (cons 42 (* 2.0 pi))))) (princ "\nDrawn.")))))) ) (MoonPhase:Restore) (princ) ) (princ)