Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Convert mathematica and physics to DiscreteTraject0ry #3162

Merged
merged 3 commits into from
Oct 22, 2021
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 3 additions & 4 deletions mathematica/mathematica.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#include <string>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>

#include "astronomy/orbital_elements.hpp"
Expand Down Expand Up @@ -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,
Expand Down
8 changes: 4 additions & 4 deletions mathematica/mathematica_body.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
16 changes: 8 additions & 8 deletions mathematica/mathematica_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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>(
Expand All @@ -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{
Expand Down Expand Up @@ -339,7 +339,7 @@ TEST_F(MathematicaTest, ExpressIn) {
ToMathematica(v, ExpressIn(Metre, Second)));
}
{
DiscreteTrajectory<F> trajectory;
DiscreteTraject0ry<F> trajectory;
trajectory.Append(
Instant(),
DegreesOfFreedom<F>(
Expand All @@ -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{
Expand Down
38 changes: 20 additions & 18 deletions physics/apsides.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand All @@ -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
Expand All @@ -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
Expand Down
26 changes: 14 additions & 12 deletions physics/apsides_body.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
}
}
Expand All @@ -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(
Expand Down Expand Up @@ -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)) {
Expand All @@ -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;
}
}
Expand Down
53 changes: 32 additions & 21 deletions physics/apsides_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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(μ);
Expand Down Expand Up @@ -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(
Expand All @@ -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(),
Expand Down Expand Up @@ -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(
Expand All @@ -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(),
Expand Down Expand Up @@ -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(),
Expand All @@ -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));
}
}

Expand Down
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
Expand Up @@ -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"
Expand Down Expand Up @@ -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);
Expand Down
Loading