Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: b2489b0f1997
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 684e0eda8860
Choose a head ref
  • 13 commits
  • 7 files changed
  • 1 contributor

Commits on Jun 19, 2021

  1. Copy the full SHA
    0e8cbd3 View commit details
  2. tabulation

    eggrobin committed Jun 19, 2021
    Copy the full SHA
    003b387 View commit details
  3. export

    eggrobin committed Jun 19, 2021
    Copy the full SHA
    6f38772 View commit details

Commits on Jun 22, 2021

  1. Copy the full SHA
    11bc7a4 View commit details

Commits on Jul 4, 2021

  1. merge

    eggrobin committed Jul 4, 2021
    Copy the full SHA
    cf1be14 View commit details
  2. Copy the full SHA
    bdb328d View commit details
  3. some review comments

    eggrobin committed Jul 4, 2021
    Copy the full SHA
    3236122 View commit details
  4. Copy the full SHA
    c72ccee View commit details
  5. FORCE_INLINE

    eggrobin committed Jul 4, 2021
    Copy the full SHA
    c4888cf View commit details
  6. using

    eggrobin committed Jul 4, 2021
    Copy the full SHA
    cb4bfd9 View commit details
  7. Copy the full SHA
    17e3df9 View commit details

Commits on Jul 5, 2021

  1. lint

    eggrobin committed Jul 5, 2021
    Copy the full SHA
    0df0769 View commit details
  2. Merge pull request #3034 from eggrobin/to_date

    Legible instants
    eggrobin authored Jul 5, 2021
    Copy the full SHA
    684e0ed View commit details
Showing with 126 additions and 18 deletions.
  1. +15 −0 astronomy/time_scales.hpp
  2. +45 −2 astronomy/time_scales_body.hpp
  3. +44 −0 astronomy/time_scales_test.cpp
  4. +4 −4 geometry/point.hpp
  5. +2 −2 geometry/point_body.hpp
  6. +8 −5 numerics/double_precision.hpp
  7. +8 −5 numerics/double_precision_body.hpp
15 changes: 15 additions & 0 deletions astronomy/time_scales.hpp
Original file line number Diff line number Diff line change
@@ -3,12 +3,15 @@

#include <string>

#include "astronomy/date_time.hpp"
#include "geometry/named_quantities.hpp"

namespace principia {
namespace astronomy {
namespace internal_time_scales {

using astronomy::date_time::Date;
using astronomy::date_time::DateTime;
using geometry::Instant;
using quantities::Angle;

@@ -72,6 +75,16 @@ constexpr Instant operator""_北斗(char const* str, std::size_t size);
Instant ParseGPSTime(std::string const& s);
Instant Parse北斗Time(std::string const& s);

// The following functions effectively round their argument toward InfinitePast
// to various granularities.

// Returns the TT day containing t.
constexpr Date TTDay(Instant const& t);
// Returns a DateTime representing the beginning of the TT second containing t.
constexpr DateTime TTSecond(Instant const& t);
// TODO(egg): TTMillisecond, but this is trickier because we use binary
// floating-point. UTC would be nice, too.

} // namespace internal_time_scales

using internal_time_scales::EarthRotationAngle;
@@ -81,6 +94,8 @@ using internal_time_scales::ParseTAI;
using internal_time_scales::ParseTT;
using internal_time_scales::ParseUT1;
using internal_time_scales::ParseUTC;
using internal_time_scales::TTDay;
using internal_time_scales::TTSecond;
using internal_time_scales::operator""_北斗;
using internal_time_scales::operator""_GPS;
using internal_time_scales::operator""_TAI;
47 changes: 45 additions & 2 deletions astronomy/time_scales_body.hpp
Original file line number Diff line number Diff line change
@@ -11,20 +11,23 @@
#include "astronomy/epoch.hpp"
#include "geometry/named_quantities.hpp"
#include "glog/logging.h"
#include "numerics/double_precision.hpp"
#include "quantities/si.hpp"

namespace principia {
namespace astronomy {
namespace internal_time_scales {

using astronomy::date_time::Date;
using astronomy::date_time::DateTime;
using astronomy::date_time::IsJulian;
using astronomy::date_time::JulianDate;
using astronomy::date_time::operator""_Date;
using astronomy::date_time::operator""_DateTime;
using astronomy::date_time::operator""_Julian;
using numerics::DoublePrecision;
using numerics::TwoDifference;
using quantities::si::Day;
using quantities::si::Hour;
using quantities::si::Minute;
using quantities::si::Radian;
using quantities::si::Second;

@@ -570,6 +573,46 @@ inline Instant Parse北斗Time(std::string const& s) {
return operator""_北斗(s.c_str(), s.size());
}

constexpr Date TTDay(Instant const& t) {
// We use a cast as a constexpr version of a floor; this is only correct with
// a positive JD, but we do not support dates before JD0.5 anyway.
CONSTEXPR_CHECK(t > "JD0.5"_TT);
std::int64_t jd_minus_half =
static_cast<std::int64_t>((t - "JD0.5"_TT) / Day);
// We want operations rounded toward negative infinity, but we also don’t
// want to fiddle with rounding modes. The product and sum here should be
// exact for all reasonable times.
if (jd_minus_half * Day + "JD0.5"_TT > t) {
--jd_minus_half;
}
return Date::JD(jd_minus_half + 0.5);
}

constexpr DateTime TTSecond(Instant const& t) {
auto const date = TTDay(t);
Instant const beginning_of_day = DateTimeAsTT(DateTime::BeginningOfDay(date));
// Close to J2000, Sterbenz’s lemma can fail to apply to this subtraction. We
// compute it exactly and round by hand (obtaining the preceding number by
// bit-casting and integer arithmetic).
DoublePrecision<quantities::Time> const time_of_day =
TwoDifference(t, beginning_of_day);
quantities::Time const time_of_day_rounded_down =
time_of_day.error >= 0 * Second
? time_of_day.value
: std::bit_cast<double>(
std::bit_cast<std::uint64_t>(time_of_day.value / Second) - 1) *
Second;
int const second_of_day = static_cast<int>(time_of_day_rounded_down / Second);
return DateTime(
date,
date_time::Time(
/*hour=*/second_of_day / (Hour / Second),
/*minute=*/second_of_day % static_cast<int>(Hour / Second) /
(Minute / Second),
/*second=*/second_of_day % static_cast<int>(Minute / Second),
/*millisecond=*/0));
}

} // namespace internal_time_scales
} // namespace astronomy
} // namespace principia
44 changes: 44 additions & 0 deletions astronomy/time_scales_test.cpp
Original file line number Diff line number Diff line change
@@ -29,6 +29,7 @@ using ::testing::AllOf;
using ::testing::Eq;
using ::testing::Gt;
using ::testing::Lt;
using ::testing::Ne;

constexpr Instant j2000_week = "1999-W52-6T12:00:00"_TT;

@@ -477,6 +478,49 @@ TEST_F(TimeScalesTest, GNSS) {
EXPECT_THAT("1999-08-22T00:00:13"_GPS, Eq("1999-08-22T00:00:00"_UTC));
}

TEST_F(TimeScalesTest, DateUnparsing) {
EXPECT_THAT(TTDay("1999-12-31T00:00:00"_TT), Eq("1999-12-31"_Date));
EXPECT_THAT(TTDay("2000-01-01T00:00:00"_TT), Eq("2000-01-01"_Date));
EXPECT_THAT(TTDay(J2000), Eq("2000-01-01"_Date));
EXPECT_THAT(TTDay("2000-01-01T23:59:59,999"_TT), Eq("2000-01-01"_Date));

// The floating-point division by Day turns the carriage into a pumpkin;
// adjustment is required.
constexpr double u = 0x1p-53;
constexpr Instant pumpkin = "2000-01-01T24:00:00"_TT;
static_assert(TTDay("2000-01-02T00:00:00"_TT) == "2000-01-02"_Date);
constexpr Instant carriage =
("2000-01-01T24:00:00"_TT - (pumpkin - J2000) * u);
static_assert(carriage != pumpkin);
static_assert(TTDay(carriage) != TTDay(pumpkin));
}

TEST_F(TimeScalesTest, DateTimeUnparsing) {
EXPECT_THAT(TTSecond(J2000), Eq("2000-01-01T12:00:00"_DateTime));
EXPECT_THAT(TTSecond("MJD0"_TT), Eq("1858-11-17T00:00:00"_DateTime));

{
// This is far enough from J2000 that Sterbenz’s lemma applies.
constexpr double u = 0x1p-53;
constexpr Instant pumpkin = "2000-01-02T24:00:00"_TT;
constexpr Instant carriage =
("2000-01-02T24:00:00"_TT - (pumpkin - J2000) * u);
static_assert(carriage != pumpkin);
static_assert(TTSecond(carriage) == "2000-01-02T23:59:59"_DateTime);
}

{
// Here again adjustment is needed to avoid the carriage turning into a
// pumpkin.
constexpr double u = 0x1p-53;
constexpr Instant pumpkin = "2000-01-01T24:00:00"_TT;
constexpr Instant carriage =
("2000-01-01T24:00:00"_TT - (pumpkin - J2000) * u);
static_assert(carriage != pumpkin);
static_assert(TTSecond(carriage) == "2000-01-01T23:59:59"_DateTime);
}
}

} // namespace internal_time_scales
} // namespace astronomy
} // namespace principia
8 changes: 4 additions & 4 deletions geometry/point.hpp
Original file line number Diff line number Diff line change
@@ -63,7 +63,8 @@ class Point final {
Vector coordinates_;

template<typename V>
friend Point<V> operator+(V const& translation, Point<V> const& point);
friend constexpr Point<V> operator+(V const& translation,
Point<V> const& point);
template<typename L, typename R>
friend Point<Product<L, R>> FusedMultiplyAdd(L const& a, R const& b,
Point<Product<L, R>> const& c);
@@ -91,10 +92,9 @@ class Point final {
friend class geometry::BarycentreCalculator;
};

// TODO(egg): constexpr these operators.
template<typename Vector>
Point<Vector> operator+(Vector const& translation,
Point<Vector> const& point);
constexpr Point<Vector> operator+(Vector const& translation,
Point<Vector> const& point);

template<typename L, typename R>
Point<Product<L, R>> FusedMultiplyAdd(L const& a,
4 changes: 2 additions & 2 deletions geometry/point_body.hpp
Original file line number Diff line number Diff line change
@@ -125,8 +125,8 @@ constexpr Point<Vector>::Point(Vector const& coordinates)
: coordinates_(coordinates) {}

template<typename Vector>
Point<Vector> operator+(Vector const& translation,
Point<Vector> const& point) {
constexpr Point<Vector> operator+(Vector const& translation,
Point<Vector> const& point) {
return point + translation;
}

13 changes: 8 additions & 5 deletions numerics/double_precision.hpp
Original file line number Diff line number Diff line change
@@ -61,16 +61,17 @@ DoublePrecision<Product<T, U>> TwoProduct(T const& a, U const& b);

// Same as |TwoProduct|, but never uses FMA.
template<typename T, typename U>
DoublePrecision<Product<T, U>> VeltkampDekkerProduct(T const& a, U const& b);
constexpr DoublePrecision<Product<T, U>> VeltkampDekkerProduct(T const& a,
U const& b);

// Computes the exact sum of a and b. The arguments must be such that
// |a| >= |b| or a == 0.
template<typename T, typename U>
DoublePrecision<Sum<T, U>> QuickTwoSum(T const& a, U const& b);
constexpr DoublePrecision<Sum<T, U>> QuickTwoSum(T const& a, U const& b);

// Computes the exact sum of a and b.
template<typename T, typename U>
DoublePrecision<Sum<T, U>> TwoSum(T const& a, U const& b);
constexpr DoublePrecision<Sum<T, U>> TwoSum(T const& a, U const& b);

// |TwoDifference| may have any of the following signatures:
// 1. Point × Point → Vector;
@@ -82,10 +83,12 @@ template<typename T,
typename U,
typename = Difference<T, Difference<T, U>>,
typename = std::enable_if_t<!std::is_same<U, Difference<U>>::value>>
DoublePrecision<Difference<T, U>> TwoDifference(T const& a, U const& b);
constexpr DoublePrecision<Difference<T, U>> TwoDifference(T const& a,
U const& b);

template<typename T, typename U, typename = Difference<Difference<T, U>, T>>
DoublePrecision<Difference<T, U>> TwoDifference(T const& a, U const& b);
constexpr DoublePrecision<Difference<T, U>> TwoDifference(T const& a,
U const& b);

DoublePrecision<Angle> Mod2π(DoublePrecision<Angle> const& θ);

13 changes: 8 additions & 5 deletions numerics/double_precision_body.hpp
Original file line number Diff line number Diff line change
@@ -211,7 +211,8 @@ DoublePrecision<Product<T, U>> Scale(T const & scale,
}

template<typename T, typename U>
DoublePrecision<Product<T, U>> VeltkampDekkerProduct(T const& a, U const& b) {
constexpr DoublePrecision<Product<T, U>> VeltkampDekkerProduct(T const& a,
U const& b) {
DoublePrecision<Product<T, U>> result;
auto const& x = a;
auto const& y = b;
@@ -251,7 +252,7 @@ DoublePrecision<Product<T, U>> TwoProduct(T const& a, U const& b) {
}

template<typename T, typename U>
FORCE_INLINE(inline)
FORCE_INLINE(constexpr)
DoublePrecision<Sum<T, U>> QuickTwoSum(T const& a, U const& b) {
#if _DEBUG
using quantities::DebugString;
@@ -269,7 +270,7 @@ DoublePrecision<Sum<T, U>> QuickTwoSum(T const& a, U const& b) {
}

template<typename T, typename U>
DoublePrecision<Sum<T, U>> TwoSum(T const& a, U const& b) {
constexpr DoublePrecision<Sum<T, U>> TwoSum(T const& a, U const& b) {
// [HLB07].
DoublePrecision<Sum<T, U>> result;
auto& s = result.value;
@@ -282,7 +283,8 @@ DoublePrecision<Sum<T, U>> TwoSum(T const& a, U const& b) {

// Point × Point → Vector.
template<typename T, typename U, typename, typename>
DoublePrecision<Difference<T, U>> TwoDifference(T const& a, U const& b) {
constexpr DoublePrecision<Difference<T, U>> TwoDifference(T const& a,
U const& b) {
static_assert(std::is_same<T, U>::value,
"Template metaprogramming went wrong");
using Point = T;
@@ -299,7 +301,8 @@ DoublePrecision<Difference<T, U>> TwoDifference(T const& a, U const& b) {

// Point × Vector → Point, or Vector × Vector → Vector.
template<typename T, typename U, typename>
DoublePrecision<Difference<T, U>> TwoDifference(T const& a, U const& b) {
constexpr DoublePrecision<Difference<T, U>> TwoDifference(T const& a,
U const& b) {
return TwoSum(a, -b);
}