still not working

DSOlab · Oct 31, 2023 · 19b9703 · 19b9703
1 parent bf93b0e
commit 19b9703
Show file tree

Hide file tree

Showing 4 changed files with 118 additions and 29 deletions.
diff --git a/src/datetime_tops.hpp b/src/datetime_tops.hpp
@@ -391,7 +391,6 @@ template <typename S, typename T = double,
           typename = std::enable_if_t<std::is_floating_point<T>::value>>
 #endif
 T to_fractional_days(S nsec) noexcept {
-  //printf("\tto_fractional_days got %ld\n", nsec.as_underlying_type());
   const T sec = static_cast<T>(nsec.__member_ref__());
   //return sec / static_cast<T>(S::max_in_day);
   return sec / S::max_in_day;

diff --git a/src/tpdate.hpp b/src/tpdate.hpp
@@ -66,7 +66,7 @@ TwoPartJulianDate jd_split(double mjd, double fday) noexcept {
  */
 class TwoPartDate {
 private:
-  using FDOUBLE = long double;
+  using FDOUBLE = /*long*/ double;
   double _mjd;  /** Mjd */
   FDOUBLE _fday; /** fractional days */
 
@@ -132,21 +132,24 @@ class TwoPartDate {
    * @warning Does not take into account leap seconds.
    */
   void add_seconds(FDOUBLE sec) noexcept {
-    // first approach
-    // _fday += sec / SEC_PER_DAY;
-    // this->normalize();
-    //
-    // second approach
     FDOUBLE dsec = _fday * SEC_PER_DAY;
     dsec += sec;
     _fday = dsec / SEC_PER_DAY;
     this->normalize();
-    //constexpr const long double SPD = 86400e0;
-    //long double dsec = _fday * SPD;
-    //dsec += sec;
-    //long double fday = dsec / SPD;
-    //_fday = fday;
-    //this->normalize();
+  }
+
+  /** Add seconds to instance.
+   * @warning Does not take into account leap seconds.
+   */
+  void add_seconds(FDOUBLE fsec, FDOUBLE &err) noexcept {
+    /* result: 23:59:59.000'000'008 */
+    const FDOUBLE sec = fsec + err;
+    const FDOUBLE dsec = _fday * SEC_PER_DAY;
+    const FDOUBLE s = sec + dsec;
+    const FDOUBLE z = s - dsec;
+    err = sec - z;
+    _fday = s / SEC_PER_DAY;
+    this->normalize();
   }
 
   /** Difference between two dates as integral number of days and seconds of
@@ -296,23 +299,19 @@ class TwoPartDate {
    * stored in the _fday part, while _mjd holds the integral part of date
    */
   void normalize() noexcept {
-    /* fractional part should NOT be negative */
-    while (_fday < 0e0) {
-      _fday = 1 - _fday;
+    /* split _mjd to fractional and integral part */
+    FDOUBLE intpart;
+    FDOUBLE fraction = std::modf(_mjd, &intpart);
+    _mjd = intpart;
+    /* split _fday to fractional and integral part */
+    fraction += std::modf(_fday, &intpart);
+    _mjd += intpart;
+    _fday = fraction;
+    /* fraction cannot be negative */
+    if (_fday < 0e0) {
+      _fday += 1e0;
       _mjd -= 1e0;
     }
-    FDOUBLE fmore(0e0), extra(0e0);
-    /* check if _mjd part has a fractional part */
-    if ((fmore = std::modf((FDOUBLE)_mjd, &extra)) != 0e0) {
-      /* assign fractional part to _fday and keep integral part to _mjd */
-      _fday += fmore;
-      _mjd = extra;
-    }
-    /* check if fractional part is >= 1e0 */
-    if (_fday >= 1e0) {
-      _fday = std::modf(_fday, &extra);
-      _mjd += extra;
-    }
 #ifdef DEBUG
     assert(_fday >= 0e0 && _fday < 1e0);
 #endif

diff --git a/test/unit_tests/dwrite9.cpp b/test/unit_tests/dwrite9.cpp
@@ -225,10 +225,13 @@ int main() {
 
   d1 = datetime<nanoseconds>(year(2023), month(10), day_of_month(24),
                              nanoseconds(0));
+  [[maybe_unused]]double secerr = 0e0;
   td1 = TwoPartDate(d1);
   for (long i = 0; i < 86400L-1; i++)
-    td1.add_seconds(1e0);
+    td1.add_seconds(1e0,secerr);
+  printf("Seconds reached: %.15e", td1.sec_of_day<seconds>());
   to_char<YMDFormat::YYYYMMDD, HMSFormat::HHMMSSF>(td1, buffer);
+  printf(", i.e. %.29s\n", buffer);
   assert(!std::strncmp(buffer, "2023/10/24 23:59:59.000000000", sz));
   td1.add_seconds(.999999999e0);
   to_char<YMDFormat::YYYYMMDD, HMSFormat::HHMMSSF>(td1, buffer);

diff --git a/test/unit_tests/tpdate_add.cpp b/test/unit_tests/tpdate_add.cpp
@@ -0,0 +1,88 @@
+#include "datetime_write.hpp"
+#include <cassert>
+#include <cstring>
+
+using namespace dso;
+constexpr const long SEC = nanoseconds::sec_factor<long>();
+
+int main() {
+
+  /* check dates */
+  datetime<nanoseconds> d1(year(2023), month(10), day_of_month(24),
+                            nanoseconds(0));
+  TwoPartDate td1(d1);
+  printf("Central date is                   : %.5f %.15e\n", td1.imjd(), td1.fday());
+
+  for (int i=0; i<86400; i++) {
+    td1.add_seconds(1e0);
+  }
+  printf("After adding   one day in seconds : %.5f %.15e\n", td1.imjd(), td1.fday());
+
+  td1 = datetime<nanoseconds>(year(2023), month(10), day_of_month(24), nanoseconds(0));
+  for (int i=0; i<2*86400; i++) {
+    td1.add_seconds(1e0);
+    //if (i%3600==0) printf("\t%.5f %.15e\n", td1.imjd(), td1.fday());
+  }
+  printf("After adding   two days in seconds: %.5f %.15e\n", td1.imjd(), td1.fday());
+
+  td1 = datetime<nanoseconds>(year(2023), month(10), day_of_month(24), nanoseconds(0));
+  for (int i=0; i<86400; i++) {
+    td1.add_seconds(-1e0);
+  }
+  printf("After removing one day in seconds : %.5f %.15e\n", td1.imjd(), td1.fday());
+
+  td1 = datetime<nanoseconds>(year(2023), month(10), day_of_month(24), nanoseconds(0));
+  for (int i=0; i<2*86400; i++) {
+    td1.add_seconds(-1e0);
+  }
+  printf("After removing two days in seconds: %.5f %.15e\n", td1.imjd(), td1.fday());
+
+  td1 = datetime<nanoseconds>(year(2023), month(10), day_of_month(24), nanoseconds(0));
+  for (int i=0; i<86400; i++) {
+    td1.add_seconds(2e0);
+    td1.add_seconds(-1e0);
+  }
+  printf("After adding   one day in seconds : %.5f %.15e\n", td1.imjd(), td1.fday());
+
+
+  double error = 0e0;
+  td1 = datetime<nanoseconds>(year(2023), month(10), day_of_month(24), nanoseconds(0));
+  printf("Central date is                   : %.5f %.15e\n", td1.imjd(), td1.fday());
+
+  for (int i=0; i<86400; i++) {
+    td1.add_seconds(1e0,error);
+  }
+  printf("After adding   one day in seconds : %.5f %.15e\n", td1.imjd(), td1.fday());
+
+  error = 0e0;
+  td1 = datetime<nanoseconds>(year(2023), month(10), day_of_month(24), nanoseconds(0));
+  for (int i=0; i<2*86400; i++) {
+    td1.add_seconds(1e0,error);
+    //if (i%3600==0) printf("\t%.5f %.15e\n", td1.imjd(), td1.fday());
+  }
+  printf("After adding   two days in seconds: %.5f %.15e\n", td1.imjd(), td1.fday());
+
+  error = 0e0;
+  td1 = datetime<nanoseconds>(year(2023), month(10), day_of_month(24), nanoseconds(0));
+  for (int i=0; i<86400; i++) {
+    td1.add_seconds(-1e0,error);
+  }
+  printf("After removing one day in seconds : %.5f %.15e\n", td1.imjd(), td1.fday());
+
+  error = 0e0;
+  td1 = datetime<nanoseconds>(year(2023), month(10), day_of_month(24), nanoseconds(0));
+  for (int i=0; i<2*86400; i++) {
+    td1.add_seconds(-1e0,error);
+  }
+  printf("After removing two days in seconds: %.5f %.15e\n", td1.imjd(), td1.fday());
+
+  error = 0e0;
+  td1 = datetime<nanoseconds>(year(2023), month(10), day_of_month(24), nanoseconds(0));
+  for (int i=0; i<86400; i++) {
+    td1.add_seconds(2e0,error);
+    td1.add_seconds(-1e0,error);
+  }
+  printf("After adding   one day in seconds : %.5f %.15e\n", td1.imjd(), td1.fday());
+
+  return 0;
+};