;;; -*- Mode:Common-Lisp; Package:USER; Base:10 -*-

;;; Sun and Moon position functions.
;;;
;;; Most of the code in this program is adapted from algorithms
;;; presented in "Practical Astronomy With Your Calculator" by
;;; Peter Duffet-Smith.
;;;
;;; The GST and ALT-AZIMUTH algorithms are from Sky and Telescope,
;;; June, 1984 by Roger W. Sinnott.
;;;
;;; Translated from C programs by Robert Bond and Keith E. Brandt.

(DEFCONSTANT *latitude*  32.95       "Default latitude (currently Dallas, Tx.)")
(DEFCONSTANT *longitude* 96.73333333 "Default Longitude (currently Dallas, Tx.)")
(DEFVAR      latitude    *latitude*  "Current latitude")
(DEFVAR      longitude   *longitude* "Current longitude")

(DEFUN sun (&optional (verbose nil) (lat *latitude*) (long *longitude*))
  "Print sunrise and sunset times.  If VERBOSE is non-NIL, it will also
print the sun's azimuth during sunrise and sunset, and the current
position of the sun in altitude and azimuth.  The default values for
latitude and longitude can be changed by specifying LAT and LONG."

  (SETQ latitude lat				; Set global variables
	longitude long)

  (MULTIPLE-VALUE-BIND (seconds minutes hours day month year nil dsp tz)
      (GET-DECODED-TIME)
    (LET* ((tzs (time:timezone-string))
	   (jdate (JulianDate month day year))
	   (ed (- jdate 2444238.5))		; Start in 1980
	   ð1 Æ1 ð2 Æ2
	   trise tset
	   ar as Æ tri da dt)
      (WHEN dsp					; Correct timezone when in
	(DECF tz))				;     daylight savings time

      (MULTIPLE-VALUE-SETQ (ð1 Æ1) (Lon2Eq (SolarLongitude ed)))
      (MULTIPLE-VALUE-SETQ (ð2 Æ2) (Lon2Eq (SolarLongitude (1+ ed))))

      (MULTIPLE-VALUE-BIND (st1r st1s a1r a1s) (RiseSet ð1 Æ1)
	(MULTIPLE-VALUE-BIND (st2r st2s a2r a2s) (RiseSet ð2 Æ2)

	  (WHEN (> (ABS (- st2r st1r)) 1.0)
	    (INCF st2r 24.))

	  (WHEN (> (ABS (- st2s st1s)) 1.0)
	    (INCF st2s 24.))

	  (LET* ((m1 (Adjust24 (- (gmst (- jdate 0.5) (+ 0.5 (/ tz 24.)))
				  (/ longitude 15.))))
		 (ratio (/ (Adjust24 (- st1r m1)) 24.07)))

	    (SETQ trise (Adjust24 (+ (* (- 1.0 ratio) st1r) (* ratio st2r)))
		  ratio (/ (Adjust24 (- st1s m1)) 24.07)
		  tset (Adjust24 (+ (* (- 1.0 ratio) st1s) (* ratio st2s)))))

	  (SETQ ar (/ (* a1r 360.) (+ a1r 360. (- a2r)))
		as (/ (* a1s 360.) (+ a1s 360. (- a2s))))))

      (SETQ Æ (/ (+ Æ1 Æ2) 2.0)
	    tri (acosd (/ (SIND latitude) (COSD Æ))))

      ;; Correction for refraction, parallax, size of sun

      (SETQ da (asind (/ 0.014585145 (tand tri)))
	    dt (* 240. (/ (asind (/ 0.014583594 (SIND tri)))
			  (COSD Æ) 3600.)))

      (MULTIPLE-VALUE-CALL #'FORMAT t "~&Sunrise: ~2d:~2,'0d    "
			   (LstToHm (- trise dt) jdate year tz))
      (WHEN verbose
	(MULTIPLE-VALUE-CALL #'FORMAT t " Azimuth: ~2d deg ~2,'0d min.~%"
			     (DhToHm (- ar da))))
      (MULTIPLE-VALUE-CALL #'FORMAT t "Sunset:  ~2d:~2,'0d    "
			   (LstToHm (+ tset dt) jdate year tz))
      (IF verbose
	  (MULTIPLE-VALUE-CALL #'FORMAT t " Azimuth: ~2d deg ~2,'0d min.~%"
			       (DhToHm (+ as da)))
	  (FORMAT t "(~a)" tzs))
      (WHEN verbose
	(LET ((ð (IF (< ð1 ð2)
			 (/ (+ ð1 ð2) 2)
			 (/ (+ ð1 ð2 24.) 2)))
	      (dh (/ (+ hours (/ minutes 60.0) (/ seconds 3600.0) tz) 24.)))

	  (WHEN (> ð 24.)
	    (DECF ð 24.))
	  (IF (> dh 0.5)
	      (SETQ dh (- dh 0.5)
		    jdate (+ jdate 0.5))
	      (SETQ dh (+ dh 0.5)
		    jdate (- jdate 0.5)))
	  (FORMAT t "The sun is at")
	  (MULTIPLE-VALUE-BIND (alt az) (Eq2AltAz ð Æ (gmst jdate dh))
	    (MULTIPLE-VALUE-CALL #'FORMAT t " altitude: ~2d deg ~2,'0d min," (DhToHm alt))
	    (MULTIPLE-VALUE-CALL #'FORMAT t " azimuth: ~2d deg ~2,'0d min.~%" (DhToHm az)))))
      )
    ))

(DEFUN moon ()
  "Print the current phase of the moon."
  (LET* ((days (- (/ (GET-UNIVERSAL-TIME) 86400.0) 30314))
	 (phase (moon-phase days))
	 (rphase (ROUND phase)))
    (FORMAT t "~&The Moon is ")
    (CASE rphase
      (0
       (FORMAT t "New~%"))
      (50
       (FORMAT t "at the ~a Quarter~%"
	       (IF (> (moon-phase (1+ days)) phase)
		   "First" "Last")))
      (100
       (FORMAT t "full~%"))
      (t
       (FORMAT t "~a ~a (~d% of full)~%"
	       (IF (> (moon-phase (1+ days)) phase) "Waxing" "Waning")
	       (IF (> rphase 50) "Gibbous" "Crescent")
	       rphase)))))


(DEFUN moon-phase (days)
  "Return the current phase of the moon in decimal notation."
  (LET* ((n (Adjust360 (/ (* 360 days) 365.2422)))
	 (smsol (SIND (- (+ n 279.103035) 282.648015)))
	 (Äsol (Adjust360 (+ n (* 1.9155422 smsol) 279.103035)))
	 (l (Adjust360 (+ (* 13.1763966 days) 106.306091)))
	 (mm (Adjust360 (- l (* days 0.1114041) 111.481526)))
	 (ev (* 1.2739 (SIND (- (ASH (- 1 Äsol) 1) mm))))
	 (ac (* 0.1858 smsol))
	 (mmprime (- (+ mm ev) ac (* 0.37 smsol)))
	 (lprime (+ l ev (* 6.2886 (SIND mmprime)) (- ac) (* 0.214 (SIND (ASH mmprime 1))))))
    (* 50 (- 1 (COSD (- (+ lprime (* 0.6583 (SIND (ASH (- lprime Äsol) 1))))
			Äsol))))))

(DEFUN Adjust360 (deg)
  (let ((x (rem deg 360.)))
    (when (minusp x) (incf x 360.))
    x))

(DEFUN Adjust24 (deg)
  (let ((x (rem deg 24.)))
    (when (minusp x) (incf x 24.))
    x))

(DEFUN JulianDate (month day year)
  (WHEN (OR (= month 3) (= month 2))
    (DECF year)
    (INCF month 12.))
  (WHEN (< year 1583)
    (ERROR "Can't handle dates before 1583"))
  (LET ((a (TRUNCATE year 100.)))
    (+ day 1720994.5
       (+ (+ 2 (- a) (TRUNCATE a 4))
	  (TRUNCATE (* year 365.25))
	  (TRUNCATE (* 30.6001 (1+ month)))))))

(DEFUN SolarLongitude (ed)
  (LET* ((m (* (Adjust360
		 (- (Adjust360 (/ (* ed 360.0) 365.2422)) 3.7628784))
	       0.017453292))
	 errt
	 (e (DO ((e m (- e (/ errt (- 1 (* 0.016718 (COS e)))))))
		((<= (SETQ errt (- e (* 0.016718 (SIN e)) m)) 0.0000001)
		 e))))
    (Adjust360 (+ (* 114.59155 (ATAN (* 1.0168601 (TAN (/ e 2)))))
		  282.596403))))

(DEFUN atanqd (y x)
  (LET ((rv 0.0))
    (UNLESS (ZEROP y)
      (SETQ rv (IF (ZEROP x)
		   (IF (> y 0) 90.0 -90.0)
		   (* 57.295776 (ATAN (/ y x))))))
    (IF (MINUSP x)
	(+ rv 180.)
	(IF (MINUSP y)
	    (+ rv 360.)
	    rv))))

(DEFUN Lon2Eq (lambda)
  (LET ((stlam (SIND lambda)))
    (VALUES
      (/ (atanqd (* stlam 0.9174641) (COSD lambda)) 15.0)
      (asind (* 0.39781868 stlam)))))

(DEFUN tand (number)
  "Tangent of an angle measured in degrees."
  (TAN (* number (/ (FLOAT pi number) 180.))))

(DEFUN acosd (number)
  "Arc cosinge of an angle measured in degrees."
  (* (/ 180. (FLOAT pi number)) (ACOS number)))

(DEFUN asind (number)
  "Arc cosinge of an angle measured in degrees."
  (* (/ 180. (FLOAT pi number)) (ASIN number)))

(DEFUN atand (number)
  "Arc tangent of an angle measured in degrees."
  (* (/ 180. (FLOAT pi number)) (ATAN number)))

(DEFUN RiseSet (ð Æ)
  (LET ((tar (/ (SIND Æ) (COSD latitude)))
	h)
    (UNLESS (<= -1.0 tar 1.0)
      (ERROR "The object is circumpolar~%"))
    (SETQ h (/ (acosd (* (- (tand latitude)) (tand Æ))) 15.))
    (VALUES
      (LET ((x (- ð h)))
	(IF (> x 0)
	    x
	    (+ 24. x)))
      (LET ((x (+ ð h)))
	(IF (> x 24.)
	    (- x 24.)
	    x))
      (acosd tar)
      (- 360. (acosd tar)))))

(DEFUN LstToHm (lst jd year tz)
  (LET* ((gmt (+ lst (/ longitude 15.)))
	 (jzjd (JulianDate 1 0 year))
	 (t1 (/ (- jzjd 2415020.0) 36525.0))
	 (t0 (- (* (- jd jzjd) 0.0657098)
		(- 24 (- (+ 6.6460656 (* 2400.05126 t1)
			    (* 2.58e-5 t1 t1))
			 (* 24 (- year 1900)))))))
    (WHEN (> gmt 24.)
      (DECF gmt 24.))
    (WHEN (MINUSP t0)
      (INCF t0 24.))
    (SETQ gmt (- gmt t0))
    (WHEN (MINUSP gmt)
      (INCF gmt 24.))
    (SETQ gmt (- (* gmt 0.99727) tz))
    (WHEN (MINUSP gmt)
      (INCF gmt 24.))
    (DhToHm gmt)))

(DEFUN DhToHm (dh)
  "Convert DH from decimal hour notation to hours and minutes"
  (MULTIPLE-VALUE-BIND (hour minute) (TRUNCATE (+ dh 4.166667e-4))
    (SETQ minute (TRUNCATE (* 60. minute)))
    (VALUES hour minute)))

(DEFUN Eq2AltAz (r d t0)
  (SETQ r  (* r 0.2617994)
	d  (* d 0.017453292)
	t0 (* t0 0.2617994))
  (LET* ((b (* latitude 0.017453292))
	 (t5 (+ (- t0 r) (* (- 360 longitude) 0.017453292)))
	 (s1 (+ (* (SIN b) (SIN d)) (* (COS b) (COS d) (COS t5))))
	 (c1 (- 1 (* s1 s1)))
	 (s2 (- (* (COS d) (SIN t5))))
	 (c2 (- (* (COS b) (SIN d)) (* (SIN b) (COS d) (COS t5))))
	 a h)
    (IF (> c1 0)
	(SETQ c1 (SQRT c1)
	      h (ATAN (/ s1 c1)))
	(SETQ h (* (SIGNUM s1) 1.5707964)))
    (IF (ZEROP c2)
	(SETQ a (* (SIGNUM s2) 1.5707964))
	(SETQ a (ATAN (/ s2 c2)))
	(WHEN (MINUSP c2)
	  (INCF a 3.1415927)))
    (WHEN (MINUSP a)
      (INCF a 6.2831855))
    (VALUES
      (/ h 0.017453292)
      (/ a 0.017453292))))

(DEFUN gmst (j f)
  (LET* ((t0 (/ (- j 2451545.0) 36525.0))
	 (t1 (FLOOR t0))
	 (t2 (/ (- j (+ (* t1 36525.0) 2451545.0) -0.5) 36525.0))
	 (s  (* 24. (+ (* (- f 0.5) 1.002737909)
		       (REM (/ (+ 24110.54841 (* 184.812866 t1)
				  (* 8640184.812866 t2)
				  (* t0 t0
				     (- 0.093104 (* 0.0000062 t0))))
			       86400.) 1)))))
    (IF (MINUSP s)
	(+ s 24.)
	(IF (> s 24.)
	    (- s 24.)
	    s))))
