clean me

src/dtfund.hpp

@@ -148,10 +148,13 @@ inline constexpr double jd2epj(double dj1, double dj2) noexcept {
  * @param[in] mjd1 Second part of MJD (if any), such that MJD = mjd0 + mjd1
  * @return The input date as Julian Epoch.
  */
-inline constexpr double mjd2epj(double mjd0, double mjd1 = 0e0) noexcept {
+inline constexpr double mjd2epj1(double mjd0, double mjd1 = 0e0) noexcept {
   return 2000e0 + ((mjd0 - J2000_MJD) / DAYS_IN_JULIAN_YEAR +
                    mjd1 / DAYS_IN_JULIAN_YEAR);
 }
+inline constexpr double mjd2epj2(double mjd0, double mjd1 = 0e0) noexcept {
+  return 2000e0 + ((mjd0 - J2000_MJD) + mjd1) / DAYS_IN_JULIAN_YEAR;
+}
 
 /** @brief Julian Epoch to Modified Julian Date
  *
@@ -163,6 +166,31 @@ inline constexpr double mjd2epj(double mjd0, double mjd1 = 0e0) noexcept {
 inline constexpr double epj2mjd(double epj) noexcept {
   return J2000_MJD + (epj - 2000e0) * DAYS_IN_JULIAN_YEAR;
 }
+
+/** @brief Julian Epoch to two-part Modified Julian Date
+ *
+ * This function returns the correponding MJD the input Julian Epoch as a 
+ * two-part MJD, where the first, "big" part is always J2000_MJD and the 
+ * rest ("small" part) is returned in the parameter \p mjd1. So that the 
+ * actual MJD = BigPart (i.e. J200_MJD) + SmallPart (i.e. mjd1).
+ * This is meant to better preserve accuracy and convieniently place the 
+ * result in a TwoPartDate.
+ *
+ * @param[in] epj The Julian Epoch to convert
+ * @param[out] mjd1 Small part of MJD
+ * @return Big part of MJD, i.e. J2000_MJD such the the MJD corresponding to 
+ *         the input Julian Epoch is MJD = BigPart + SmallPart
+ *
+ * @see IAU SOFA epj2jd
+ */
+inline double epj2mjd(double epj, double &mjd1) noexcept {
+  /* lose the .5 part */
+  const double ipart = static_cast<int>(J2000_MJD);
+  mjd1 = (epj - 2000e0) * DAYS_IN_JULIAN_YEAR + .5e0;
+  double iextra;
+  mjd1 = std::modf(mjd1, &iextra);
+  return ipart+iextra;
+}
 } /* namespace core */
 
 /** @brief For a given UTC date, calculate delta(AT) = TAI-UTC.

src/tpdate.hpp

@@ -9,6 +9,7 @@
 #define __DSO_DATETIME_TWOPARTDATES2_HPP__
 
 #include "dtcalendar.hpp"
+#include "dtfund.hpp"
 
 namespace dso {
 
@@ -470,6 +471,18 @@ class TwoPartDate {
            _fsec / SEC_PER_DAY / DAYS_IN_JULIAN_CENT;
   }
 
+  FDOUBLE epj1() const noexcept {
+    // return core::mjd2epj((double)imjd(), seconds()/SEC_PER_DAY);
+    const auto d = this->operator-(TwoPartDate(51544, 86400e0/2));
+    return (d.imjd() + d.seconds()/SEC_PER_DAY)/DAYS_IN_JULIAN_YEAR + 2e3;
+  }
+  FDOUBLE epj2() const noexcept {
+    return core::mjd2epj1((double)imjd(), seconds()/SEC_PER_DAY);
+  }
+  FDOUBLE epj3() const noexcept {
+    return core::mjd2epj2((double)imjd(), seconds()/SEC_PER_DAY);
+  }
+
   bool operator>(const TwoPartDate &d) const noexcept {
     return (_mjd > d._mjd) || ((_mjd == d._mjd) && (_fsec > d._fsec));
   }
@@ -536,6 +549,28 @@ class TwoPartDate {
 
 }; /* class TwoPartDate */
 
+/** @brief Julian Epoch to two-part Modified Julian Date
+ *
+ * @param[in] epj The Julian Epoch to convert
+ * @return The corresponding MJD as a TwoPartDate
+ */
+inline TwoPartDate epjtp1(double epj) noexcept {
+  double smallpart;
+  const double bigpart = core::epj2mjd(epj, smallpart);
+  return TwoPartDate(bigpart, smallpart*SEC_PER_DAY);
+}
+inline TwoPartDate epjtp2(double epj) noexcept {
+  // return J2000_MJD + (epj - 2000e0) * DAYS_IN_JULIAN_YEAR;
+  double epjBig;
+  const double epjSmall = std::modf(epj, &epjBig) * DAYS_IN_JULIAN_YEAR;
+  epjBig = (epjBig - 2000e0) * DAYS_IN_JULIAN_YEAR + J2000_MJD;
+  /* Now MJD is BigPart + SmallPart */
+  double i1,i2,i3;
+  double d = std::modf(epjBig, &i1) + std::modf(epjSmall, &i2);
+  d = std::modf(d,&i3);
+  return TwoPartDate(i1+i2+i3, d*SEC_PER_DAY);
+}
+
 } /* namespace dso */
 
 #endif

test/sofa/epj_date.cpp

@@ -0,0 +1,116 @@
+#include "dtfund.hpp"
+#include "sofa.h"
+#include "tpdate.hpp"
+#include <random>
+
+using namespace dso;
+
+/* number of tests to perform (pre template parameter) */
+long num_tests = 1'000'000;
+
+int main() {
+  /* Generators for random numbers ... */
+  std::random_device rd;
+  std::mt19937 gen(rd());
+  std::uniform_int_distribution<> ydstr(1972, 2050); /* range for years */
+  std::uniform_int_distribution<> mdstr(1, 12);      /* range for months */
+  std::uniform_int_distribution<> ddstr(1, 31); /* range for day of month */
+  std::uniform_int_distribution<> msstr(
+      0, 86400 * 1000); /* range for milliseconds of day */
+
+  double maxdiffs[10], avediffs[10];
+  for (int i=0; i<10; i++) {
+    maxdiffs[i] = 0e0;
+    avediffs[i] = 0e0;
+  }
+  int n=0;
+  for (long i = 0; i < num_tests; i++) {
+    const int iy = ydstr(gen);
+    const int im = mdstr(gen);
+    const int id = ddstr(gen);
+    const int ms = msstr(gen);
+    try {
+      long imjd = core::cal2mjd(iy, im, id);
+      const TwoPartDate d(imjd, (ms / 1e3));
+
+      {
+        const double je = d.epj1();
+        TwoPartDate di = epjtp1(je);
+        //printf("DeltaMjd = %+d DeltaSec=%+.6f[sec]\n", d.imjd() - di.imjd(),
+        //       d.seconds() - di.seconds());
+        assert(d.imjd() - di.imjd() == 0);
+        if (std::abs(d.seconds()-di.seconds())>maxdiffs[0]) maxdiffs[0] = std::abs(d.seconds()-di.seconds());
+        avediffs[0] += std::abs(d.seconds()-di.seconds());
+      }
+      {
+        const double je = d.epj2();
+        TwoPartDate di = epjtp1(je);
+        //printf("DeltaMjd = %+d DeltaSec=%+.6f[sec]\n", d.imjd() - di.imjd(),
+        //       d.seconds() - di.seconds());
+        assert(d.imjd() - di.imjd() == 0);
+        if (std::abs(d.seconds()-di.seconds())>maxdiffs[1]) maxdiffs[1] = std::abs(d.seconds()-di.seconds());
+        avediffs[1] += std::abs(d.seconds()-di.seconds());
+      }
+      {
+        const double je = d.epj3();
+        TwoPartDate di = epjtp1(je);
+        //printf("DeltaMjd = %+d DeltaSec=%+.6f[sec]\n", d.imjd() - di.imjd(),
+        //       d.seconds() - di.seconds());
+        assert(d.imjd() - di.imjd() == 0);
+        if (std::abs(d.seconds()-di.seconds())>maxdiffs[2]) maxdiffs[2] = std::abs(d.seconds()-di.seconds());
+        avediffs[2] += std::abs(d.seconds()-di.seconds());
+      }
+      {
+        const double je = d.epj1();
+        TwoPartDate di = epjtp2(je);
+        //printf("DeltaMjd = %+d DeltaSec=%+.6f[sec]\n", d.imjd() - di.imjd(),
+        //       d.seconds() - di.seconds());
+        assert(d.imjd() - di.imjd() == 0);
+        if (std::abs(d.seconds()-di.seconds())>maxdiffs[3]) maxdiffs[3] = std::abs(d.seconds()-di.seconds());
+        avediffs[3] += std::abs(d.seconds()-di.seconds());
+      }
+      {
+        const double je = d.epj2();
+        TwoPartDate di = epjtp2(je);
+        //printf("DeltaMjd = %+d DeltaSec=%+.6f[sec]\n", d.imjd() - di.imjd(),
+        //       d.seconds() - di.seconds());
+        assert(d.imjd() - di.imjd() == 0);
+        if (std::abs(d.seconds()-di.seconds())>maxdiffs[4]) maxdiffs[4] = std::abs(d.seconds()-di.seconds());
+        avediffs[4] += std::abs(d.seconds()-di.seconds());
+      }
+      {
+        const double je = d.epj3();
+        TwoPartDate di = epjtp2(je);
+        //printf("DeltaMjd = %+d DeltaSec=%+.6f[sec]\n", d.imjd() - di.imjd(),
+        //       d.seconds() - di.seconds());
+        assert(d.imjd() - di.imjd() == 0);
+        if (std::abs(d.seconds()-di.seconds())>maxdiffs[5]) maxdiffs[5] = std::abs(d.seconds()-di.seconds());
+        avediffs[5] += std::abs(d.seconds()-di.seconds());
+      }
+
+      /* do the same with SOFA */
+      {
+        const double jd1 = d.seconds() / SEC_PER_DAY;
+        const double jd2 = d.imjd() + MJD0_JD;
+        const double sje = iauEpj(jd2, jd1);
+        double jd11, jd21;
+        iauEpj2jd(sje, &jd11, &jd21);
+        //printf("Sofa DeltaSec=%+.6f[sec]\n", (jd21 - d.as_mjd()) * 86400e0);
+        const double dsec = ((jd2-jd11) + (jd1-jd21)) * SEC_PER_DAY;
+        if (std::abs(dsec) > maxdiffs[6]) maxdiffs[6] = std::abs(dsec);
+        avediffs[6] += std::abs(dsec);
+      }
+      ++i;
+    } catch (std::exception &) {
+      ;
+    }
+    if (i % 10)
+      printf("%ld/%ld\r", i, num_tests);
+  }
+
+  for (int i=0; i<6; i++) {
+    printf("[%d]: MaxDiff=%+.6f[sec] AveDiff=%.6f[sec]\n", i, maxdiffs[i], avediffs[i]/n);
+  }
+
+  return 0;
+}