more

src/calendar.hpp

@@ -0,0 +1,18 @@
+#ifndef __DSO_CALENDAR_GENINC_HPP__
+#define __DSO_CALENDAR_GENINC_HPP__
+
+#include "dtdatetime.hpp"
+#include "tpdate.hpp"
+#include "tpdate2.hpp"
+
+namespace dso {
+
+/** alias for datetime<T> */
+template <typename T> using Datetime = datetime<T>;
+
+/** alias for TwoPartDate */
+using MjdEpoch = TwoPartDate;
+
+} /* namespace dso */
+
+#endif

src/datetime_read.hpp

@@ -7,7 +7,7 @@
 #define __DSO_DATETIME_IO_READ_HPP__
 
 #include "datetime_io_core.hpp"
-#include "dtcalendar.hpp"
+#include "dtdatetime.hpp"
 #include "hms_time.hpp"
 #include "tpdate.hpp"
 #include <stdexcept>

src/datetime_write.hpp

@@ -7,7 +7,7 @@
 #define __DSO_DATETIME_IO_WRITE_HPP__
 
 #include "datetime_io_core.hpp"
-#include "dtcalendar.hpp"
+#include "dtdatetime.hpp"
 #include "hms_time.hpp"
 #include "tpdate.hpp"
 #include <cstdio>

src/dtcalendar.hpp → src/dtdatetime.hpp

@@ -282,9 +282,20 @@ class datetime {
     return datetime(modified_julian_day::min(), S(0));
   }
 
+  /** Reference epoch (J2000.0), as a Modified Julian Date. */
+  constexpr static datetime j2000_mjd() noexcept {
+    return datetime(modified_julian_day(51544), S(S::max_in_day / 2L));
+  }
+
   /** Default constructor. */
   explicit constexpr datetime() noexcept : m_mjd(dso::J2000_MJD), m_sec(0){};
 
+  double jcenturies_sinceJ2000() const noexcept {
+    const double d_mjd = (double)(m_mjd.as_underlying_type());
+    const double fdays = fractional_days(m_sec);
+    return ((d_mjd - J2000_MJD) + fdays) / DAYS_IN_JULIAN_CENT;
+  }
+
   /** Constructor from year, month, day of month and sec type.
    * If an invalid date is passed-in, the constructor will throw.
    */
@@ -521,14 +532,6 @@ class datetime {
     return m_sec.fractional_days() + jd;
   }
 
-  /** @brief compute Julian centuries since J2000 */
-  constexpr double jcenturies_sinceJ2000() const noexcept {
-    const double jd = m_mjd.to_julian_day();
-    const double jc =
-        (jd - J2000_JD) / 36525e0 + m_sec.fractional_days() / 36525e0;
-    return jc;
-  }
-
   /** @brief Cast to year, month, day of month */
   constexpr ymd_date as_ymd() const noexcept { return m_mjd.to_ymd(); }
 

src/dtfund.hpp

@@ -899,20 +899,20 @@ class modified_julian_day {
    */
   constexpr underlying_type &__member_ref__() noexcept { return m_mjd; }
 
-  /** @brief Max possible modified_julian_day 
-   * 
-   * Note that we are return the maximum allowed integer here (not 
-   * long/unsigned, etc..). This is for easy comparisson (i.e. guarding 
+  /** @brief Max possible modified_julian_day
+   *
+   * Note that we are return the maximum allowed integer here (not
+   * long/unsigned, etc..). This is for easy comparisson (i.e. guarding
    * against overflow when comparing with ints).
    */
   constexpr static modified_julian_day max() noexcept {
     return modified_julian_day{std::numeric_limits<int>::max()};
   }
 
   /** @brief Min possible modified_julian_day
-   * 
-   * Note that we are return the minimum allowed integer here (not 
-   * long/unsigned, etc..). This is for easy comparisson (i.e. guarding 
+   *
+   * Note that we are return the minimum allowed integer here (not
+   * long/unsigned, etc..). This is for easy comparisson (i.e. guarding
    * against overflow when comparing with ints).
    */
   constexpr static modified_julian_day min() noexcept {

src/tpdate.hpp

@@ -8,8 +8,9 @@
 #ifndef __DSO_DATETIME_TWOPARTDATES2_HPP__
 #define __DSO_DATETIME_TWOPARTDATES2_HPP__
 
-#include "dtcalendar.hpp"
+#include "dtdatetime.hpp"
 #include "dtfund.hpp"
+#include <random>
 
 namespace dso {
 
@@ -276,6 +277,21 @@ class TwoPartDate {
     return TwoPartDate(51544, 86400e0 / 2e0, 'y');
   }
 
+  /** Random Date within some MJD limits
+   * @todo transfer this into a .cpp file
+   */
+  static TwoPartDate
+  random(modified_julian_day from = modified_julian_day::min(),
+         modified_julian_day to = modified_julian_day::max()) noexcept {
+    int istart = (int)from.as_underlying_type();
+    int istop = (int)to.as_underlying_type();
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_int_distribution<int> distr(istart, istop);
+    std::uniform_real_distribution<double> unif(0, 86400e0);
+    return TwoPartDate(distr(gen), unif(gen), 'y');
+  }
+
   /** Min date */
   static constexpr TwoPartDate min() noexcept {
     return TwoPartDate(datetime<nanoseconds>::min());
@@ -294,9 +310,12 @@ class TwoPartDate {
   /** Get the MJD as an intgral number, i.e. no fractional part */
   int imjd() const noexcept { return _mjd; }
 
-  /** Get the fractional seconds of the MJD */
+  /** Get the (fractional) seconds of the MJD */
   FDOUBLE seconds() const noexcept { return _fsec; }
 
+  /** Get the seconds of the MJD as fractional day. Always in range [0,1) */
+  FDOUBLE fractional_days() const noexcept { return _fsec / 86400e0; }
+
   /** @brief Get the fractional seconds of day as some multiple of seconds.
    *
    * 'Some multiple of seconds' can be any type defined with a member
@@ -495,8 +514,8 @@ class TwoPartDate {
   FDOUBLE as_mjd() const noexcept { return _fsec / SEC_PER_DAY + _mjd; }
 
   FDOUBLE jcenturies_sinceJ2000() const noexcept {
-    return (_mjd - J2000_MJD) / DAYS_IN_JULIAN_CENT +
-           _fsec / SEC_PER_DAY / DAYS_IN_JULIAN_CENT;
+    return ((static_cast<FDOUBLE>(_mjd) - J2000_MJD) + _fsec / SEC_PER_DAY) /
+           DAYS_IN_JULIAN_CENT;
   }
 
   /** @brief Convert to Julian Epoch, assuming the TT time-scale. */

src/tpdate2.hpp

@@ -0,0 +1,180 @@
+/** @file
+ * A utility class to hold datetime instances, in a continuous time-scale
+ * (e.g. TT, TAI, etc). In construst to datetime<S>, this is not a template
+ * class and uses a storage method of two numeric values (one for MJD and one
+ * for seconds of day) to represent datetime epochs.
+ */
+
+#ifndef __DSO_DATETIME_TWOPARTDATES_ALTERNATIVE_HPP__
+#define __DSO_DATETIME_TWOPARTDATES_ALTERNATIVE_HPP__
+
+#include "dtdatetime.hpp"
+#include "dtfund.hpp"
+
+namespace dso {
+
+/** A datetime class to represent epochs in any continuous system.
+ *
+ * A TwoPartDate instance conviniently splits a datetime into two numeric
+ * values:
+ * - the Modified Julian Day (i.e. an integral value), and
+ * - the time of day, which is stored in fractional seconds of day
+ *
+ * The methods of the class, invluding constructors, take special care to
+ * always keep the seconds as seconds of day, i.e. in the range [0,86400) and
+ * correspondingly increase/decrease the day count.
+ *
+ */
+class TwoPartDate2 {
+private:
+  using FDOUBLE = /*long*/ double;
+  int _mjd;      /** Mjd */
+  FDOUBLE _fday; /** fractional part of day */
+
+public:
+  /** Constructor from datetime<T> */
+#if __cplusplus >= 202002L
+  template <gconcepts::is_sec_dt T>
+#else
+  template <typename T, typename = std::enable_if_t<T::is_of_sec_type>>
+#endif
+  TwoPartDate2(const datetime<T> &d) noexcept
+      : _mjd(d.imjd().as_underlying_type()),
+        _fday(to_fractional_days<T, FDOUBLE>(d.sec())) {
+    this->normalize();
+  }
+
+  /** Constructor from a pair of doubles, such that MJD = a + b */
+  explicit TwoPartDate2(int b = 0, FDOUBLE s = 0) noexcept : _mjd(b), _fday(s) {
+    this->normalize();
+  }
+
+  /** Get the MJD as an intgral number, i.e. no fractional part */
+  int imjd() const noexcept { return _mjd; }
+
+  /** Get the fractional seconds of the MJD */
+  FDOUBLE seconds() const noexcept { return _fday * SEC_PER_DAY; }
+
+  FDOUBLE fractional_days() const noexcept { return _fday; }
+
+  /** @brief Transform the (integral part of the) date to Year Month Day */
+  ymd_date to_ymd() const noexcept { return core::mjd2ymd((long)_mjd); }
+
+  /** Add seconds to instance.
+   * @warning Does not take into account leap seconds.
+   */
+  void add_seconds(FDOUBLE sec) noexcept {
+    _fday += sec / SEC_PER_DAY;
+    this->normalize();
+  }
+  void add_seconds2(FDOUBLE sec) noexcept {
+    _fday *= SEC_PER_DAY;
+    _fday += sec;
+    _fday /= SEC_PER_DAY;
+    this->normalize();
+  }
+
+  /** Difference between two dates as integral number of days and seconds of
+   * day
+   *
+   * This is not an 'actual date' but rather a datetime interval, but can be
+   * represented by a TwoPartDate instance. If the calling instance is prior to
+   * the operand (i.e. d1-d2 with d2>d1) the interval is signed as 'negative'.
+   * This means that the number of days can be negative, but the fractional
+   * day will always be positive
+   */
+  TwoPartDate2 operator-(const TwoPartDate2 &d) const noexcept {
+    return TwoPartDate2(_mjd - d._mjd, _fday - d._fday);
+  }
+
+  /** Add two instances.
+   *
+   * In this case, the second instance can be though of a datetime interval
+   * (rather than an actual datetime instance). The right operand can be a
+   * negative interval, which means that we are going backwards in time.
+   */
+  TwoPartDate2 operator+(const TwoPartDate2 &d) const noexcept {
+    return TwoPartDate2(_mjd + d._mjd, _fday + d._fday);
+  }
+
+  /** Get the difference between two datetime instances an an arithmetic value.
+   *
+   * The difference can be obtained as a fractional days or fractional seconds,
+   * depending on the template parameter \p DT.
+   * If called as diff(d1,d2), the computation is d1-d2; the difference can be
+   * negative if d2>d1.
+   *
+   * @warning Does not take into account leap seconds.
+   */
+  template <DateTimeDifferenceType DT>
+  FDOUBLE diff(const TwoPartDate2 &d) const noexcept {
+    if constexpr (DT == DateTimeDifferenceType::FractionalDays) {
+      /* difference as fractional days */
+      return (_mjd - d._mjd) + (_fday - d._fday);
+    } else if constexpr (DT == DateTimeDifferenceType::FractionalSeconds) {
+      /* difference as fractional seconds */
+      return (_fday - d._fday) * SEC_PER_DAY + (_mjd - d._mjd) * SEC_PER_DAY;
+    } else {
+      /* difference as fractional (julian) years */
+      return this->diff<DateTimeDifferenceType::FractionalDays>(d) /
+             DAYS_IN_JULIAN_YEAR;
+    }
+  }
+
+  /** Get the date as (fractional) Julian Date */
+  FDOUBLE julian_date() const noexcept { return _fday + (_mjd + dso::MJD0_JD); }
+
+  FDOUBLE as_mjd() const noexcept { return _fday + _mjd; }
+
+  FDOUBLE jcenturies_sinceJ2000() const noexcept {
+    return (_mjd - J2000_MJD) / DAYS_IN_JULIAN_CENT +
+           _fday / DAYS_IN_JULIAN_CENT;
+  }
+
+  bool operator>(const TwoPartDate2 &d) const noexcept {
+    return (_mjd > d._mjd) || ((_mjd == d._mjd) && (_fday > d._fday));
+  }
+  bool operator>=(const TwoPartDate2 &d) const noexcept {
+    return (_mjd > d._mjd) || ((_mjd == d._mjd) && (_fday >= d._fday));
+  }
+  bool operator<(const TwoPartDate2 &d) const noexcept {
+    return (_mjd < d._mjd) || ((_mjd == d._mjd) && (_fday < d._fday));
+  }
+  bool operator<=(const TwoPartDate2 &d) const noexcept {
+    return (_mjd < d._mjd) || ((_mjd == d._mjd) && (_fday <= d._fday));
+  }
+
+  /** @brief Normalize an instance.
+   */
+  void normalize() noexcept {
+    double extra_days;
+    _fday = std::modf(_fday, &extra_days);
+    _mjd += extra_days;
+#ifdef DEBUG
+    assert(_fday >= 0e0 && _fday < 1e0);
+#endif
+    /* all done */
+    return;
+  }
+
+  TwoPartDate2 normalized() const noexcept {
+    TwoPartDate2 d(*this);
+    d.normalize();
+    return d;
+  }
+
+  /** @brief Equality operator */
+  bool operator==(const TwoPartDate2 &d) const noexcept {
+    return (_mjd == d._mjd) && (_fday == d._fday);
+  }
+
+  /** @brief In-equality operator */
+  bool operator!=(const TwoPartDate2 &d) const noexcept {
+    return !(this->operator==(d));
+  }
+
+}; /* class TwoPartDate2 */
+
+} /* namespace dso */
+
+#endif

src/twopartdates.cpp

@@ -1 +1 @@
-#include "dtcalendar.hpp"
+#include "dtdatetime.hpp"

test/precision/tp2dt.cpp

@@ -1,4 +1,4 @@
-#include "dtcalendar.hpp"
+#include "dtdatetime.hpp"
 #include "tpdate.hpp"
 #include <cassert>
 #include <random>

test/timer/test_tpdate_normalization.cpp

@@ -1,4 +1,4 @@
-#include "dtcalendar.hpp"
+#include "dtdatetime.hpp"
 #include "tpdate.hpp"
 #include <cassert>
 #include <chrono>

test/unit_tests/datetime.cpp

@@ -1,4 +1,4 @@
-#include "dtcalendar.hpp"
+#include "dtdatetime.hpp"
 #include <cassert>
 #include <random>
 

test/unit_tests/datetime_addition.cpp

@@ -1,4 +1,4 @@
-#include "dtcalendar.hpp"
+#include "dtdatetime.hpp"
 #include <cassert>
 #include <random>
 

test/unit_tests/datetime_addition_sec.cpp

@@ -1,4 +1,4 @@
-#include "dtcalendar.hpp"
+#include "dtdatetime.hpp"
 #include <cassert>
 #include <random>
 

test/unit_tests/datetime_addition_sec2.cpp

@@ -1,4 +1,4 @@
-#include "dtcalendar.hpp"
+#include "dtdatetime.hpp"
 #include <cassert>
 #include <random>
 

test/unit_tests/datetime_diff.cpp

@@ -1,4 +1,4 @@
-#include "dtcalendar.hpp"
+#include "dtdatetime.hpp"
 #include <cassert>
 #include <limits>
 #include <random>

test/unit_tests/interval.cpp

@@ -1,4 +1,4 @@
-#include "dtcalendar.hpp"
+#include "dtdatetime.hpp"
 #include <cassert>
 #include <random>
 

test/unit_tests/tpdates1.cpp

@@ -1,4 +1,4 @@
-#include "dtcalendar.hpp"
+#include "dtdatetime.hpp"
 #include "tpdate.hpp"
 #include <cassert>
 #include <random>