;; Mathematics on timepoints.
;; Liam Healy Fri Oct 12 2001 - 16:32
;; Time-stamp: <2014-10-06 23:10:23EDT dtmath.lisp>
;; Copyright 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program. If not, see .
(in-package :antik)
(export '(timepoint-j2000sec make-timepoint-J2000sec
datime+ datime- days+ days-
start-of-day seconds-per-day datime-J2000day))
;;; Mathematics can only be done on timepoints, not on dtspecs,
;;; because the scale must be known first.
;;; Keep in mind that CL's *-universal-time are handy for
;;; linearizing and delinearizing a specification to do math,
;;; but they carry no scale; they are simply a dtspec.
(defvar *j2000-epoch* (read-timepoint "2000-01-01T12:00:00"))
(defvar *clut-j2000-tai*
(timeparse-to-clut
(timeparse (convert-time-scale *j2000-epoch* :tai)))
"J2000 epoch in CL's `universal-time' time linearization, in the TAI scale.")
(defun timepoint-j2000sec (tp)
"Seconds from the J2000 epoch, 2000-01-01T12:00:00UTC."
(flet ((tp-to-clut (tp)
(timeparse-to-clut (timeparse tp))))
(check-type tp timepoint)
(cl:- (tp-to-clut (convert-time-scale tp :tai)) *clut-j2000-tai*)))
(defun make-timepoint-J2000sec (J2000sec &optional day-only)
"Make a timepoint from J2000sec."
(convert-time-scale
(make-timepoint
:timeparse
(clut-to-timeparse (cl:+ *clut-j2000-tai* J2000sec))
:scale :tai
:day-only day-only)
:utc))
;;;****************************************************************************
;;; Whole days
;;;****************************************************************************
(defun start-of-day
(datime &optional (relative-day 0) (day-only (day-only datime)))
"From the given datime, return the start of the given day,
relative-day = 0: this day, -1: day before, +1 next day, etc."
(let ((answer
(dtspec-from-julian-day-number
(cl:+ (julian-day-number datime t) relative-day)
t)))
(change-class answer (class-of datime))
(if (typep datime 'timepoint)
(setf (slot-value answer 'scale) (scale datime)))
(if (or day-only (day-only datime))
(setf (slot-value answer 'day-only) day-only))
answer))
;;; (seconds-per-day #d1998-12-31) -> 86401
(defun seconds-per-day (datime)
"Compute the number of seconds in the given day."
(pq-magnitude
(datime- (start-of-day datime 1) (start-of-day datime))))
(defun datime-J2000day (datime)
"Julian Day from J2000 epoch, and seconds from the start of the day."
(values
;; Days from J2000 epoch to the start of the UTC day; always integer + 0.5.
(cl:- (days- datime (MAKE-TIMEPOINT :TIMEPARSE '(0 0 12 1 1 2000) :DAY-ONLY NIL :SCALE :UTC))
1/2)
;; Seconds from the start of the UTC day to the current timepoint.
;; No timescale is required because this is
;; only a difference from the start of the day.
(cl:- (timeparse-to-clut (timeparse datime)) ;(timepoint-J2000sec datime)
(timeparse-to-clut (timeparse (start-of-day datime))))))
(defun days- (datime1 datime2)
"Integer number of days separating the days. Differs from datime-
in that the time of day doesn't matter, just the day on which the two
datimes occur."
(cl:- (julian-day-number datime1) (julian-day-number datime2)))
(defun days+ (datime &optional (increment 1))
"Add the number of days to the given datime."
(start-of-day datime increment))
;;;****************************************************************************
;;; Addition, subtraction
;;;****************************************************************************
(defun datime+ (tp time)
"Add the time to the timepoint. If time is a number,
it is assumed to represent seconds."
(check-type tp timepoint)
(with-units
(with-pq ((time time))
(make-timepoint-j2000sec
(cl:+ (timepoint-J2000sec tp)
(if (numberp time)
time
(pq-magnitude time)))
(day-only tp)))))
;;; N.B.:
;;; (datime- #d1999-01-01T12:00:00 #d1998-12-31T12:00:00)
;;; 86401.0 second
(defun datime- (datime1 b)
"The difference of the datimes, or subtract a time from a datime."
(if (typep b 'timepoint)
;; add conditional to check if either are day-only
;; if so, just count the number of days.
(if (or (day-only datime1) (day-only b))
(make-pq (days- datime1 b) 'day)
(make-pq ; time difference between two datimes
(cl:- (timepoint-J2000sec datime1)
(timepoint-J2000sec b))
'second))
(datime+ datime1 (antik:- b)))) ; subtract a time
(defmethod +i ((a dtspec) b)
(datime+ a b))
(defmethod +i (a (b dtspec))
(datime+ b a))
(defmethod -i ((a dtspec) b)
(datime- a b))
;;;****************************************************************************
;;; Relational
;;;****************************************************************************
(defun timepoint-relation (relation x y)
"Do the two timepoints obey the relation?"
(datime-relation
relation
(multiple-value-bind (converted error)
(ignore-errors (convert-time-scale y (scale x)))
(if (typep error 'error) ; can't convert y to x type?
;; then convert x to y type
(list (convert-time-scale x (scale y)) y)
(list x (if (eq (scale x) (scale y)) y converted))))))
;;; define <, etc.
(defmethod i ((x timepoint) (y timepoint)) (timepoint-relation '> x y))
(defmethod <=i ((x timepoint) (y timepoint)) (timepoint-relation '<= x y))
(defmethod >=i ((x timepoint) (y timepoint)) (timepoint-relation '>= x y))
(defmethod =i ((x timepoint) (y timepoint))
(if (timepoint-relation '= x y) x))