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: ae77ba9698ff
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 96a485ba3576
Choose a head ref
  • 4 commits
  • 11 files changed
  • 1 contributor

Commits on Oct 22, 2021

  1. Copy the full SHA
    c14d7ee View commit details
  2. DiscreteTraject0ry in physics.

    pleroy committed Oct 22, 2021
    Copy the full SHA
    2e895ef View commit details
  3. Copy the full SHA
    6d41161 View commit details
  4. Merge pull request #3162 from pleroy/Physics

    Convert mathematica and physics to DiscreteTraject0ry
    pleroy authored Oct 22, 2021
    Copy the full SHA
    96a485b View commit details
7 changes: 3 additions & 4 deletions mathematica/mathematica.hpp
Original file line number Diff line number Diff line change
@@ -6,6 +6,7 @@
#include <string>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>

#include "astronomy/orbital_elements.hpp"
@@ -221,11 +222,9 @@ template<typename Scalar, typename OptionalExpressIn = std::nullopt_t>
std::string ToMathematica(UnboundedVector<Scalar> const& vector,
OptionalExpressIn express_in = std::nullopt);

template<typename R,
typename = std::void_t<decltype(std::declval<R>().time)>,
typename = std::void_t<decltype(std::declval<R>().degrees_of_freedom)>,
template<typename T1, typename T2,
typename OptionalExpressIn = std::nullopt_t>
std::string ToMathematica(R ref,
std::string ToMathematica(std::pair<T1, T2> const& pair,
OptionalExpressIn express_in = std::nullopt);

template<typename V, typename A, int d,
8 changes: 4 additions & 4 deletions mathematica/mathematica_body.hpp
Original file line number Diff line number Diff line change
@@ -461,13 +461,13 @@ std::string ToMathematica(UnboundedVector<Scalar> const& vector,
return Apply("List", elements);
}

template<typename R, typename, typename, typename OptionalExpressIn>
std::string ToMathematica(R const ref,
template<typename T1, typename T2, typename OptionalExpressIn>
std::string ToMathematica(std::pair<T1, T2> const& pair,
OptionalExpressIn express_in) {
return Apply("List",
std::vector<std::string>{
ToMathematica(ref.time, express_in),
ToMathematica(ref.degrees_of_freedom, express_in)});
ToMathematica(pair.first, express_in),
ToMathematica(pair.second, express_in)});
}

template<typename V, typename A, int d,
16 changes: 8 additions & 8 deletions mathematica/mathematica_test.cpp
Original file line number Diff line number Diff line change
@@ -22,7 +22,7 @@
#include "numerics/polynomial_evaluators.hpp"
#include "numerics/unbounded_arrays.hpp"
#include "physics/degrees_of_freedom.hpp"
#include "physics/discrete_trajectory.hpp"
#include "physics/discrete_traject0ry.hpp"
#include "quantities/quantities.hpp"
#include "quantities/si.hpp"

@@ -52,7 +52,7 @@ using numerics::UnboundedLowerTriangularMatrix;
using numerics::UnboundedUpperTriangularMatrix;
using numerics::UnboundedVector;
using physics::DegreesOfFreedom;
using physics::DiscreteTrajectory;
using physics::DiscreteTraject0ry;
using quantities::Infinity;
using quantities::Length;
using quantities::Speed;
@@ -204,7 +204,7 @@ TEST_F(MathematicaTest, ToMathematica) {
EXPECT_EQ(ToMathematica(std::tuple{1.0, 2.0}), ToMathematica(v2));
}
{
DiscreteTrajectory<F> trajectory;
DiscreteTraject0ry<F> trajectory;
trajectory.Append(
Instant(),
DegreesOfFreedom<F>(
@@ -213,9 +213,9 @@ TEST_F(MathematicaTest, ToMathematica) {
Velocity<F>({-1.0 * Metre / Second,
-5.0 * Metre / Second,
8.0 * Metre / Second})));
EXPECT_EQ(ToMathematica(std::tuple{trajectory.front().time,
trajectory.front().degrees_of_freedom}),
ToMathematica(trajectory.front()));
EXPECT_EQ(ToMathematica(std::pair{trajectory.begin()->first,
trajectory.begin()->second}),
ToMathematica(*trajectory.begin()));
}
{
OrbitalElements::EquinoctialElements elements{
@@ -339,7 +339,7 @@ TEST_F(MathematicaTest, ExpressIn) {
ToMathematica(v, ExpressIn(Metre, Second)));
}
{
DiscreteTrajectory<F> trajectory;
DiscreteTraject0ry<F> trajectory;
trajectory.Append(
Instant(),
DegreesOfFreedom<F>(
@@ -352,7 +352,7 @@ TEST_F(MathematicaTest, ExpressIn) {
ToMathematica(std::tuple{0.0,
std::tuple{std::tuple{2.0, 3.0, -4.0},
std::tuple{-1.0, -5.0, 8.0}}}),
ToMathematica(trajectory.front(), ExpressIn(Metre, Second)));
ToMathematica(*trajectory.begin(), ExpressIn(Metre, Second)));
}
{
OrbitalElements::EquinoctialElements elements{
38 changes: 20 additions & 18 deletions physics/apsides.hpp
Original file line number Diff line number Diff line change
@@ -4,7 +4,7 @@

#include "absl/status/status.h"
#include "base/constant_function.hpp"
#include "physics/discrete_trajectory.hpp"
#include "physics/discrete_traject0ry.hpp"
#include "physics/trajectory.hpp"

namespace principia {
@@ -15,29 +15,31 @@ using base::ConstantFunction;
using base::Identically;
using geometry::Vector;

// Computes the apsides with respect to |reference| for the discrete trajectory
// segment given by |begin| and |end|. Appends to the given trajectories one
// point for each apsis.
// Computes the apsides with respect to |reference| for the section given by
// |begin| and |end| of |trajectory|. Appends to the given output trajectories
// one point for each apsis.
template<typename Frame>
void ComputeApsides(Trajectory<Frame> const& reference,
typename DiscreteTrajectory<Frame>::Iterator begin,
typename DiscreteTrajectory<Frame>::Iterator end,
DiscreteTraject0ry<Frame> const& trajectory,
typename DiscreteTraject0ry<Frame>::iterator begin,
typename DiscreteTraject0ry<Frame>::iterator end,
int max_points,
DiscreteTrajectory<Frame>& apoapsides,
DiscreteTrajectory<Frame>& periapsides);
DiscreteTraject0ry<Frame>& apoapsides,
DiscreteTraject0ry<Frame>& periapsides);

// Computes the crossings of the discrete trajectory segment given by |begin|
// and |end| with the xy plane. Appends the crossings that go towards the
// Computes the crossings of the section given by |begin| and |end| of
// |trajectory| with the xy plane. Appends the crossings that go towards the
// |north| side of the xy plane to |ascending|, and those that go away from the
// |north| side to |descending|.
// Nodes for which |predicate| returns false are excluded.
template<typename Frame, typename Predicate = ConstantFunction<bool>>
absl::Status ComputeNodes(typename DiscreteTrajectory<Frame>::Iterator begin,
typename DiscreteTrajectory<Frame>::Iterator end,
absl::Status ComputeNodes(DiscreteTraject0ry<Frame> const& trajectory,
typename DiscreteTraject0ry<Frame>::iterator begin,
typename DiscreteTraject0ry<Frame>::iterator end,
Vector<double, Frame> const& north,
int max_points,
DiscreteTrajectory<Frame>& ascending,
DiscreteTrajectory<Frame>& descending,
DiscreteTraject0ry<Frame>& ascending,
DiscreteTraject0ry<Frame>& descending,
Predicate predicate = Identically(true));

// TODO(egg): when we can usefully iterate over an arbitrary |Trajectory|, move
@@ -46,10 +48,10 @@ absl::Status ComputeNodes(typename DiscreteTrajectory<Frame>::Iterator begin,
template<typename Frame>
void ComputeApsides(Trajectory<Frame> const& trajectory1,
Trajectory<Frame> const& trajectory2,
DiscreteTrajectory<Frame>& apoapsides1,
DiscreteTrajectory<Frame>& periapsides1,
DiscreteTrajectory<Frame>& apoapsides2,
DiscreteTrajectory<Frame>& periapsides2);
DiscreteTraject0ry<Frame>& apoapsides1,
DiscreteTraject0ry<Frame>& periapsides1,
DiscreteTraject0ry<Frame>& apoapsides2,
DiscreteTraject0ry<Frame>& periapsides2);
#endif

} // namespace internal_apsides
26 changes: 14 additions & 12 deletions physics/apsides_body.hpp
Original file line number Diff line number Diff line change
@@ -29,11 +29,12 @@ using quantities::Variation;

template<typename Frame>
void ComputeApsides(Trajectory<Frame> const& reference,
typename DiscreteTrajectory<Frame>::Iterator const begin,
typename DiscreteTrajectory<Frame>::Iterator const end,
DiscreteTraject0ry<Frame> const& trajectory,
typename DiscreteTraject0ry<Frame>::iterator const begin,
typename DiscreteTraject0ry<Frame>::iterator const end,
int const max_points,
DiscreteTrajectory<Frame>& apoapsides,
DiscreteTrajectory<Frame>& periapsides) {
DiscreteTraject0ry<Frame>& apoapsides,
DiscreteTraject0ry<Frame>& periapsides) {
std::optional<Instant> previous_time;
std::optional<DegreesOfFreedom<Frame>> previous_degrees_of_freedom;
std::optional<Square<Length>> previous_squared_distance;
@@ -112,13 +113,13 @@ void ComputeApsides(Trajectory<Frame> const& reference,
// 3rd-degree polynomial would yield |squared_distance_approximation|, so
// we shouldn't be far from the truth.
DegreesOfFreedom<Frame> const apsis_degrees_of_freedom =
begin.trajectory()->EvaluateDegreesOfFreedom(apsis_time);
trajectory.EvaluateDegreesOfFreedom(apsis_time);
if (Sign(squared_distance_derivative).is_negative()) {
apoapsides.Append(apsis_time, apsis_degrees_of_freedom);
} else {
periapsides.Append(apsis_time, apsis_degrees_of_freedom);
}
if (apoapsides.Size() >= max_points && periapsides.Size() >= max_points) {
if (apoapsides.size() >= max_points && periapsides.size() >= max_points) {
break;
}
}
@@ -131,12 +132,13 @@ void ComputeApsides(Trajectory<Frame> const& reference,
}

template<typename Frame, typename Predicate>
absl::Status ComputeNodes(typename DiscreteTrajectory<Frame>::Iterator begin,
typename DiscreteTrajectory<Frame>::Iterator end,
absl::Status ComputeNodes(DiscreteTraject0ry<Frame> const& trajectory,
typename DiscreteTraject0ry<Frame>::iterator begin,
typename DiscreteTraject0ry<Frame>::iterator end,
Vector<double, Frame> const& north,
int const max_points,
DiscreteTrajectory<Frame>& ascending,
DiscreteTrajectory<Frame>& descending,
DiscreteTraject0ry<Frame>& ascending,
DiscreteTraject0ry<Frame>& descending,
Predicate predicate) {
static_assert(
std::is_convertible<decltype(predicate(
@@ -185,7 +187,7 @@ absl::Status ComputeNodes(typename DiscreteTrajectory<Frame>::Iterator begin,
}

DegreesOfFreedom<Frame> const node_degrees_of_freedom =
begin.trajectory()->EvaluateDegreesOfFreedom(node_time);
trajectory.EvaluateDegreesOfFreedom(node_time);
if (predicate(node_degrees_of_freedom)) {
if (Sign(InnerProduct(north, Vector<double, Frame>({0, 0, 1}))) ==
Sign(z_speed)) {
@@ -195,7 +197,7 @@ absl::Status ComputeNodes(typename DiscreteTrajectory<Frame>::Iterator begin,
} else {
descending.Append(node_time, node_degrees_of_freedom);
}
if (ascending.Size() >= max_points && descending.Size() >= max_points) {
if (ascending.size() >= max_points && descending.size() >= max_points) {
break;
}
}
53 changes: 32 additions & 21 deletions physics/apsides_test.cpp
Original file line number Diff line number Diff line change
@@ -56,7 +56,7 @@ class ApsidesTest : public ::testing::Test {

#if !defined(_DEBUG)

TEST_F(ApsidesTest, ComputeApsidesDiscreteTrajectory) {
TEST_F(ApsidesTest, ComputeApsidesDiscreteTraject0ry) {
Instant const t0;
GravitationalParameter const μ = SolarGravitationalParameter;
auto const b = new MassiveBody(μ);
@@ -89,7 +89,7 @@ TEST_F(ApsidesTest, ComputeApsidesDiscreteTrajectory) {
Pow<3>(r_norm * Pow<2>(v_norm) - 2 * μ)));
Length const a = -r_norm * μ / (r_norm * Pow<2>(v_norm) - 2 * μ);

DiscreteTrajectory<World> trajectory;
DiscreteTraject0ry<World> trajectory;
trajectory.Append(t0, DegreesOfFreedom<World>(World::origin + r, v));

ephemeris.FlowWithAdaptiveStep(
@@ -105,9 +105,10 @@ TEST_F(ApsidesTest, ComputeApsidesDiscreteTrajectory) {
1e-3 * Metre / Second),
Ephemeris<World>::unlimited_max_ephemeris_steps);

DiscreteTrajectory<World> apoapsides;
DiscreteTrajectory<World> periapsides;
DiscreteTraject0ry<World> apoapsides;
DiscreteTraject0ry<World> periapsides;
ComputeApsides(*ephemeris.trajectory(b),
trajectory,
trajectory.begin(),
trajectory.end(),
/*max_points=*/std::numeric_limits<int>::max(),
@@ -182,7 +183,7 @@ TEST_F(ApsidesTest, ComputeNodes) {
*ephemeris.bodies()[0], MasslessBody{}, elements, t0};
elements = orbit.elements_at_epoch();

DiscreteTrajectory<World> trajectory;
DiscreteTraject0ry<World> trajectory;
trajectory.Append(t0, initial_state[0] + orbit.StateVectors(t0));

ephemeris.FlowWithAdaptiveStep(
@@ -200,9 +201,10 @@ TEST_F(ApsidesTest, ComputeNodes) {

Vector<double, World> const north({0, 0, 1});

DiscreteTrajectory<World> ascending_nodes;
DiscreteTrajectory<World> descending_nodes;
ComputeNodes(trajectory.begin(),
DiscreteTraject0ry<World> ascending_nodes;
DiscreteTraject0ry<World> descending_nodes;
ComputeNodes(trajectory,
trajectory.begin(),
trajectory.end(),
north,
/*max_points=*/std::numeric_limits<int>::max(),
@@ -236,20 +238,21 @@ TEST_F(ApsidesTest, ComputeNodes) {
previous_time = time;
}

EXPECT_THAT(ascending_nodes.Size(), Eq(10));
EXPECT_THAT(descending_nodes.Size(), Eq(10));
EXPECT_THAT(ascending_nodes.size(), Eq(10));
EXPECT_THAT(descending_nodes.size(), Eq(10));

DiscreteTrajectory<World> south_ascending_nodes;
DiscreteTrajectory<World> south_descending_nodes;
DiscreteTraject0ry<World> south_ascending_nodes;
DiscreteTraject0ry<World> south_descending_nodes;
Vector<double, World> const mostly_south({1, 1, -1});
ComputeNodes(trajectory.begin(),
ComputeNodes(trajectory,
trajectory.begin(),
trajectory.end(),
mostly_south,
/*max_points=*/std::numeric_limits<int>::max(),
south_ascending_nodes,
south_descending_nodes);
EXPECT_THAT(south_ascending_nodes.Size(), Eq(10));
EXPECT_THAT(south_descending_nodes.Size(), Eq(10));
EXPECT_THAT(south_ascending_nodes.size(), Eq(10));
EXPECT_THAT(south_descending_nodes.size(), Eq(10));

for (auto south_ascending_it = south_ascending_nodes.begin(),
ascending_it = ascending_nodes.begin(),
@@ -260,12 +263,20 @@ TEST_F(ApsidesTest, ComputeNodes) {
++ascending_it,
++south_descending_it,
++descending_it) {
EXPECT_THAT(south_ascending_it->degrees_of_freedom,
Eq(descending_it->degrees_of_freedom));
EXPECT_THAT(south_ascending_it->time, Eq(descending_it->time));
EXPECT_THAT(south_descending_it->degrees_of_freedom,
Eq(ascending_it->degrees_of_freedom));
EXPECT_THAT(south_descending_it->time, Eq(ascending_it->time));
const auto& [ascending_time,
ascending_degrees_of_freedom] = *ascending_it;
const auto& [south_ascending_time,
south_ascending_degrees_of_freedom] = *south_ascending_it;
const auto& [descending_time,
descending_degrees_of_freedom] = *descending_it;
const auto& [south_descending_time,
south_descending_degrees_of_freedom] = *south_descending_it;
EXPECT_THAT(south_ascending_degrees_of_freedom,
Eq(descending_degrees_of_freedom));
EXPECT_THAT(south_ascending_time, Eq(descending_time));
EXPECT_THAT(south_descending_degrees_of_freedom,
Eq(ascending_degrees_of_freedom));
EXPECT_THAT(south_descending_time, Eq(ascending_time));
}
}

3 changes: 2 additions & 1 deletion physics/body_centred_body_direction_dynamic_frame_test.cpp
Original file line number Diff line number Diff line change
@@ -12,6 +12,7 @@
#include "gtest/gtest.h"
#include "integrators/methods.hpp"
#include "integrators/symplectic_runge_kutta_nyström_integrator.hpp"
#include "physics/discrete_traject0ry.hpp"
#include "physics/ephemeris.hpp"
#include "physics/mock_continuous_trajectory.hpp"
#include "physics/mock_ephemeris.hpp"
@@ -452,7 +453,7 @@ TEST_F(BodyCentredBodyDirectionDynamicFrameTest, ConstructFromOneBody) {
// A discrete trajectory that remains motionless at the barycentre. Since
// both bodies don't have the same mass, this means it has an intrinsic
// acceleration.
DiscreteTrajectory<ICRS> barycentre_trajectory;
DiscreteTraject0ry<ICRS> barycentre_trajectory;
for (Time t; t <= period_; t += period_ / 16) {
auto const big_dof =
ephemeris_->trajectory(big_)->EvaluateDegreesOfFreedom(t0_ + t);
Loading