epj ok

README.md

@@ -102,3 +102,9 @@ by running the test program (unit_tests/tpdates3.cpp)[unit_tests/tpdates3.cpp]
 Transforming from a `datetime<seconds>` instance to an instance of type 
 `TwoPartDate` preserves a precision better than 2e-12 seconds. This is verified 
 by running the test program (unit_tests/tpdates4.cpp)[unit_tests/tpdates4.cpp]
+
+**Transforming from/to Julian Epoch**
+
+Transforming from a `TwoPartDate` to a *Julian Epoch* and back, preserves an 
+accuracy of `~1e-5` seconds. This is verified by running the test program 
+(test/sofa/epj_date.cpp)[test/sofa/epj_date.cpp].

src/dtfund.hpp

@@ -113,6 +113,8 @@ inline constexpr long ydoy2mjd(long iyr, long idoy) {
 }
 
 /* @brief Julian Date to Julian Epoch
+ *
+ * The function assumes the TT time-scale (for input and output dates).
  *
  * Julian epoch uses the Julian year of exactly 365.25 days, and the TT time
  * scale; Julian epoch 2000.0 is defined to be 2000 January 1.5, which is
@@ -136,27 +138,23 @@ inline constexpr double jd2epj(double dj1, double dj2) noexcept {
 }
 
 /** @brief Modified Julian Date to Julian Epoch
+ * The function assumes the TT time-scale (for input and output dates).
  *
  * Convert a Modified Julian date to a Julian Epoch. The MJD can be given as a
  * single value (i.e. in parameter \p mjd0) or can be split into two parts,
  * (e.g. the first being the intergal part of MJD and the second be fractional
  * day).
  *
- * @see jd2epj
- *
  * @param[in] mjd0 The Modified Julian Date (first part or whole number)
  * @param[in] mjd1 Second part of MJD (if any), such that MJD = mjd0 + mjd1
  * @return The input date as Julian Epoch.
  */
-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 {
+inline constexpr double mjd2epj(double mjd0, double mjd1 = 0e0) noexcept {
   return 2000e0 + ((mjd0 - J2000_MJD) + mjd1) / DAYS_IN_JULIAN_YEAR;
 }
 
 /** @brief Julian Epoch to Modified Julian Date
+ * The function assumes the TT time-scale (for input and output dates).
  *
  * @param[in] epj The Julian Epoch to convert
  * @return The corresponding (fractional) Modified Julian Date
@@ -170,25 +168,26 @@ inline constexpr double epj2mjd(double epj) noexcept {
 /** @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).
+ * two-part MJD, where the first, "big" part is the MJDay and the 
+ * rest ("small" part) is returned in the parameter \p fday, representing the 
+ * fractional part of the MJday.  
+ * So, the actual MJD = BigPart (i.e. Day) + SmallPart (i.e. fraction of day).
  * This is meant to better preserve accuracy and convieniently place the 
  * result in a TwoPartDate.
+ * The function assumes the TT time-scale (for input and output dates).
  *
  * @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
+ * @param[out] fday Fractional part of MJDay
+ * @return Integral part of MJD, i.e. the MJDay, such the the MJD 
+ *         corresponding to the input Julian Epoch is 
+ *         MJD = BigPart + SmallPart
  */
-inline double epj2mjd(double epj, double &mjd1) noexcept {
+inline double epj2mjd(double epj, double &fday) noexcept {
   /* lose the .5 part */
   const double ipart = static_cast<int>(J2000_MJD);
-  mjd1 = (epj - 2000e0) * DAYS_IN_JULIAN_YEAR + .5e0;
+  fday = (epj - 2000e0) * DAYS_IN_JULIAN_YEAR + .5e0;
   double iextra;
-  mjd1 = std::modf(mjd1, &iextra);
+  fday = std::modf(fday, &iextra);
   return ipart+iextra;
 }
 } /* namespace core */

src/tpdate.hpp

@@ -471,16 +471,9 @@ 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);
+  /** @brief Convert to Julian Epoch, assuming the TT time-scale. */
+  FDOUBLE epj() const noexcept {
+    return core::mjd2epj((double)imjd(), seconds()/SEC_PER_DAY);
   }
 
   bool operator>(const TwoPartDate &d) const noexcept {
@@ -549,28 +542,18 @@ class TwoPartDate {
 
 }; /* class TwoPartDate */
 
-/** @brief Julian Epoch to two-part Modified Julian Date
+/** @brief Julian Epoch to two-part Modified Julian Date (TT).
+ *
+ * The function assume the TT time-scale.
  *
  * @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 epj2tpd(double epj) noexcept {
+  double fday;
+  const double mjd = core::epj2mjd(epj, fday);
+  return TwoPartDate(mjd, fday*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

@@ -5,6 +5,14 @@
 
 using namespace dso;
 
+constexpr const double EPSILON = 1e-12; /* fractional days */
+template <typename T>
+inline bool fequal(const T &a, const T &b, double eps = EPSILON) noexcept {
+  if (std::abs(a - b) <= eps)
+    return true;
+  return std::abs(a - b) <= (eps * std::max(std::abs(a), std::abs(b)));
+}
+
 /* number of tests to perform (pre template parameter) */
 long num_tests = 1'000'000;
 
@@ -18,12 +26,14 @@ int main() {
   std::uniform_int_distribution<> msstr(
       0, 86400 * 1000); /* range for milliseconds of day */
 
+#ifdef DEBUG
   double maxdiffs[10], avediffs[10];
   for (int i=0; i<10; i++) {
     maxdiffs[i] = 0e0;
     avediffs[i] = 0e0;
   }
   int n=0;
+#endif
   for (long i = 0; i < num_tests; i++) {
     const int iy = ydstr(gen);
     const int im = mdstr(gen);
@@ -33,84 +43,55 @@ int main() {
       long imjd = core::cal2mjd(iy, im, id);
       const TwoPartDate d(imjd, (ms / 1e3));
 
+      /* TwoPartDate to Julian Epoch and back */
+      double epj_lib;
       {
-        const double je = d.epj1();
-        TwoPartDate di = epjtp1(je);
-        //printf("DeltaMjd = %+d DeltaSec=%+.6f[sec]\n", d.imjd() - di.imjd(),
-        //       d.seconds() - di.seconds());
+        const double je = d.epj();
+        TwoPartDate di = epj2tpd(je);
         assert(d.imjd() - di.imjd() == 0);
+#ifdef DEBUG
         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());
+#endif
+        epj_lib = je;
+
+        /* compare initial to resulting dates */
+        assert(fequal(d.seconds(),di.seconds(),1e-5));
       }
 
       /* do the same with SOFA */
+      double epj_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);
+#ifdef DEBUG
         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);
+        if (std::abs(dsec) > maxdiffs[1]) maxdiffs[1] = std::abs(dsec);
+        avediffs[1] += std::abs(dsec);
+#endif
+        epj_sofa = sje;
       }
-      ++i;
+
+      /* compare SOFA's epj to lib's epj */
+      assert(fequal(epj_lib, epj_sofa));
+#ifdef DEBUG
+      ++n;
+#endif
     } catch (std::exception &) {
       ;
     }
     if (i % 10)
       printf("%ld/%ld\r", i, num_tests);
   }
 
-  for (int i=0; i<6; i++) {
+#ifdef DEBUG
+  for (int i=0; i<2; i++) {
     printf("[%d]: MaxDiff=%+.6f[sec] AveDiff=%.6f[sec]\n", i, maxdiffs[i], avediffs[i]/n);
   }
+#endif
 
   return 0;
 }