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

Commits on Oct 23, 2021

  1. Copy the full SHA
    14978f4 View commit details
  2. Merge pull request #3167 from mockingbirdnest/revert-3162-Physics

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

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

template<typename T1, typename T2,
template<typename R,
typename = std::void_t<decltype(std::declval<R>().time)>,
typename = std::void_t<decltype(std::declval<R>().degrees_of_freedom)>,
typename OptionalExpressIn = std::nullopt_t>
std::string ToMathematica(std::pair<T1, T2> const& pair,
std::string ToMathematica(R ref,
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 T1, typename T2, typename OptionalExpressIn>
std::string ToMathematica(std::pair<T1, T2> const& pair,
template<typename R, typename, typename, typename OptionalExpressIn>
std::string ToMathematica(R const ref,
OptionalExpressIn express_in) {
return Apply("List",
std::vector<std::string>{
ToMathematica(pair.first, express_in),
ToMathematica(pair.second, express_in)});
ToMathematica(ref.time, express_in),
ToMathematica(ref.degrees_of_freedom, 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_traject0ry.hpp"
#include "physics/discrete_trajectory.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::DiscreteTraject0ry;
using physics::DiscreteTrajectory;
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));
}
{
DiscreteTraject0ry<F> trajectory;
DiscreteTrajectory<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::pair{trajectory.begin()->first,
trajectory.begin()->second}),
ToMathematica(*trajectory.begin()));
EXPECT_EQ(ToMathematica(std::tuple{trajectory.front().time,
trajectory.front().degrees_of_freedom}),
ToMathematica(trajectory.front()));
}
{
OrbitalElements::EquinoctialElements elements{
@@ -339,7 +339,7 @@ TEST_F(MathematicaTest, ExpressIn) {
ToMathematica(v, ExpressIn(Metre, Second)));
}
{
DiscreteTraject0ry<F> trajectory;
DiscreteTrajectory<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.begin(), ExpressIn(Metre, Second)));
ToMathematica(trajectory.front(), ExpressIn(Metre, Second)));
}
{
OrbitalElements::EquinoctialElements elements{
38 changes: 18 additions & 20 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_traject0ry.hpp"
#include "physics/discrete_trajectory.hpp"
#include "physics/trajectory.hpp"

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

// 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.
// 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.
template<typename Frame>
void ComputeApsides(Trajectory<Frame> const& reference,
DiscreteTraject0ry<Frame> const& trajectory,
typename DiscreteTraject0ry<Frame>::iterator begin,
typename DiscreteTraject0ry<Frame>::iterator end,
typename DiscreteTrajectory<Frame>::Iterator begin,
typename DiscreteTrajectory<Frame>::Iterator end,
int max_points,
DiscreteTraject0ry<Frame>& apoapsides,
DiscreteTraject0ry<Frame>& periapsides);
DiscreteTrajectory<Frame>& apoapsides,
DiscreteTrajectory<Frame>& periapsides);

// Computes the crossings of the section given by |begin| and |end| of
// |trajectory| with the xy plane. Appends the crossings that go towards the
// Computes the crossings of the discrete trajectory segment given by |begin|
// and |end| 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(DiscreteTraject0ry<Frame> const& trajectory,
typename DiscreteTraject0ry<Frame>::iterator begin,
typename DiscreteTraject0ry<Frame>::iterator end,
absl::Status ComputeNodes(typename DiscreteTrajectory<Frame>::Iterator begin,
typename DiscreteTrajectory<Frame>::Iterator end,
Vector<double, Frame> const& north,
int max_points,
DiscreteTraject0ry<Frame>& ascending,
DiscreteTraject0ry<Frame>& descending,
DiscreteTrajectory<Frame>& ascending,
DiscreteTrajectory<Frame>& descending,
Predicate predicate = Identically(true));

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

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

template<typename Frame>
void ComputeApsides(Trajectory<Frame> const& reference,
DiscreteTraject0ry<Frame> const& trajectory,
typename DiscreteTraject0ry<Frame>::iterator const begin,
typename DiscreteTraject0ry<Frame>::iterator const end,
typename DiscreteTrajectory<Frame>::Iterator const begin,
typename DiscreteTrajectory<Frame>::Iterator const end,
int const max_points,
DiscreteTraject0ry<Frame>& apoapsides,
DiscreteTraject0ry<Frame>& periapsides) {
DiscreteTrajectory<Frame>& apoapsides,
DiscreteTrajectory<Frame>& periapsides) {
std::optional<Instant> previous_time;
std::optional<DegreesOfFreedom<Frame>> previous_degrees_of_freedom;
std::optional<Square<Length>> previous_squared_distance;
@@ -113,13 +112,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 =
trajectory.EvaluateDegreesOfFreedom(apsis_time);
begin.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;
}
}
@@ -132,13 +131,12 @@ void ComputeApsides(Trajectory<Frame> const& reference,
}

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

DegreesOfFreedom<Frame> const node_degrees_of_freedom =
trajectory.EvaluateDegreesOfFreedom(node_time);
begin.trajectory()->EvaluateDegreesOfFreedom(node_time);
if (predicate(node_degrees_of_freedom)) {
if (Sign(InnerProduct(north, Vector<double, Frame>({0, 0, 1}))) ==
Sign(z_speed)) {
@@ -197,7 +195,7 @@ absl::Status ComputeNodes(DiscreteTraject0ry<Frame> const& trajectory,
} 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: 21 additions & 32 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, ComputeApsidesDiscreteTraject0ry) {
TEST_F(ApsidesTest, ComputeApsidesDiscreteTrajectory) {
Instant const t0;
GravitationalParameter const μ = SolarGravitationalParameter;
auto const b = new MassiveBody(μ);
@@ -89,7 +89,7 @@ TEST_F(ApsidesTest, ComputeApsidesDiscreteTraject0ry) {
Pow<3>(r_norm * Pow<2>(v_norm) - 2 * μ)));
Length const a = -r_norm * μ / (r_norm * Pow<2>(v_norm) - 2 * μ);

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

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

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

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

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

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

DiscreteTraject0ry<World> ascending_nodes;
DiscreteTraject0ry<World> descending_nodes;
ComputeNodes(trajectory,
trajectory.begin(),
DiscreteTrajectory<World> ascending_nodes;
DiscreteTrajectory<World> descending_nodes;
ComputeNodes(trajectory.begin(),
trajectory.end(),
north,
/*max_points=*/std::numeric_limits<int>::max(),
@@ -238,21 +236,20 @@ 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));

DiscreteTraject0ry<World> south_ascending_nodes;
DiscreteTraject0ry<World> south_descending_nodes;
DiscreteTrajectory<World> south_ascending_nodes;
DiscreteTrajectory<World> south_descending_nodes;
Vector<double, World> const mostly_south({1, 1, -1});
ComputeNodes(trajectory,
trajectory.begin(),
ComputeNodes(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(),
@@ -263,20 +260,12 @@ TEST_F(ApsidesTest, ComputeNodes) {
++ascending_it,
++south_descending_it,
++descending_it) {
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));
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));
}
}

3 changes: 1 addition & 2 deletions physics/body_centred_body_direction_dynamic_frame_test.cpp
Original file line number Diff line number Diff line change
@@ -12,7 +12,6 @@
#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"
@@ -453,7 +452,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.
DiscreteTraject0ry<ICRS> barycentre_trajectory;
DiscreteTrajectory<ICRS> barycentre_trajectory;
for (Time t; t <= period_; t += period_ / 16) {
auto const big_dof =
ephemeris_->trajectory(big_)->EvaluateDegreesOfFreedom(t0_ + t);
Loading