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: 24827ed95742
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: eb45b55f5f02
Choose a head ref
  • 13 commits
  • 6 files changed
  • 2 contributors

Commits on Aug 23, 2020

  1. A first implementation.

    pleroy committed Aug 23, 2020
    Copy the full SHA
    20ba225 View commit details

Commits on Aug 24, 2020

  1. Copy the full SHA
    e138aa9 View commit details
  2. Copy the full SHA
    6441b3e View commit details
  3. Proper cast.

    pleroy committed Aug 24, 2020
    Copy the full SHA
    3f51e3e View commit details

Commits on Aug 25, 2020

  1. Copy the full SHA
    16d02ab View commit details
  2. A test that fails.

    pleroy committed Aug 25, 2020
    Copy the full SHA
    420feab View commit details
  3. A test that passes.

    pleroy committed Aug 25, 2020
    Copy the full SHA
    2cb8770 View commit details
  4. Merge.

    pleroy committed Aug 25, 2020
    Copy the full SHA
    541385c View commit details
  5. Cleanup.

    pleroy committed Aug 25, 2020
    Copy the full SHA
    9cd3647 View commit details
  6. A failed merge.

    pleroy committed Aug 25, 2020
    Copy the full SHA
    1633413 View commit details
  7. Cleanups.

    pleroy committed Aug 25, 2020
    Copy the full SHA
    8305052 View commit details
  8. Embarrassing.

    pleroy committed Aug 25, 2020
    Copy the full SHA
    228760a View commit details

Commits on Sep 5, 2020

  1. Merge pull request #2695 from pleroy/ToPWPS

    Extraction of a piecewise Poisson series from a continuous trajectory
    pleroy authored Sep 5, 2020
    Copy the full SHA
    eb45b55 View commit details
Showing with 156 additions and 12 deletions.
  1. +5 −3 base/not_null.hpp
  2. +2 −2 numerics/polynomial.hpp
  3. +4 −5 numerics/polynomial_body.hpp
  4. +18 −0 physics/continuous_trajectory.hpp
  5. +111 −2 physics/continuous_trajectory_body.hpp
  6. +16 −0 physics/continuous_trajectory_test.cpp
8 changes: 5 additions & 3 deletions base/not_null.hpp
Original file line number Diff line number Diff line change
@@ -318,9 +318,11 @@ template<typename Pointer>
std::ostream& operator<<(std::ostream& stream,
not_null<Pointer> const& pointer);

template<typename Result,
typename Pointer>
Result dynamic_cast_not_null(Pointer pointer);
template<typename Result, typename T>
not_null<Result> dynamic_cast_not_null(not_null<T*> const pointer);

template<typename Result, typename T>
not_null<Result> dynamic_cast_not_null(not_null<std::unique_ptr<T>>&& pointer);

} // namespace base
} // namespace principia
4 changes: 2 additions & 2 deletions numerics/polynomial.hpp
Original file line number Diff line number Diff line change
@@ -71,8 +71,8 @@ class Polynomial {
virtual Derivative<Value, Argument> EvaluateDerivative(
Argument const& argument) const = 0;

// Only useful for benchmarking or analyzing performance. Do not use in real
// code.
// Only useful for benchmarking, analyzing performance or for downcasting. Do
// not use in other circumstances.
virtual int degree() const = 0;

// Only useful for logging. Do not use in real code.
9 changes: 4 additions & 5 deletions numerics/polynomial_body.hpp
Original file line number Diff line number Diff line change
@@ -390,8 +390,8 @@ Derivative() const {
template<typename Value, typename Argument, int degree_,
template<typename, typename, int> class Evaluator>
template<typename, typename>
PolynomialInMonomialBasis<
Primitive<Value, Argument>, Argument, degree_ + 1, Evaluator>
PolynomialInMonomialBasis<Primitive<Value, Argument>, Argument,
degree_ + 1, Evaluator>
PolynomialInMonomialBasis<Value, Argument, degree_, Evaluator>::
Primitive() const {
return PolynomialInMonomialBasis<
@@ -543,9 +543,8 @@ Derivative() const {
template<typename Value, typename Argument, int degree_,
template<typename, typename, int> class Evaluator>
template<typename, typename>
PolynomialInMonomialBasis<
Primitive<Value, Argument>, Point<Argument>,
degree_ + 1, Evaluator>
PolynomialInMonomialBasis<Primitive<Value, Argument>, Point<Argument>,
degree_ + 1, Evaluator>
PolynomialInMonomialBasis<Value, Point<Argument>, degree_, Evaluator>::
Primitive() const {
return PolynomialInMonomialBasis<
18 changes: 18 additions & 0 deletions physics/continuous_trajectory.hpp
Original file line number Diff line number Diff line change
@@ -10,7 +10,9 @@
#include "base/not_null.hpp"
#include "base/status.hpp"
#include "geometry/named_quantities.hpp"
#include "numerics/poisson_series.hpp"
#include "numerics/polynomial.hpp"
#include "numerics/polynomial_evaluators.hpp"
#include "physics/checkpointer.hpp"
#include "physics/degrees_of_freedom.hpp"
#include "physics/trajectory.hpp"
@@ -29,6 +31,8 @@ using geometry::Position;
using geometry::Velocity;
using quantities::Length;
using quantities::Time;
using numerics::EstrinEvaluator;
using numerics::PiecewisePoissonSeries;
using numerics::Polynomial;

template<typename Frame>
@@ -88,6 +92,19 @@ class ContinuousTrajectory : public Trajectory<Frame> {

// End of the implementation of the interface.

// Returns the degree for a piecewise Poisson series covering the given time
// interval.
int PiecewisePoissonSeriesDegree(Instant const& t_min,
Instant const& t_max) const;

// Computes a piecewise Poisson series covering the given time interval. The
// degree must be at least the one returned by the preceding function.
template<int degree>
PiecewisePoissonSeries<Displacement<Frame>, degree, EstrinEvaluator>
ToPiecewisePoissonSeries(Instant const& t_min,
Instant const& t_max,
Instant const& origin) const;

void WriteToMessage(not_null<serialization::ContinuousTrajectory*> message)
const EXCLUDES(lock_);
template<typename F = Frame,
@@ -115,6 +132,7 @@ class ContinuousTrajectory : public Trajectory<Frame> {
// never need to extract their |t_min|. Logically, the |t_min| for a
// polynomial is the |t_max| of the previous one. The first polynomial has a
// |t_min| which is |*first_time_|.
// TODO(phl): These should be polynomials returning Position<Frame>.
struct InstantPolynomialPair {
InstantPolynomialPair(
Instant t_max,
113 changes: 111 additions & 2 deletions physics/continuous_trajectory_body.hpp
Original file line number Diff line number Diff line change
@@ -11,6 +11,7 @@
#include <vector>

#include "astronomy/epoch.hpp"
#include "geometry/interval.hpp"
#include "glog/stl_logging.h"
#include "numerics/newhall.hpp"
#include "numerics/polynomial_evaluators.hpp"
@@ -22,18 +23,21 @@ namespace principia {
namespace physics {
namespace internal_continuous_trajectory {

using base::dynamic_cast_not_null;
using base::Error;
using base::make_not_null_unique;
using geometry::Interval;
using numerics::EstrinEvaluator;
using numerics::PoissonSeries;
using numerics::ULPDistance;
using numerics::ЧебышёвSeries;
using quantities::DebugString;
using quantities::si::Metre;
using quantities::si::Second;
namespace si = quantities::si;

int const max_degree = 17;
int const min_degree = 3;
constexpr int max_degree = 17;
constexpr int min_degree = 3;
int const max_degree_age = 100;

// Only supports 8 divisions for now.
@@ -212,6 +216,111 @@ DegreesOfFreedom<Frame> ContinuousTrajectory<Frame>::EvaluateDegreesOfFreedom(
polynomial->EvaluateDerivative(time));
}

template<typename Frame>
int ContinuousTrajectory<Frame>::PiecewisePoissonSeriesDegree(
Instant const& t_min,
Instant const& t_max) const {
absl::ReaderMutexLock l(&lock_);
CHECK_LE(t_min_locked(), t_min);
CHECK_GE(t_max_locked(), t_max);
auto const it_min = FindPolynomialForInstant(t_min);
auto const it_max = FindPolynomialForInstant(t_max);
int degree = min_degree;
for (auto it = it_min;; ++it) {
degree = std::max(degree, it->polynomial->degree());
if (it == it_max) {
break;
}
}
return degree;
}


// Casts the polynomial to one of degree d1, and then increases the degree to
// d2.
#define PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, d1, d2) \
case d1: { \
if constexpr (d1 <= d2) { \
return PolynomialInMonomialBasis<Displacement<Frame>, Instant, \
d2, EstrinEvaluator>( \
*dynamic_cast_not_null< \
PolynomialInMonomialBasis<Displacement<Frame>, Instant, \
d1, EstrinEvaluator> const*>( \
polynomial)); \
} else { \
LOG(FATAL) << "Inconsistent degrees " << d1 << " and " << d2; \
} \
}

template<typename Frame>
template<int degree>
PiecewisePoissonSeries<Displacement<Frame>, degree, EstrinEvaluator>
ContinuousTrajectory<Frame>::ToPiecewisePoissonSeries(
Instant const& t_min,
Instant const& t_max,
Instant const& origin) const {
static_assert(degree >= min_degree && degree <= max_degree);
CHECK(!polynomials_.empty());
using PiecewisePoisson =
PiecewisePoissonSeries<Displacement<Frame>, degree, EstrinEvaluator>;
using Poisson = PoissonSeries<Displacement<Frame>, degree, EstrinEvaluator>;

auto cast_to_degree =
[](not_null<Polynomial<Displacement<Frame>, Instant> const*> const
polynomial)
-> PolynomialInMonomialBasis<Displacement<Frame>, Instant,
degree, EstrinEvaluator> {
switch (polynomial->degree()) {
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 3, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 4, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 5, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 6, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 7, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 8, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 9, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 10, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 11, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 12, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 13, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 14, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 15, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 16, degree);
PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS(polynomial, 17, degree);
default:
LOG(FATAL) << "Unexpected degree " << polynomial->degree();
}
};

std::unique_ptr<PiecewisePoisson> result;

absl::ReaderMutexLock l(&lock_);
auto const it_min = FindPolynomialForInstant(t_min);
auto const it_max = FindPolynomialForInstant(t_max);
Instant current_t_min = t_min;
for (auto it = it_min;; ++it) {
Instant const current_t_max = std::min(t_max, it->t_max);
Interval<Instant> interval;
interval.Include(current_t_min);
interval.Include(current_t_max);
auto const polynomial_cast_to_degree_at_origin =
cast_to_degree(it->polynomial.get()).AtOrigin(origin);
if (result == nullptr) {
result = std::make_unique<PiecewisePoisson>(
interval, Poisson(polynomial_cast_to_degree_at_origin, {{}}));
} else {
result->Append(interval,
Poisson(polynomial_cast_to_degree_at_origin, {{}}));
}
current_t_min = current_t_max;
if (it == it_max) {
break;
}
}
return *result;
}

#undef PRINCIPIA_CAST_TO_POLYNOMIAL_IN_MONOMIAL_BASIS

template<typename Frame>
void ContinuousTrajectory<Frame>::WriteToMessage(
not_null<serialization::ContinuousTrajectory*> const message) const {
16 changes: 16 additions & 0 deletions physics/continuous_trajectory_test.cpp
Original file line number Diff line number Diff line change
@@ -424,6 +424,8 @@ TEST_F(ContinuousTrajectoryTest, Polynomial) {
EXPECT_EQ(t0_ + (((number_of_steps - 1) / 8) * 8 + 1) * step,
trajectory->t_max());

// Check that the positions and velocities match the ones given by the
// functions above.
for (Instant time = trajectory->t_min();
time <= trajectory->t_max();
time += step / number_of_substeps) {
@@ -435,6 +437,20 @@ TEST_F(ContinuousTrajectoryTest, Polynomial) {
DegreesOfFreedom<World>(trajectory->EvaluatePosition(time),
trajectory->EvaluateVelocity(time)));
}

// Now check that it can be converted to a piecewise Poisson series.
Instant const t_min = trajectory->t_min() + 2 * step / number_of_substeps;
Instant const t_max = trajectory->t_max() - 3 * step / number_of_substeps;
EXPECT_EQ(3, trajectory->PiecewisePoissonSeriesDegree(t_min, t_max));
auto const piecewise_poisson_series =
trajectory->ToPiecewisePoissonSeries<3>(t_min, t_max, t0_);
EXPECT_EQ(t_min, piecewise_poisson_series.t_min());
EXPECT_EQ(t_max, piecewise_poisson_series.t_max());
for (Instant time = t_min; time <= t_max; time += step / number_of_substeps) {
EXPECT_THAT(
piecewise_poisson_series(time),
AlmostEquals(trajectory->EvaluatePosition(time) - World::origin, 0, 5));
}
}

// An approximation to the trajectory of Io.