(in-package "COMMON-LISP-USER")

(defpackage #:cl-libtai
  (:use #:common-lisp)
  (:export :leapsecs-add :*leapsecs* :*leapsecs-strings*
	   :tai-unpack
	   :taia-unpack
	   :tai-now
	   :caldate-scan
	   :caldate-mjd
	   :caldate-frommjd
	   :caltime-tai
	   :caltime-utc
	   :leapsecs-gen
	   :leapsecs-add
	   :leapsecs-sub
	   :make-caldate
	   :make-caltime
	   :make-tai64-internal
	   :make-tai64na-internal))

(in-package "CL-LIBTAI")

(defvar times365 (list 0 365 730 1095))
(defvar times36524 (list 0 36524 73048 109572))
(defvar montab (list 0 31 61 92 122 153 184 214 245 275 306 337))

(defun make-tai64-internal(s)
  (pairlis '(x) (list s)))

(defun make-tai64na-internal(s a n)
  (pairlis '(sec atto nano) (list s a n)))

(defun make-caldate(y m d)
  (list :year y :month m :day d))

(defun make-caltime(cd h m s o)
  (pairlis '(date hour minute second offset) (list cd h m s o)))

; http://maia.usno.navy.mil/leapsec.html
; http://maia.usno.navy.mil/ser7/leapsec.dat

(defvar *leapsecs-strings* (list 
"+1972-06-30"
"+1972-12-31"
"+1973-12-31"
"+1974-12-31"
"+1975-12-31"
"+1976-12-31"
"+1977-12-31"
"+1978-12-31"
"+1979-12-31"
"+1981-06-30"
"+1982-06-30"
"+1983-06-30"
"+1985-06-30"
"+1987-12-31"
"+1989-12-31"
"+1990-12-31"
"+1992-06-30"
"+1993-06-30"
"+1994-06-30"
"+1995-12-31"
"+1997-06-30"
"+1998-12-31"
"+2005-12-31"))

(defun leapsecs-add(t2 hit leapsecs)
  (let ((u (cdr (assoc 'x t2))))
    (loop for a in leapsecs for m = (cdr (assoc 'x a))
	do (if (< u m) (loop-finish)
		           (if (or (= hit 0) (> u m)) (setf u (1+ u))))
	finally (return (make-tai64-internal u)))))

(defun leapsecs-sub(t2 leapsecs)
  (let ((u (cdr (assoc 'x t2))))
    (loop for a in leapsecs for s = 0 then (1+ s)
	do (if (< u (cdr (assoc 'x a))) (return (values (make-tai64-internal (- u s)) 0)))
	   (if (= u (cdr (assoc 'x a))) (return (values (make-tai64-internal (- u (1+ s))) 1)))
	finally (return (values (make-tai64-internal (- u s)) 0)))))

(defun tai-now()
  (make-tai64-internal (+ (- 4611686018427387914 2208988800) (get-universal-time))))

(defun tai-pack(tai64i)
  (format nil "~x" (cdr (assoc 'x tai64i))))

(defun tai-unpack(s)
  (make-tai64-internal (parse-integer s :radix 16)))

(defun taia-pack(tai64nai)
  (format nil "~x~x~x" (cdr (assoc 'sec tai64nai))
			    (cdr (assoc 'atto tai64nai))
				 (cdr (assoc 'nano tai64nai))))

(defun taia-unpack(s)
  (let ((tai-s (subseq s 0 8)) (atto-s (subseq s 8 12)) (nano-s (subseq s 12 16 )))
    (make-tai64na-internal (parse-integer tai-s :radix 16)
		 (parse-integer atto-s :radix 16)
		 (parse-integer nano-s :radix 16))))

(defun display-numeric(tai-s atto-s nano-s)
  (format t "~a ~a ~a" (parse-integer tai-s :radix 16) 
	  (parse-integer atto-s :radix 16)
	  (parse-integer nano-s :radix 16)))

(defun caldate-mjd (cd)
  (defun final-yd(y d)
    	(+ d 
	   (nth (logand #x3 y) times365) 
	   (* 1461 (rem (floor (/ y 4)) 25))
	   (nth (logand #x3 (floor (/ y 100))) times36524)))
  (defun mjd-ycheck1(year m day)
    (let ((y (rem year 400)) (d (+ (* (floor (/ year 400)) 146097) (+ day (+ (nth m montab))))))
      (if (< y 0) (final-yd (+ y 400) (- d 146097)) (final-yd y d))))
  (defun mcheck2(year month d)
    (let ((y (+ year (floor (/ month 12)))) (m (rem month 12)))
      (if (< m 0) (mjd-ycheck1 (- y 1) (+ m 12) d) (mjd-ycheck1 y m d))))
  (defun mcheck1(y m d)
	(if (>= m 2) (mcheck2 y (- m 2) d) (mcheck2 (- y 1) (+ m 10) d)))
  (destructuring-bind (&key year month day) cd 
    (mcheck1 (rem year 400) (- month 1) (+ (* 146097 (floor (/ year 400))) (- day 678882)))))

(defun caldate-frommjd(day)
  (let ((pwday 0))
  (defun daydec(n acc)
    (if (>= n 146097) (daydec (- n 146097) (1+ acc)) (values n acc)))
  (defun check5(yday y m d)
    (values (make-caldate y (+ m 1) (+ d 1)) pwday yday))
  (defun check4(yday year day)
    (let ((da (* 10 day)))
      (let ((m (floor (/ (+ 5 da) 306))) (d (floor (/ (rem (+ da 5) 306) 10))))
	(if (>= m 10) (check5 (- yday 306) (+ year 1) (- m 10) d) (check5 (+ yday 59) year (+ m 2) d)))))
  (defun check3(year day)
    (let ((yday (if (< day 306) 1 0)))
      (if (eq day 1460) (check4 yday (+ year 3) 365) (check4 yday (+ year (floor (/ day 365))) (rem day 365)))))
  (defun check2(year day)
    (let ((y (* 4 (+ (floor (/ day 1461)) (* year 25)))) (d (rem day 1461)))
      (check3 y d)))
  (defun fmjd-ycheck1(year day)
    (let ((y (* year 4)))
      (if (eq day 146096) (check2 (+ y 3) 36524) (check2 (+ y (floor (/ day 36524))) (rem day 36524)))))
  (let ((year (floor (/ day 146097))) (d (+ 678881 (rem day 146097))))
    (multiple-value-bind (newday yeardiff)
	  (daydec d 0)
      (progn (setf pwday (mod (+ 3 newday) 7))
	(fmjd-ycheck1 (+ year yeardiff) newday))))))

(defun caldate-normalize(cd)
  (caldate-frommjd (caldate-mjd cd)))

(defun caldate-scan(s)
  (let* ((sign (if (eq (char s 0) #\-) -1 1)) (r (string-trim "-+" s)) (p1 (position #\- r)) (p2 (position #\- r :start (1+ p1))))
    (make-caldate (* (parse-integer r :end p1) sign) (parse-integer r :start (1+ p1) :end p2) (parse-integer r :start (1+ p2)))))

(defun leapsecs-gen(ls)
  (loop for a in ls for leaps = 0 then (1+ leaps) 
	collect (make-tai64-internal (+ leaps (+ 4611686014920671114 (* 86400 (1+ (caldate-mjd (caldate-scan a)))))))))

(defparameter *leapsecs* (leapsecs-gen *leapsecs-strings*))

(defun caltime-tai(ct)
  (defun m60p(mul syp)
    (+ (cdr (assoc syp ct)) (* 60 mul)))
  (defun inner() (m60p (cdr (assoc 'hour ct)) 'minute))
  (defun xval(d)
    (+ 4611686014920671114 (* 86400 d) (m60p (- (inner) (cdr (assoc 'offset ct))) 'second)))
  (let ((day (caldate-mjd (cdr (assoc 'date ct)))))
    (leapsecs-add (make-tai64-internal (xval day)) (if (= 60 (cdr (assoc 'second ct))) 1 0) *leapsecs*)))

(defun caltime-utc(tai64i)
  (multiple-value-bind (tn leap)
      (leapsecs-sub tai64i *leapsecs*)
  (let ((u (+ 58486 (cdr (assoc 'x tn)))))
    (let ((s (rem u 86400)))
      (make-caltime (caldate-frommjd (logand #xFFFFFFFF (- (floor (/ u 86400)) 53375995543064))) (floor (/ s 3600)) (rem (floor (/ s 60)) 60) (+ leap (rem s 60)) 0)))))