diff --git a/clang-format b/.clang-format similarity index 100% rename from clang-format rename to .clang-format diff --git a/CMakeLists.txt b/CMakeLists.txt index 748d6f81..5297f584 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -228,6 +228,7 @@ endif() if(kep3_BUILD_BENCHMARKS) add_subdirectory(benchmark) + endif() if(kep3_BUILD_PYTHON_BINDINGS) @@ -247,3 +248,4 @@ if(kep3_BUILD_PYTHON_BINDINGS) # Build directory add_subdirectory("${CMAKE_SOURCE_DIR}/pykep") endif() + diff --git a/include/kep3/core_astro/convert_julian_dates.hpp b/include/kep3/core_astro/convert_julian_dates.hpp index b26c509e..946dc6d3 100644 --- a/include/kep3/core_astro/convert_julian_dates.hpp +++ b/include/kep3/core_astro/convert_julian_dates.hpp @@ -7,16 +7,21 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. -#ifndef kep3_CONVERT_JULIAN_DATES_H -#define kep3_CONVERT_JULIAN_DATES_H +#ifndef kep3_CONVERT_JULIAN_DATES_HPP +#define kep3_CONVERT_JULIAN_DATES_HPP + +#include "kep3/epoch.hpp" +#include +namespace kep3 +{ + inline double jd2mjd(double in) { return (in - 2400000.5); } + inline double jd2mjd2000(double in) { return (in - 2451544.5); } + inline double mjd2jd(double in) { return (in + 2400000.5); } + inline double mjd2mjd2000(double in) { return (in - 51544); } + inline double mjd20002jd(double in) { return (in + 2451544.5); } + inline double mjd20002mjd(double in){ return (in + 51544); } + -namespace kep3 { -inline double jd2mjd(double in) { return (in - 2400000.5); } -inline double jd2mjd2000(double in) { return (in - 2451544.5); } -inline double mjd2jd(double in) { return (in + 2400000.5); } -inline double mjd2mjd2000(double in) { return (in - 51544); } -inline double mjd20002jd(double in) { return (in + 2451544.5); } -inline double mjd20002mjd(double in) { return (in + 51544); } } // namespace kep3 -#endif // kep3_CONVERT_JULIAN_DATES_H +#endif // kep3_CONVERT_JULIAN_DATES_HPP diff --git a/include/kep3/detail/s11n.hpp b/include/kep3/detail/s11n.hpp index 5608b07e..6878cc30 100644 --- a/include/kep3/detail/s11n.hpp +++ b/include/kep3/detail/s11n.hpp @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include diff --git a/include/kep3/epoch.hpp b/include/kep3/epoch.hpp index 093349cd..eaec0fe5 100644 --- a/include/kep3/epoch.hpp +++ b/include/kep3/epoch.hpp @@ -7,100 +7,240 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. -#ifndef kep3_EPOCH_H -#define kep3_EPOCH_H - -#include +#ifndef kep3_EPOCH_HPP +#define kep3_EPOCH_HPP +#include +#include +#include #include - -#include -#include -#include +#include #include #include +#include +#include +#include /// Keplerian Toolbox /** * This namespace contains astrodynamics and space flight mechanics routines * that are related to keplerian motions or models. */ + namespace kep3 { +using namespace std::literals; +namespace chr = std::chrono; +using uint = unsigned int; + +template struct is_duration : std::false_type {}; + +template +struct is_duration> : std::true_type {}; + +template +using enable_if_is_duration = std::enable_if_t::value>; + +struct kep_clock : public chr::system_clock { + + /** + * @brief Custom clock. + * Used for constructing epochs with a custom reference point (1 Jan 2000). + * By defining a custom clock, we avoid the overflow + * that std::chrono::system_clock suffers at +/- 292 years. + * To do that, we lower the resolution to 1 us (microsecond), + * which gives the clock a range of +/- 292 thousand years. + * + * NOTE: Adding durations of less than 1 us to an epoch (defined below) + * would not be registered. + * + * NOTE: As of C++20, the standard guarantees that std::chrono::system_clock + * uses the UNIX time reference point, which is midnight on 1 January 1970 + * (1970-01-01T00:00:00). We correct for that here in order to bring the + * reference point forward to midnight on 1 January 2000 + * (2000-01-01T00:00:00), which is 0 MJD2000. + */ + using rep = int_fast64_t; + // Resolution of (1 / 1'000'000)s = 1 us + using period = std::ratio<1, 1'000'000>; + using duration = chr::duration; + using time_point = chr::time_point; + static constexpr bool is_steady = false; + // Number of seconds from midnight on 1 Jan 1970 to midnight on 1 Jan 2000. + static constexpr chr::seconds y2k_offset{946684800s}; + + static constexpr time_point ref_epoch{kep_clock::time_point{} + y2k_offset}; + + static constexpr std::time_t to_time_t(const time_point &t) noexcept { + return static_cast( + chr::duration_cast(t.time_since_epoch() + y2k_offset) + .count()); + } + + static constexpr time_point from_time_t(std::time_t t) noexcept { + return chr::time_point_cast( + time_point(chr::seconds(t) - y2k_offset)); + } +}; /// epoch class. /** - * This class defines and contains a non-gregorian date (i.e. a date expressed - * in julian form). It also provides the user with an interface to boost - * gregorian dates (see boost documentation at - * http://www.boost.org/doc/libs/1_42_0/doc/html/date_time.html) - * using the posix time. - * The date is defined in MJD2000 (double) as a - * private member - * - * @author Dario Izzo (dario.izzo _AT_ googlemail.com) + * This class defines and contains a non-Gregorian date (i.e., a date expressed + * in Julian format). The date is defined in MJD2000 format as a + * kep_clock::time_point private member. Types of non-Gregorian dates supported: + * - Julian Date (JD): the number of days passed since January 1, 4713 BC + * at 12:00 (noon). + * - Modified Julian Date (MJD): the number of days passed since November + * 17, 1858 at 00:00 (midnight). + * - Modified Julian Date 2000 (MJD2000): the number of days passed since + * Juanuary 1, 2000 at 00:00 (midnight). */ class kep3_DLL_PUBLIC epoch { public: - /** Types of non gregorian dates supported. Julian Date (JD) is the number of - * days passed since January 1, 4713 BC at noon. Modified Julian Date (MJD) is - * the number of days passed since November 17, 1858 at 00:00 am. The Modified - * Julian Date 2000 (MJD2000) is the number of days passed since Juanuary 1, - * 2000 at 00:00am. - */ - enum julian_type { MJD2000, MJD, JD }; + enum class julian_type { MJD2000, MJD, JD }; /** Constructors */ - explicit epoch(const double & = 0., julian_type = MJD2000); - epoch(const boost::gregorian::greg_day &day, - const boost::gregorian::greg_month &month, - const boost::gregorian::greg_year &year); - explicit epoch(const boost::posix_time::ptime &posix_time); - - /** Computing non-gregorian dates */ - [[nodiscard]] double mjd2000() const; - [[nodiscard]] double jd() const; - [[nodiscard]] double mjd() const; - - /** Interface to boost::posix_time::ptime */ - [[nodiscard]] boost::posix_time::ptime get_posix_time() const; - void set_posix_time(const boost::posix_time::ptime &); - - /** operators overloads for sum diff (epoch-days) and the comparison operators + // Default constructor + epoch(); + + // Constructor for days (as a floating-point value) + explicit epoch(double epoch_in, + julian_type epoch_type = julian_type::MJD2000); + + /** + * Constructs an epoch from a std::chrono::duration. + * The reference point is assumed to be MJD2000. + * \param[in] time The time as a duration. + */ + template > + explicit epoch(const Duration &duration) + : tp{kep_clock::time_point{} + duration} {} + + // Constructor from duration&&) + template > + explicit epoch(Duration &&duration) + : tp{kep_clock::time_point{} + std::forward(duration)} {} + + + // Constructor for datetime broken down into its constituents. + explicit epoch(const int y, const uint mon, const uint d, const int h = 0, + const int min = 0, const int s = 0, const int ms = 0, const int us = 0); + + /* Computing non-Gregorian dates */ + + /** + * @return Number of days since 0 JD (including fractional days). + */ + [[nodiscard]] constexpr double jd() const { + return chr::duration>(tp.time_since_epoch() - kep_clock::y2k_offset + 211813444800s).count(); + } + + /** + * @return Number of days since 0 MJD (including fractional days). + */ + [[nodiscard]] constexpr double mjd() const { + return chr::duration>(tp.time_since_epoch() - kep_clock::y2k_offset + 4453401600s).count(); + } + + /** + * @return Number of days since 0 MJD2000 (including fractional days). + */ + [[nodiscard]] constexpr double mjd2000() const { + return chr::duration>(tp.time_since_epoch() - kep_clock::y2k_offset).count(); + } + + /* Helper functions for constructors */ + static kep_clock::time_point make_tp(const int y, const uint mon, const uint d, const int h = 0, + const int min = 0, const int s = 0, const int ms = 0, const int us = 0); + + static kep_clock::time_point make_tp(double epoch_in, julian_type epoch_type); + + // Conversions + static constexpr kep_clock::time_point tp_from_days(double days); + + // Duration conversions + static constexpr double as_sec(kep_clock::duration &&d) { + return std::chrono::duration(d) + .count(); + } + + // Printing + std::string as_utc_string() const; + static std::string as_utc_string(const kep_clock::time_point& tp); + + /** operator overloads for sum and diff (epoch-days) and comparison + * operators * **/ - epoch &operator+=(double rhs); - epoch &operator-=(double rhs); + + kep3_DLL_PUBLIC friend std::ostream &operator<<(std::ostream &s, + epoch const &epoch_in); + + template > + epoch &operator+=(const Duration &duration) { + tp += chr::duration_cast(duration); + return *this; + } + + template > + epoch &operator-=(const Duration &duration) { + tp -= chr::duration_cast(duration); + return *this; + } + kep3_DLL_PUBLIC friend bool operator>(const epoch &c1, const epoch &c2); kep3_DLL_PUBLIC friend bool operator<(const epoch &c1, const epoch &c2); kep3_DLL_PUBLIC friend bool operator>=(const epoch &c1, const epoch &c2); kep3_DLL_PUBLIC friend bool operator<=(const epoch &c1, const epoch &c2); kep3_DLL_PUBLIC friend bool operator==(const epoch &c1, const epoch &c2); kep3_DLL_PUBLIC friend bool operator!=(const epoch &c1, const epoch &c2); - kep3_DLL_PUBLIC friend epoch operator+(epoch lhs, double rhs); - kep3_DLL_PUBLIC friend epoch operator-(epoch lhs, double rhs); - kep3_DLL_PUBLIC friend double operator-(const epoch &lhs, const epoch &rhs); -private: + template > + epoch operator+(const Duration &duration) { + return epoch(tp + chr::duration_cast(duration)); + } + + template > + epoch operator-(const Duration &duration) { + return epoch(tp - chr::duration_cast(duration)); + } + + static constexpr auto days(const double value) { + return chr::duration_cast( + chr::duration>(value)); + } + + static constexpr auto sec(const double value) { + return chr::duration_cast( + chr::duration>(value)); + } + + kep3_DLL_PUBLIC friend kep_clock::duration operator-(const epoch &lhs, + const epoch &rhs); + + private: + + // Constructor for const time_point&) + explicit epoch(const kep_clock::time_point &time_point); + + // Constructor for const time_point&&) + explicit epoch(kep_clock::time_point &&time_point); + // Serialization code friend class boost::serialization::access; - template void serialize(Archive &ar, const unsigned int) { - ar &m_mjd2000; + template void serialize(Archive &ar, const uint) { + ar &boost::serialization::make_binary_object(&tp, sizeof(tp)); } // Serialization code (END) - /// the modified julian date 2000 stored in a double - double m_mjd2000; + // Time point relative to 1 Jan 2000 (MJD2000) + kep_clock::time_point tp; }; -kep3_DLL_PUBLIC epoch epoch_from_string(const std::string &date); -kep3_DLL_PUBLIC epoch epoch_from_iso_string(const std::string &date); - kep3_DLL_PUBLIC std::ostream &operator<<(std::ostream &s, const epoch &epoch_in); - } // end of namespace kep3 -template <> struct fmt::formatter : ostream_formatter {}; +template <> struct fmt::formatter : fmt::ostream_formatter {}; -#endif // kep3_EPOCH_H +#endif // kep3_EPOCH_HPP diff --git a/include/kep3/planets/keplerian.hpp b/include/kep3/planets/keplerian.hpp index 66cf120c..0025e4f3 100644 --- a/include/kep3/planets/keplerian.hpp +++ b/include/kep3/planets/keplerian.hpp @@ -25,7 +25,6 @@ namespace kep3::udpla { class kep3_DLL_PUBLIC keplerian { kep3::epoch m_ref_epoch; - std::array, 2> m_pos_vel_0; std::string m_name; double m_mu_central_body; double m_mu_self; @@ -33,11 +32,11 @@ class kep3_DLL_PUBLIC keplerian { double m_safe_radius; double m_period; bool m_ellipse; + std::array, 2> m_pos_vel_0; friend class boost::serialization::access; template void serialize(Archive &ar, unsigned) { ar &m_ref_epoch; - ar &m_pos_vel_0; ar &m_name; ar &m_mu_central_body; ar &m_mu_self; @@ -45,15 +44,19 @@ class kep3_DLL_PUBLIC keplerian { ar &m_safe_radius; ar &m_period; ar &m_ellipse; + ar &m_safe_radius; + ar &m_pos_vel_0; } public: // NOTE: added_param contains mu_self, radius and safe_radius explicit keplerian(const epoch &ref_epoch, const std::array &par, double mu_central_body = 1., std::string name = "Unknown", + std::array added_params = {-1., -1., -1.}, kep3::elements_type el_t = kep3::elements_type::KEP_F); // Constructor from pos_vel explicit keplerian(const epoch &ref_epoch = kep3::epoch(0), + const std::array, 2> &pos_vel = {{{1.0, 0.0, 0.0}, {0., 1.0, 0.0}}}, double mu_central_body = 1., std::string name = "Unknown", diff --git a/src/epoch.cpp b/src/epoch.cpp index 40ca411f..ab6c7592 100644 --- a/src/epoch.cpp +++ b/src/epoch.cpp @@ -7,250 +7,148 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. -#include +#include #include +#include +#include +#include #include #include +#include +#include +#include +#include +#include -#include -#include -// This is set to the precision of the boost date library (microseconds is -// default, nanoseconds can be set when compiling boosts. Note that the code has -// not been checked in that case) -#define BOOST_DATE_PRECISION 1e-6 +#include "kep3/core_astro/convert_julian_dates.hpp" +#include "kep3/epoch.hpp" -namespace kep3 { -using ::boost::gregorian::date; -using ::boost::gregorian::greg_day; -using ::boost::gregorian::greg_month; -using ::boost::gregorian::greg_year; -using ::boost::posix_time::ptime; -using ::boost::posix_time::time_duration; -using ::boost::posix_time::time_from_string; -using ::boost::posix_time::from_iso_string; +namespace kep3 +{ - -/// Constructor. /** - * Constructs an epoch from a non-gregorian date. - * \param[in] epoch_in A double indicating the non-gregorian date - * \param[in] epoch_type One of [epoch::MJD2000, epoch::MJD, epoch::JD] + * @brief Constructs a default epoch . */ -epoch::epoch(const double &epoch_in, julian_type epoch_type) - : m_mjd2000(epoch_in) { - switch (epoch_type) { - case MJD2000: - break; - case MJD: - m_mjd2000 = mjd2mjd2000(epoch_in); - break; - case JD: - m_mjd2000 = jd2mjd2000(epoch_in); - break; - } -} +epoch::epoch() + : tp{kep_clock::ref_epoch} {} -/// Constructor. /** - * Constructs an epoch from a gregorian date. The time of the day is assumed to - * be midnight. \param[in] year The gregorian year \param[in] month The month of - * the year \param[in] day The day of the month - */ -epoch::epoch(const greg_day &day, const greg_month &month, - const greg_year &year) - : m_mjd2000(0) { - // member m_mjd2000 will be here assigned. - set_posix_time(ptime(date(year, month, day))); -} - -/// Constructor. -/** - * Constructs an epoch from a boost ptime object (posix time) - * \param[in] posix_time The posix_time - */ -epoch::epoch(const ptime &posix_time) { - time_duration dt = posix_time - ptime(date(2000, 1, 1)); - bool flag = false; - if (dt.is_negative()) { - flag = true; - dt = dt.invert_sign(); - } - double fr_secs = - static_cast(dt.fractional_seconds()) * BOOST_DATE_PRECISION; - m_mjd2000 = static_cast(dt.hours()) / 24.0 + - static_cast(dt.minutes()) / 1440.0 + - (static_cast(dt.seconds()) + fr_secs) / 86400.0; - if (flag) { - m_mjd2000 = -m_mjd2000; - } -} - -/// jd getter. -/** - * Returns the julian date - * - * @return double containing the julian date + * @brief Constructs an epoch from a non-Gregorian date. * + * @param[in] epoch_in A double indicating the number of days + since day 0 in the specified calendar. + * @param[in] epoch_type epoch::julian_type */ -double epoch::jd() const { return mjd20002jd(m_mjd2000); } +epoch::epoch(const double epoch_in, const julian_type epoch_type) + : tp{make_tp(epoch_in, epoch_type)} {} -/// mjd getter. /** - * Returns the modified julian date - * - * @return double containing the modified julian date + * @brief Constructs an epoch from offsets relative to 0 MJD2000. * + * @param[in] y The number of years. + * @param[in] d The number of days. + * @param[in] h The number of hours. + * @param[in] min The number of minutes. + * @param[in] s The number of seconds. + * @param[in] ms The number of milliseconds. + * @param[in] us The number of microseconds. */ -double epoch::mjd() const { return mjd20002mjd(m_mjd2000); } +epoch::epoch(const int y, const uint mon, const uint d, const int h, + const int min, const int s, const int ms, const int us) + : tp{make_tp(y, mon, d, h, min, s, ms, us)} {} -/// mjd2000 getter. /** - * Gets the modified julian date 2000 - * @return const reference to mjd2000 - */ -double epoch::mjd2000() const { return m_mjd2000; } - -/// Extracts the posix time -/** - * Returns the posix_time representation of the epoch. The method evaluates - * from the mjd2000 the number of days, months, seconds and - * micro/nano seconds passed since the 1st of January 2000 and uses this - * information to build the posix time - * - * @return ptime containing the posix time + * @brief Constructs an epoch from a const reference to a time point. * + * @param[in] time_point Self-explanatory. */ -ptime epoch::get_posix_time() const { - long hrs = 0., min = 0., sec = 0., fsec = 0.; - bool flag = false; - double copy = m_mjd2000; - if (copy < 0) { - copy = -copy; - flag = true; - } - hrs = static_cast(copy * 24); - min = static_cast((copy * 24 - static_cast(hrs)) * 60); - sec = static_cast((((copy * 24 - static_cast(hrs)) * 60) - - static_cast(min)) * - 60); - double dblfsec = ((((copy * 24 - static_cast(hrs)) * 60) - - static_cast(min)) * - 60) - - static_cast(sec); - std::ostringstream fsecstr; - fsecstr << std::setiosflags(std::ios::fixed) - << std::setprecision( - -static_cast(std::log10(BOOST_DATE_PRECISION))) - << dblfsec; - fsec = boost::lexical_cast(fsecstr.str().substr( - 2, static_cast(-std::log10(BOOST_DATE_PRECISION) + 1))); - ptime retval; - if (flag) { - retval = ptime(date(2000, 1, 1), time_duration(-hrs, -min, -sec, -fsec)); - } else { - retval = ptime(date(2000, 1, 1), time_duration(hrs, min, sec, fsec)); - } - return retval; -} +epoch::epoch(const kep_clock::time_point &time_point) : tp{time_point} {} -/// Sets the epoch from a posix time /** - * Sets the epoch to a particular posix_time. - * - * \param[in] posix_time containing the posix time + * @brief Constructs an epoch from an rvalue reference to a time point. * + * @param[in] time_point Self-explanatory. */ -void epoch::set_posix_time(const ptime &posix_time) { +epoch::epoch(kep_clock::time_point &&time_point) : tp{std::move(time_point)} {} - m_mjd2000 = epoch(posix_time).mjd2000(); -} +kep_clock::time_point epoch::make_tp(const int y, const uint mon, const uint d, const int h, + const int min, const int s, const int ms, const int us) -epoch &epoch::operator+=(double rhs) { - /* addition of rhs to *this takes place here */ - m_mjd2000 += rhs; - return *this; // return the result by reference +{ + return kep_clock::time_point{} + chr::sys_days(chr::year_month_day{chr::year(y) / chr::month(mon) / chr::day(d)} + chr::months{0}).time_since_epoch() + + chr::hours(h) + chr::minutes(min) + chr::seconds(s) + + chr::milliseconds(ms) + chr::microseconds(us); } -epoch &epoch::operator-=(double rhs) { - /* addition of rhs to *this takes place here */ - m_mjd2000 -= rhs; - return *this; // return the result by reference +kep_clock::time_point epoch::make_tp(const double epoch_in, + const julian_type epoch_type) { + switch (epoch_type) { + case julian_type::MJD2000: + return epoch::tp_from_days(epoch_in); + case julian_type::MJD: + return epoch::tp_from_days(epoch_in) - 4453401600s; + case julian_type::JD: + return epoch::tp_from_days(epoch_in) - 211813444800s; + default: + throw; + } } -/// Returns an epoch constructed from a delimited string containing a date /** - * Builds an epoch from a delimited string. Excess digits in fractional seconds - * will be dropped. Ex: "1:02:03.123456999" => "1:02:03.123456". This behavior - * depends on the precision defined in astro_constant.h used to compile - * - * Example: - * std::string ts("2002-01-20 23:59:54.003"); - * epoch e(epoch_from_string(ts)) + * @brief Creates time point from the number of days since 0 MJD2000. * + * @return A time point */ -epoch epoch_from_string(const std::string &date) { - return epoch(time_from_string(date)); +constexpr kep_clock::time_point epoch::tp_from_days(const double days) { + return kep_clock::ref_epoch + + chr::duration_cast( + chr::duration>(days)); } -/// Returns an epoch constructed from a non delimited iso string containing a -/// date /** - * Builds an epoch from a non delimited iso string containing a date. - * - * Example: - * std::string ts("20020131T235959"); - * epoch e(epoch_from_iso_string(ts)) + * @brief Returns a time point formatted as a date/time string + * in the in the format 2000-12-31T12:34:56.123456. * + * @param tp The time point. + * @return A formatted date/time string. */ -epoch epoch_from_iso_string(const std::string &date) { - return epoch( - from_iso_string(date)); +std::string epoch::as_utc_string(const kep_clock::time_point& tp) { + std::stringstream iss; + const auto tse{tp.time_since_epoch()}; + const auto dp{std::chrono::duration_cast(tse)}; + const auto hms{std::chrono::duration_cast(tse - dp)}; + const auto us{std::chrono::duration_cast(tse - dp - hms)}; + iss << fmt::format("{:%F}", chr::sys_days(dp)) << "T" << fmt::format("{:%T}", hms) << "." << fmt::format("{:06}", us.count()); + return iss.str(); +} + +std::string epoch::as_utc_string() const { + return epoch::as_utc_string(tp); } -/// Overload the stream operator for kep_toolbox::epoch /** - * Streams out a date in the format 2000-Jan-01 00:12:30.123457 - * - * \param[in] s stream to which the epoch will be sent - * \param[in] now epoch to be sent to stream + * @brief Streams out an epoch as a UTC string. * - * \return reference to s + * @param[in] s Stream to which the epoch will be sent. + * @param[in] ep Epoch to be sent to the stream. * + * @return Reference to s. */ -std::ostream &operator<<(std::ostream &s, const epoch &epoch_in) { - s << epoch_in.get_posix_time(); +std::ostream &operator<<(std::ostream &s, const epoch &ep) { + s << ep.as_utc_string(); return s; } -bool operator>(const epoch &c1, const epoch &c2) { - return c1.m_mjd2000 > c2.m_mjd2000; -} -bool operator<(const epoch &c1, const epoch &c2) { - return c1.m_mjd2000 < c2.m_mjd2000; -} -bool operator>=(const epoch &c1, const epoch &c2) { - return c1.m_mjd2000 >= c2.m_mjd2000; -} -bool operator<=(const epoch &c1, const epoch &c2) { - return c1.m_mjd2000 <= c2.m_mjd2000; -} -bool operator==(const epoch &c1, const epoch &c2) { - return c1.m_mjd2000 == c2.m_mjd2000; -} -bool operator!=(const epoch &c1, const epoch &c2) { - return c1.m_mjd2000 != c2.m_mjd2000; -} -epoch operator+(epoch lhs, double rhs) { - lhs += rhs; // reuse compound assignment - return lhs; // return the result by value (uses move constructor) -} -epoch operator-(epoch lhs, double rhs) { - lhs -= rhs; // reuse compound assignment - return lhs; // return the result by value (uses move constructor) -} -double operator-(const epoch &lhs, const epoch &rhs) { - return lhs.mjd2000() - - rhs.mjd2000(); // return the result by value (uses move constructor) + +bool operator>(const epoch &c1, const epoch &c2) { return c1.tp > c2.tp; } +bool operator<(const epoch &c1, const epoch &c2) { return c1.tp < c2.tp; } +bool operator>=(const epoch &c1, const epoch &c2) { return c1.tp >= c2.tp; } +bool operator<=(const epoch &c1, const epoch &c2) { return c1.tp <= c2.tp; } +bool operator==(const epoch &c1, const epoch &c2) { return c1.tp == c2.tp; } +bool operator!=(const epoch &c1, const epoch &c2) { return c1.tp != c2.tp; } + +kep_clock::duration operator-(const epoch &lhs, const epoch &rhs) { + return lhs.tp - rhs.tp; } -} // namespace kep3 \ No newline at end of file +} // namespace kep3 diff --git a/src/planets/jpl_lp.cpp b/src/planets/jpl_lp.cpp index 44d40999..1ba83b31 100644 --- a/src/planets/jpl_lp.cpp +++ b/src/planets/jpl_lp.cpp @@ -179,7 +179,7 @@ double jpl_lp::get_radius() const { return m_radius; } double jpl_lp::get_safe_radius() const { return m_safe_radius; } std::string jpl_lp::get_extra_info() const { - kep3::epoch ep{0., kep3::epoch::MJD2000}; + kep3::epoch ep{0., kep3::epoch::julian_type::MJD2000}; auto par = elements(ep); auto pos_vel = eph(ep); diff --git a/src/planets/keplerian.cpp b/src/planets/keplerian.cpp index f926ae1b..e905923c 100644 --- a/src/planets/keplerian.cpp +++ b/src/planets/keplerian.cpp @@ -9,6 +9,8 @@ #include #include + +#include #include #include @@ -31,10 +33,10 @@ keplerian::keplerian(const epoch &ref_epoch, const std::array, 2> &pos_vel, double mu_central_body, std::string name, std::array added_params) - : m_ref_epoch(ref_epoch), m_pos_vel_0(pos_vel), m_name(std::move(name)), + : m_ref_epoch(ref_epoch), m_name(std::move(name)), m_mu_central_body(mu_central_body), m_mu_self(added_params[0]), m_radius(added_params[1]), m_safe_radius(added_params[2]), m_period(), - m_ellipse() { + m_ellipse(), m_pos_vel_0(pos_vel) { double R = std::sqrt(pos_vel[0][0] * pos_vel[0][0] + pos_vel[0][1] * pos_vel[0][1] + pos_vel[0][2] * pos_vel[0][2]); @@ -56,10 +58,10 @@ keplerian::keplerian(const epoch &ref_epoch, double mu_central_body, std::string name, std::array added_params, kep3::elements_type el_type) - : m_ref_epoch(ref_epoch), m_pos_vel_0(), m_name(std::move(name)), + : m_ref_epoch(ref_epoch), m_name(std::move(name)), m_mu_central_body(mu_central_body), m_mu_self(added_params[0]), m_radius(added_params[1]), m_safe_radius(added_params[2]), m_period(), - m_ellipse() { + m_ellipse(), m_pos_vel_0() { // orbital parameters a,e,i,W,w,f will be stored here std::array par(par_in); // we convert according to the chosen input @@ -103,7 +105,7 @@ keplerian::keplerian(const epoch &ref_epoch, std::array, 2> keplerian::eph(const kep3::epoch &ep) const { // 1 - We compute the dt - double dt = (ep.mjd2000() - m_ref_epoch.mjd2000()) * DAY2SEC; + auto dt = epoch::as_sec(ep - m_ref_epoch); // 2 - We propagate (make a copy as we do not want to change m_pos_vel_0) auto retval(m_pos_vel_0); kep3::propagate_lagrangian(retval, dt, m_mu_central_body); @@ -163,11 +165,11 @@ std::string keplerian::get_extra_info() const { retval += fmt::format("Mean anomly (deg.): {}\n", kep3::f2m(par[5], par[1]) * kep3::RAD2DEG); } - retval += fmt::format("Elements reference epoch (MJD2000): {}\n", - m_ref_epoch.mjd2000()) + - fmt::format("Elements reference epoch (date): {}\n", m_ref_epoch) + - fmt::format("r at ref. = {}\n", m_pos_vel_0[0]) + - fmt::format("v at ref. = {}\n", m_pos_vel_0[1]); + retval += + fmt::format("Elements reference epoch (MJD2000): {}\n", m_ref_epoch) + + fmt::format("Elements reference epoch (date): {}\n", m_ref_epoch) + + fmt::format("r at ref. = {}\n", m_pos_vel_0[0]) + + fmt::format("v at ref. = {}\n", m_pos_vel_0[1]); return retval; } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e077a613..9e53b038 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -20,6 +20,7 @@ function(ADD_kep3_TESTCASE arg1) "$<$:${kep3_CXX_FLAGS_RELEASE}>" "$<$:${kep3_CXX_FLAGS_RELEASE}>" ) + # Setup the C++ standard. target_compile_features(${arg1} PRIVATE cxx_std_20) set_property(TARGET ${arg1} PROPERTY CXX_EXTENSIONS NO) diff --git a/test/epoch_test.cpp b/test/epoch_test.cpp index 218ba098..3bac923f 100644 --- a/test/epoch_test.cpp +++ b/test/epoch_test.cpp @@ -7,7 +7,8 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. -#include +#include +#include #include #include #include @@ -17,65 +18,76 @@ #include "catch.hpp" using kep3::epoch; -using kep3::epoch_from_iso_string; -using kep3::epoch_from_string; -using boost::gregorian::date; -using boost::posix_time::ptime; +using kep3::kep_clock; +using namespace std::literals; +namespace chr = std::chrono; TEST_CASE("construct") { // test syntax + + // // > 2000 REQUIRE_NOTHROW(epoch()); REQUIRE_NOTHROW(epoch(123.456)); - REQUIRE_NOTHROW(epoch(123.456, epoch::MJD2000)); - REQUIRE_NOTHROW(epoch(123.456, epoch::JD)); - REQUIRE_NOTHROW(epoch(123.456, epoch::MJD)); - REQUIRE_NOTHROW(epoch(31, 12, 2034)); - // > 2000 - boost::posix_time::ptime posix_time_test(date(2034, 12, 31)); - REQUIRE_NOTHROW(epoch(posix_time_test)); - REQUIRE(epoch(posix_time_test).get_posix_time() == posix_time_test); - // < 2000 - boost::posix_time::ptime posix_time_test2(date(1980, 12, 31)); - REQUIRE_NOTHROW(epoch(posix_time_test2)); - REQUIRE(epoch(posix_time_test2).get_posix_time() == posix_time_test2); + REQUIRE_NOTHROW(epoch(123.456, epoch::julian_type::MJD2000)); + REQUIRE_NOTHROW(epoch(0.0, epoch::julian_type::JD)); + REQUIRE_NOTHROW(epoch(123.456, epoch::julian_type::JD)); + REQUIRE_NOTHROW(epoch(0.0, epoch::julian_type::MJD)); + REQUIRE_NOTHROW(epoch(123.456, epoch::julian_type::MJD)); + REQUIRE_NOTHROW(epoch(2000, 10, 17, 11, 36, 21, 121, 841)); + REQUIRE_NOTHROW(epoch(2064, 10, 17, 11, 36, 21, 121, 841).as_utc_string() == "2064-10-17T11:36:21.121841"); + + // // < 2000 + REQUIRE_NOTHROW(epoch(-123.456)); + REQUIRE_NOTHROW(epoch(-123.456, epoch::julian_type::MJD2000)); + REQUIRE_NOTHROW(epoch(-0.0, epoch::julian_type::JD)); + REQUIRE_NOTHROW(epoch(-123.456, epoch::julian_type::JD)); + REQUIRE_NOTHROW(epoch(-0.0, epoch::julian_type::MJD)); + REQUIRE_NOTHROW(epoch(-123.456, epoch::julian_type::MJD)); + REQUIRE_NOTHROW(epoch(1980, 10, 17, 11, 36, 21, 121, 841)); + REQUIRE_NOTHROW(epoch(1980, 10, 17, 11, 36, 21, 121, 841).as_utc_string() == "1980-10-17T11:36:21.121841"); + + // Epoch from lvalue and rvalue references + epoch ep{2000, 1, 1}; + REQUIRE(epoch(ep) == ep); + REQUIRE(epoch(epoch{2000,1,1}) == ep); // test conversions - REQUIRE(epoch(123.456).mjd2000() == epoch(123.456, epoch::MJD2000).mjd2000()); - REQUIRE(epoch(0.).mjd() == epoch(51544, epoch::MJD).mjd()); - REQUIRE(epoch(0.).jd() == epoch(2451544.5, epoch::JD).jd()); - REQUIRE(epoch(31, 12, 2034).jd() == epoch(posix_time_test).jd()); + REQUIRE(epoch(123.456).mjd2000() == + epoch(123.456, epoch::julian_type::MJD2000).mjd2000()); + REQUIRE(epoch(0.).mjd() == epoch(51544., epoch::julian_type::MJD).mjd()); + REQUIRE(epoch(0.).jd() == epoch(2451544.5, epoch::julian_type::JD).jd()); } TEST_CASE("epoch_operators") { - REQUIRE(epoch(0.) == epoch(0.)); - REQUIRE(epoch(0.) != epoch(1.)); + REQUIRE(epoch(2034, 10, 17) == epoch(2034, 10, 17)); + REQUIRE(epoch(2034, 10, 17) != epoch(2034, 11, 17)); + // Testing us precision + REQUIRE(epoch(2034, 10, 17) != epoch(2034, 10, 17, 0, 0, 0, 0, 1)); + // Check that ns precision is not supported + REQUIRE(epoch(2000, 10, 17) == epoch(2000, 10, 17, 0, 0, 0, 0, 0) + chr::nanoseconds(100)); + + // Conversion from double (defaults to days) REQUIRE(epoch(1.) > epoch(0.)); REQUIRE(epoch(1.) >= epoch(1.)); REQUIRE(epoch(1.) >= epoch(0.)); REQUIRE(epoch(0.) < epoch(1.)); REQUIRE(epoch(1.) <= epoch(1.)); epoch today(0.); - today += 100.; - REQUIRE(today == epoch(100.)); - today -= 100.; - REQUIRE(today == epoch(0.)); - auto yesterday = today - 1.; - REQUIRE(yesterday == epoch(-1.)); - today = yesterday + 1; - REQUIRE(today == epoch(0.)); - REQUIRE(today - yesterday == 1.); - REQUIRE_NOTHROW((std::cout << epoch())); -} + auto offset{chr::days(10963)}; + today += offset; + std::cout << "TODAY: " << today << "\n"; + REQUIRE(today == epoch(2030, 1, 6)); + today -= chr::duration_cast(offset); + REQUIRE(today == epoch()); + auto oneday{chr::days(1)}; + auto yesterday{today - chr::duration_cast(oneday)}; + auto yesterday1{today - oneday}; + REQUIRE(yesterday == yesterday1); -TEST_CASE("conversions_from_string") { - { - std::string ts("20020131T000000"); - epoch e(epoch_from_iso_string(ts)); - REQUIRE(e == epoch(boost::posix_time::ptime(date(2002, 01, 31)))); - } - { - std::string ts("2002-01-20 00:00:00.000"); - epoch e(epoch_from_string(ts)); - REQUIRE(e == epoch(boost::posix_time::ptime(date(2002, 01, 20)))); - } + REQUIRE(yesterday == epoch(1999, 12, 31)); + today = yesterday + chr::duration_cast(chr::days(1)); + REQUIRE(today == epoch()); + auto diff{today - yesterday}; + REQUIRE(diff == chr::duration_cast(chr::days(1))); + REQUIRE_NOTHROW((std::cout << epoch())); } diff --git a/test/planet_jpl_lp_test.cpp b/test/planet_jpl_lp_test.cpp index fc2b829d..52d71f47 100644 --- a/test/planet_jpl_lp_test.cpp +++ b/test/planet_jpl_lp_test.cpp @@ -24,7 +24,7 @@ using kep3::udpla::jpl_lp; TEST_CASE("constructor") { REQUIRE_NOTHROW(jpl_lp{}); - kep3::epoch ref_epoch{12.22, kep3::epoch::MJD2000}; + kep3::epoch ref_epoch{12.22, kep3::epoch::julian_type::MJD2000}; // From name REQUIRE_NOTHROW(jpl_lp{"Mars"}); REQUIRE_NOTHROW(jpl_lp{"mars"}); @@ -39,7 +39,7 @@ TEST_CASE("constructor") { TEST_CASE("eph") { // We use 2030-01-01 as a reference epoch for all these tests - kep3::epoch ref_epoch{2458849.5, kep3::epoch::JD}; + kep3::epoch ref_epoch{2458849.5, kep3::epoch::julian_type::JD}; { // This is Mercury w.r.t. the Sun queried from JPL Horizon at // 2020-01-01 @@ -136,12 +136,12 @@ TEST_CASE("eph") { REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 0.01); } jpl_lp udpla{"uranus"}; - REQUIRE_THROWS_AS(udpla.eph(kep3::epoch(5347534, kep3::epoch::MJD2000)), + REQUIRE_THROWS_AS(udpla.eph(kep3::epoch(5347534., kep3::epoch::julian_type::MJD2000)), std::domain_error); } TEST_CASE("elements") { - kep3::epoch ref_epoch{12.22, kep3::epoch::MJD2000}; + kep3::epoch ref_epoch{12.22, kep3::epoch::julian_type::MJD2000}; // We use Neptune jpl_lp udpla{"nePTUne"}; // casing is not important auto pos_vel = udpla.eph(ref_epoch); @@ -192,7 +192,7 @@ TEST_CASE("stream_operator") { TEST_CASE("serialization_test") { // Instantiate a generic jpl_lp - kep3::epoch ref_epoch{2423.4343, kep3::epoch::MJD2000}; + kep3::epoch ref_epoch{2423.4343, kep3::epoch::julian_type::MJD2000}; jpl_lp udpla{"neptune"}; // Store the string representation. diff --git a/test/planet_keplerian_test.cpp b/test/planet_keplerian_test.cpp index 9d94a045..a63837a6 100644 --- a/test/planet_keplerian_test.cpp +++ b/test/planet_keplerian_test.cpp @@ -18,13 +18,15 @@ #include #include "catch.hpp" +#include "kep3/epoch.hpp" #include "test_helpers.hpp" using kep3::udpla::keplerian; +using kep3::epoch; TEST_CASE("constructor") { REQUIRE_NOTHROW(keplerian{}); - kep3::epoch ref_epoch{12.22, kep3::epoch::MJD2000}; + kep3::epoch ref_epoch{12.22, kep3::epoch::julian_type::MJD2000}; // From posvel REQUIRE_NOTHROW( keplerian{ref_epoch, {{{0.3, 1., 0.2}, {0.0, 1.12, 0.}}}, 1., "unknown"}); @@ -124,7 +126,7 @@ TEST_CASE("constructor") { TEST_CASE("eph") { // We use 2000-01-01 as a reference epoch for all these tests - kep3::epoch ref_epoch{0., kep3::epoch::MJD2000}; + kep3::epoch ref_epoch{0., kep3::epoch::julian_type::MJD2000}; // This is a circular orbit at 1 AU. std::array, 2> pos_vel_0{ {{kep3::AU, 0., 0.}, {0., kep3::EARTH_VELOCITY, 0.}}}; @@ -132,13 +134,13 @@ TEST_CASE("eph") { keplerian udpla1{ref_epoch, pos_vel_0, kep3::MU_SUN}; double period_in_days = (2 * kep3::pi * kep3::AU / kep3::EARTH_VELOCITY) * kep3::SEC2DAY; - auto [r, v] = udpla1.eph(ref_epoch + period_in_days); + auto [r, v] = udpla1.eph(ref_epoch + epoch::days(period_in_days)); REQUIRE(kep3_tests::floating_point_error_vector(r, pos_vel_0[0]) < 1e-13); REQUIRE(kep3_tests::floating_point_error_vector(v, pos_vel_0[1]) < 1e-13); } TEST_CASE("elements") { - kep3::epoch ref_epoch{12.22, kep3::epoch::MJD2000}; + kep3::epoch ref_epoch{12.22, kep3::epoch::julian_type::MJD2000}; // Non singular elements std::array, 2> pos_vel{ {{1., 0.1, 0.1}, {0.1, 1., 0.1}}}; @@ -181,7 +183,7 @@ TEST_CASE("stream_operator") { TEST_CASE("serialization_test") { // Instantiate a generic udpla - kep3::epoch ref_epoch{2423.4343, kep3::epoch::MJD2000}; + kep3::epoch ref_epoch{2423.4343, kep3::epoch::julian_type::MJD2000}; keplerian udpla{ref_epoch, {{{0.33, 1.3, 0.12}, {0.01, 1.123, 0.2}}}, 1.12, @@ -196,6 +198,7 @@ TEST_CASE("serialization_test") { boost::archive::binary_oarchive oarchive(ss); oarchive << udpla; } + // Deserialize // Create a new udpla object keplerian udpla2{}; diff --git a/test/planet_test.cpp b/test/planet_test.cpp index 0a725b12..68035075 100644 --- a/test/planet_test.cpp +++ b/test/planet_test.cpp @@ -12,6 +12,10 @@ #include + +#include "catch.hpp" +#include + #include #include #include