Skip to content

Commit

Permalink
Showing 14 changed files with 138 additions and 150 deletions.
37 changes: 20 additions & 17 deletions astronomy/geodesy_test.cpp
Original file line number Diff line number Diff line change
@@ -7,7 +7,6 @@
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "physics/body_surface_dynamic_frame.hpp"
#include "physics/discrete_traject0ry.hpp"
#include "physics/solar_system.hpp"
#include "quantities/si.hpp"
#include "testing_utilities/approximate_quantity.hpp"
@@ -33,7 +32,7 @@ using integrators::methods::QuinlanTremaine1990Order12;
using physics::BodySurfaceDynamicFrame;
using physics::ContinuousTrajectory;
using physics::DegreesOfFreedom;
using physics::DiscreteTraject0ry;
using physics::DiscreteTrajectory;
using physics::Ephemeris;
using physics::KeplerianElements;
using physics::KeplerOrbit;
@@ -108,29 +107,33 @@ TEST_F(GeodesyTest, DISABLED_LAGEOS2) {
StandardProduct3::SatelliteIdentifier const lageos2_id{
StandardProduct3::SatelliteGroup::General, 52};

CHECK_EQ(initial_ilrsa.orbit(lageos2_id).front()->begin()->first,
initial_ilrsb.orbit(lageos2_id).front()->begin()->first);
CHECK_EQ(initial_ilrsa.orbit(lageos2_id).front()->front().time,
initial_ilrsb.orbit(lageos2_id).front()->front().time);

auto const& [initial_time, initial_dof_ilrsa] =
*initial_ilrsa.orbit(lageos2_id).front()->begin();
Instant const initial_time =
initial_ilrsa.orbit(lageos2_id).front()->front().time;
DegreesOfFreedom<ITRS> const initial_dof_ilrsa =
initial_ilrsa.orbit(lageos2_id).front()->front().degrees_of_freedom;

auto const& [_, initial_dof_ilrsb] =
*initial_ilrsb.orbit(lageos2_id).front()->begin();
DegreesOfFreedom<ITRS> const initial_dof_ilrsb =
initial_ilrsb.orbit(lageos2_id).front()->front().degrees_of_freedom;

auto const& [final_time, expected_final_dof] =
*final_ilrsa.orbit(lageos2_id).front()->begin();
Instant const final_time =
final_ilrsa.orbit(lageos2_id).front()->front().time;
DegreesOfFreedom<ITRS> const expected_final_dof =
final_ilrsa.orbit(lageos2_id).front()->front().degrees_of_freedom;

ephemeris_->Prolong(final_time);

DiscreteTraject0ry<ICRS> primary_lageos2_trajectory;
DiscreteTrajectory<ICRS> primary_lageos2_trajectory;
primary_lageos2_trajectory.Append(
initial_time, itrs_.FromThisFrameAtTime(initial_time)(initial_dof_ilrsa));
DiscreteTraject0ry<ICRS> secondary_lageos2_trajectory;
DiscreteTrajectory<ICRS> secondary_lageos2_trajectory;
secondary_lageos2_trajectory.Append(
initial_time, itrs_.FromThisFrameAtTime(initial_time)(initial_dof_ilrsb));
auto flow_lageos2 =
[this, final_time](
DiscreteTraject0ry<ICRS>& lageos2_trajectory) -> absl::Status {
DiscreteTrajectory<ICRS>& lageos2_trajectory) -> absl::Status {
return ephemeris_->FlowWithAdaptiveStep(
&lageos2_trajectory,
Ephemeris<ICRS>::NoIntrinsicAcceleration,
@@ -152,13 +155,13 @@ TEST_F(GeodesyTest, DISABLED_LAGEOS2) {
return flow_lageos2(secondary_lageos2_trajectory);
});
bundle.Join();
EXPECT_THAT(primary_lageos2_trajectory.rbegin()->first, Eq(final_time));
EXPECT_THAT(secondary_lageos2_trajectory.rbegin()->first, Eq(final_time));
EXPECT_THAT(primary_lageos2_trajectory.back().time, Eq(final_time));
EXPECT_THAT(secondary_lageos2_trajectory.back().time, Eq(final_time));

auto const primary_actual_final_dof = itrs_.ToThisFrameAtTime(final_time)(
primary_lageos2_trajectory.rbegin()->second);
primary_lageos2_trajectory.back().degrees_of_freedom);
auto const secondary_actual_final_dof = itrs_.ToThisFrameAtTime(final_time)(
secondary_lageos2_trajectory.rbegin()->second);
secondary_lageos2_trajectory.back().degrees_of_freedom);

// Absolute error in position.
EXPECT_THAT(AbsoluteError(primary_actual_final_dof.position(),
26 changes: 12 additions & 14 deletions astronomy/lunar_orbit_test.cpp
Original file line number Diff line number Diff line change
@@ -18,7 +18,7 @@
#include "mathematica/mathematica.hpp"
#include "physics/apsides.hpp"
#include "physics/body_surface_dynamic_frame.hpp"
#include "physics/discrete_traject0ry.hpp"
#include "physics/discrete_trajectory.hpp"
#include "physics/kepler_orbit.hpp"
#include "physics/massless_body.hpp"
#include "physics/oblate_body.hpp"
@@ -57,7 +57,7 @@ using physics::BodySurfaceDynamicFrame;
using physics::ComputeApsides;
using physics::ComputeNodes;
using physics::DegreesOfFreedom;
using physics::DiscreteTraject0ry;
using physics::DiscreteTrajectory;
using physics::Ephemeris;
using physics::KeplerianElements;
using physics::KeplerOrbit;
@@ -363,7 +363,7 @@ TEST_P(LunarOrbitTest, NearCircularRepeatGroundTrackOrbit) {
IsNear(4.7e-13_⑴));
}

DiscreteTraject0ry<ICRS> trajectory;
DiscreteTrajectory<ICRS> trajectory;
trajectory.Append(J2000, initial_state);
auto const instance = ephemeris_->NewInstance(
{&trajectory},
@@ -377,7 +377,7 @@ TEST_P(LunarOrbitTest, NearCircularRepeatGroundTrackOrbit) {

// To find the nodes, we need to convert the trajectory to a reference frame
// whose xy plane is the Moon's equator.
DiscreteTraject0ry<LunarSurface> surface_trajectory;
DiscreteTrajectory<LunarSurface> surface_trajectory;
for (auto const& [time, degrees_of_freedom] : trajectory) {
surface_trajectory.Append(
time, lunar_frame_.ToThisFrameAtTime(time)(degrees_of_freedom));
@@ -413,20 +413,18 @@ TEST_P(LunarOrbitTest, NearCircularRepeatGroundTrackOrbit) {
mathematica::ExpressIn(Metre));
}

DiscreteTraject0ry<LunarSurface> ascending_nodes;
DiscreteTraject0ry<LunarSurface> descending_nodes;
ComputeNodes(surface_trajectory,
surface_trajectory.begin(),
DiscreteTrajectory<LunarSurface> ascending_nodes;
DiscreteTrajectory<LunarSurface> descending_nodes;
ComputeNodes(surface_trajectory.begin(),
surface_trajectory.end(),
/*north=*/Vector<double, LunarSurface>({0, 0, 1}),
/*max_points=*/std::numeric_limits<int>::max(),
ascending_nodes,
descending_nodes);

DiscreteTraject0ry<ICRS> apoapsides;
DiscreteTraject0ry<ICRS> periapsides;
DiscreteTrajectory<ICRS> apoapsides;
DiscreteTrajectory<ICRS> periapsides;
ComputeApsides(*ephemeris_->trajectory(moon_),
trajectory,
trajectory.begin(),
trajectory.end(),
/*max_points=*/std::numeric_limits<int>::max(),
@@ -435,12 +433,12 @@ TEST_P(LunarOrbitTest, NearCircularRepeatGroundTrackOrbit) {

struct Nodes {
std::string_view const name;
DiscreteTraject0ry<LunarSurface> const& trajectory;
DiscreteTrajectory<LunarSurface> const& trajectory;
};

struct Apsides {
std::string_view const name;
DiscreteTraject0ry<ICRS> const& trajectory;
DiscreteTrajectory<ICRS> const& trajectory;
};

std::vector<double> descending_node_eccentricities;
@@ -532,7 +530,7 @@ TEST_P(LunarOrbitTest, NearCircularRepeatGroundTrackOrbit) {
{
EccentricityVectorRange actual_period_ends;
for (int orbit = 0;
orbit < descending_nodes.size();
orbit < descending_nodes.Size();
orbit += orbits_per_period) {
auto& actual = actual_period_ends;
auto const e = descending_node_eccentricities[orbit];
2 changes: 2 additions & 0 deletions astronomy/mercury_perihelion_test.cpp
Original file line number Diff line number Diff line change
@@ -12,6 +12,7 @@
#include "integrators/symmetric_linear_multistep_integrator.hpp"
#include "mathematica/mathematica.hpp"
#include "physics/degrees_of_freedom.hpp"
#include "physics/discrete_trajectory.hpp"
#include "physics/ephemeris.hpp"
#include "physics/kepler_orbit.hpp"
#include "physics/massive_body.hpp"
@@ -32,6 +33,7 @@ using geometry::Position;
using integrators::SymmetricLinearMultistepIntegrator;
using integrators::methods::QuinlanTremaine1990Order12;
using physics::ContinuousTrajectory;
using physics::DiscreteTrajectory;
using physics::Ephemeris;
using physics::KeplerianElements;
using physics::KeplerOrbit;
9 changes: 4 additions & 5 deletions astronomy/orbit_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -14,7 +14,6 @@
#include "mathematica/mathematica.hpp"
#include "numerics/polynomial.hpp"
#include "physics/body_centred_non_rotating_dynamic_frame.hpp"
#include "physics/discrete_traject0ry.hpp"
#include "physics/ephemeris.hpp"
#include "physics/solar_system.hpp"
#include "testing_utilities/approximate_quantity.hpp"
@@ -36,7 +35,7 @@ using numerics::EstrinEvaluator;
using numerics::PolynomialInMonomialBasis;
using physics::BodyCentredNonRotatingDynamicFrame;
using physics::BodySurfaceDynamicFrame;
using physics::DiscreteTraject0ry;
using physics::DiscreteTrajectory;
using physics::Ephemeris;
using physics::MasslessBody;
using physics::RotatingBody;
@@ -144,18 +143,18 @@ class OrbitAnalysisTest : public ::testing::Test {

// Returns a GCRS trajectory obtained by stitching together the trajectories
// of |sp3_orbit.satellites| in |sp3_orbit.files|.
not_null<std::unique_ptr<DiscreteTraject0ry<GCRS>>> EarthCentredTrajectory(
not_null<std::unique_ptr<DiscreteTrajectory<GCRS>>> EarthCentredTrajectory(
SP3Orbit const& sp3_orbit) {
BodyCentredNonRotatingDynamicFrame<ICRS, GCRS> gcrs{ephemeris_.get(),
&earth_};
BodySurfaceDynamicFrame<ICRS, ITRS> itrs{ephemeris_.get(), &earth_};

auto result = make_not_null_unique<DiscreteTraject0ry<GCRS>>();
auto result = make_not_null_unique<DiscreteTrajectory<GCRS>>();
for (auto const& file : sp3_orbit.files.names) {
StandardProduct3 sp3(
SOLUTION_DIR / "astronomy" / "standard_product_3" / file,
sp3_orbit.files.dialect);
std::vector<not_null<DiscreteTraject0ry<ITRS> const*>> const& orbit =
std::vector<not_null<DiscreteTrajectory<ITRS> const*>> const& orbit =
sp3.orbit(sp3_orbit.satellite);
CHECK_EQ(orbit.size(), 1);
auto const& arc = *orbit.front();
6 changes: 3 additions & 3 deletions astronomy/orbit_ground_track.hpp
Original file line number Diff line number Diff line change
@@ -7,7 +7,7 @@
#include "astronomy/orbit_recurrence.hpp"
#include "astronomy/orbital_elements.hpp"
#include "geometry/interval.hpp"
#include "physics/discrete_traject0ry.hpp"
#include "physics/discrete_trajectory.hpp"
#include "physics/rotating_body.hpp"

namespace principia {
@@ -16,7 +16,7 @@ namespace internal_orbit_ground_track {

using geometry::Instant;
using geometry::Interval;
using physics::DiscreteTraject0ry;
using physics::DiscreteTrajectory;
using physics::RotatingBody;
using quantities::Angle;
using quantities::AngularFrequency;
@@ -67,7 +67,7 @@ class OrbitGroundTrack {
// sun-synchronicity is analysed.
template<typename PrimaryCentred, typename Inertial>
static absl::StatusOr<OrbitGroundTrack> ForTrajectory(
DiscreteTraject0ry<PrimaryCentred> const& trajectory,
DiscreteTrajectory<PrimaryCentred> const& trajectory,
RotatingBody<Inertial> const& primary,
std::optional<MeanSun> const& mean_sun);

35 changes: 17 additions & 18 deletions astronomy/orbit_ground_track_body.hpp
Original file line number Diff line number Diff line change
@@ -31,13 +31,14 @@ Angle CelestialLongitude(Position<PrimaryCentred> const& q) {
// The resulting angles are neither normalized nor unwound.
template<typename PrimaryCentred, typename Inertial>
std::vector<Angle> PlanetocentricLongitudes(
DiscreteTraject0ry<PrimaryCentred> const& nodes,
DiscreteTrajectory<PrimaryCentred> const& nodes,
RotatingBody<Inertial> const& primary) {
std::vector<Angle> longitudes;
longitudes.reserve(nodes.size());
for (auto const& [time, degrees_of_freedom] : nodes) {
longitudes.push_back(CelestialLongitude(degrees_of_freedom.position()) -
primary.AngleAt(time) - π / 2 * Radian);
longitudes.reserve(nodes.Size());
for (auto const& node : nodes) {
longitudes.push_back(
CelestialLongitude(node.degrees_of_freedom.position()) -
primary.AngleAt(node.time) - π / 2 * Radian);
}
return longitudes;
}
@@ -46,16 +47,15 @@ std::vector<Angle> PlanetocentricLongitudes(
template<typename Iterator>
Angle MeanSolarTime(Iterator const& it,
OrbitGroundTrack::MeanSun const& mean_sun) {
auto const& [time, degrees_of_freedom] = *it;
Time const t = time - mean_sun.epoch;
return π * Radian + CelestialLongitude(degrees_of_freedom.position()) -
Time const t = it->time - mean_sun.epoch;
return π * Radian + CelestialLongitude(it->degrees_of_freedom.position()) -
(mean_sun.mean_longitude_at_epoch +
(2 * π * Radian * t / mean_sun.year));
}

template<typename PrimaryCentred>
Interval<Angle> MeanSolarTimesOfNodes(
DiscreteTraject0ry<PrimaryCentred> const& nodes,
DiscreteTrajectory<PrimaryCentred> const& nodes,
OrbitGroundTrack::MeanSun const& mean_sun) {
Interval<Angle> mean_solar_times;
std::optional<Angle> mean_solar_time;
@@ -120,25 +120,24 @@ inline OrbitGroundTrack::EquatorCrossingLongitudes::EquatorCrossingLongitudes(

template<typename PrimaryCentred, typename Inertial>
absl::StatusOr<OrbitGroundTrack> OrbitGroundTrack::ForTrajectory(
DiscreteTraject0ry<PrimaryCentred> const& trajectory,
DiscreteTrajectory<PrimaryCentred> const& trajectory,
RotatingBody<Inertial> const& primary,
std::optional<MeanSun> const& mean_sun) {
DiscreteTraject0ry<PrimaryCentred> ascending_nodes;
DiscreteTraject0ry<PrimaryCentred> descending_nodes;
DiscreteTrajectory<PrimaryCentred> ascending_nodes;
DiscreteTrajectory<PrimaryCentred> descending_nodes;
OrbitGroundTrack ground_track;
RETURN_IF_ERROR(ComputeNodes(trajectory,
trajectory.begin(),
RETURN_IF_ERROR(ComputeNodes(trajectory.begin(),
trajectory.end(),
Vector<double, PrimaryCentred>({0, 0, 1}),
/*max_points=*/std::numeric_limits<int>::max(),
ascending_nodes,
descending_nodes));
if (mean_sun.has_value()) {
if (!ascending_nodes.empty()) {
if (!ascending_nodes.Empty()) {
ground_track.mean_solar_times_of_ascending_nodes_ =
MeanSolarTimesOfNodes(ascending_nodes, *mean_sun);
}
if (!descending_nodes.empty()) {
if (!descending_nodes.Empty()) {
ground_track.mean_solar_times_of_descending_nodes_ =
MeanSolarTimesOfNodes(descending_nodes, *mean_sun);
}
@@ -148,8 +147,8 @@ absl::StatusOr<OrbitGroundTrack> OrbitGroundTrack::ForTrajectory(
ground_track.longitudes_of_equator_crossings_of_descending_passes_ =
PlanetocentricLongitudes(descending_nodes, primary);
ground_track.first_descending_pass_before_first_ascending_pass_ =
!ascending_nodes.empty() && !descending_nodes.empty() &&
descending_nodes.begin()->first < ascending_nodes.begin()->first;
!ascending_nodes.Empty() && !descending_nodes.Empty() &&
descending_nodes.front().time < ascending_nodes.front().time;
return ground_track;
}

10 changes: 5 additions & 5 deletions astronomy/orbital_elements.hpp
Original file line number Diff line number Diff line change
@@ -6,7 +6,7 @@
#include "geometry/interval.hpp"
#include "geometry/named_quantities.hpp"
#include "physics/body.hpp"
#include "physics/discrete_traject0ry.hpp"
#include "physics/discrete_trajectory.hpp"
#include "physics/massive_body.hpp"
#include "quantities/named_quantities.hpp"
#include "quantities/quantities.hpp"
@@ -18,7 +18,7 @@ namespace internal_orbital_elements {
using geometry::Instant;
using geometry::Interval;
using physics::Body;
using physics::DiscreteTraject0ry;
using physics::DiscreteTrajectory;
using physics::MassiveBody;
using quantities::Angle;
using quantities::AngularFrequency;
@@ -37,7 +37,7 @@ class OrbitalElements {

template<typename PrimaryCentred>
static absl::StatusOr<OrbitalElements> ForTrajectory(
DiscreteTraject0ry<PrimaryCentred> const& trajectory,
DiscreteTrajectory<PrimaryCentred> const& trajectory,
MassiveBody const& primary,
Body const& secondary);

@@ -146,13 +146,13 @@ class OrbitalElements {

template<typename PrimaryCentred>
static std::vector<EquinoctialElements> OsculatingEquinoctialElements(
DiscreteTraject0ry<PrimaryCentred> const& trajectory,
DiscreteTrajectory<PrimaryCentred> const& trajectory,
MassiveBody const& primary,
Body const& secondary);

template<typename PrimaryCentred>
static std::vector<Length> RadialDistances(
DiscreteTraject0ry<PrimaryCentred> const& trajectory);
DiscreteTrajectory<PrimaryCentred> const& trajectory);

// |equinoctial_elements| must contain at least 2 elements.
static absl::StatusOr<Time> SiderealPeriod(
18 changes: 9 additions & 9 deletions astronomy/orbital_elements_body.hpp
Original file line number Diff line number Diff line change
@@ -32,13 +32,13 @@ using quantities::si::Radian;

template<typename PrimaryCentred>
absl::StatusOr<OrbitalElements> OrbitalElements::ForTrajectory(
DiscreteTraject0ry<PrimaryCentred> const& trajectory,
DiscreteTrajectory<PrimaryCentred> const& trajectory,
MassiveBody const& primary,
Body const& secondary) {
OrbitalElements orbital_elements;
if (trajectory.size() < 2) {
if (trajectory.Size() < 2) {
return absl::InvalidArgumentError(
"trajectory.size() is " + std::to_string(trajectory.size()));
"trajectory.Size() is " + std::to_string(trajectory.Size()));
}
orbital_elements.osculating_equinoctial_elements_ =
OsculatingEquinoctialElements(trajectory, primary, secondary);
@@ -65,8 +65,8 @@ absl::StatusOr<OrbitalElements> OrbitalElements::ForTrajectory(
"trajectory does not span one sidereal period: sidereal period is " +
DebugString(orbital_elements.sidereal_period_) +
", trajectory spans " +
DebugString(trajectory.rbegin()->first -
trajectory.begin()->first));
DebugString(trajectory.back().time -
trajectory.front().time));
}
auto mean_classical_elements =
ToClassicalElements(orbital_elements.mean_equinoctial_elements_);
@@ -138,13 +138,13 @@ inline Interval<Length> OrbitalElements::radial_distance_interval() const {
template<typename PrimaryCentred>
std::vector<OrbitalElements::EquinoctialElements>
OrbitalElements::OsculatingEquinoctialElements(
DiscreteTraject0ry<PrimaryCentred> const& trajectory,
DiscreteTrajectory<PrimaryCentred> const& trajectory,
MassiveBody const& primary,
Body const& secondary) {
DegreesOfFreedom<PrimaryCentred> const primary_dof{
PrimaryCentred::origin, PrimaryCentred::unmoving};
std::vector<EquinoctialElements> result;
result.reserve(trajectory.size());
result.reserve(trajectory.Size());
for (auto const& [time, degrees_of_freedom] : trajectory) {
auto const osculating_elements =
KeplerOrbit<PrimaryCentred>(primary,
@@ -175,9 +175,9 @@ OrbitalElements::OsculatingEquinoctialElements(

template<typename PrimaryCentred>
std::vector<Length> OrbitalElements::RadialDistances(
DiscreteTraject0ry<PrimaryCentred> const& trajectory) {
DiscreteTrajectory<PrimaryCentred> const& trajectory) {
std::vector<Length> radial_distances;
radial_distances.reserve(trajectory.size());
radial_distances.reserve(trajectory.Size());
DegreesOfFreedom<PrimaryCentred> const primary_dof{PrimaryCentred::origin,
PrimaryCentred::unmoving};
for (auto const& [time, degrees_of_freedom] : trajectory) {
9 changes: 4 additions & 5 deletions astronomy/orbital_elements_test.cpp
Original file line number Diff line number Diff line change
@@ -10,7 +10,6 @@
#include "gtest/gtest.h"
#include "mathematica/mathematica.hpp"
#include "physics/body_centred_non_rotating_dynamic_frame.hpp"
#include "physics/discrete_traject0ry.hpp"
#include "physics/ephemeris.hpp"
#include "physics/kepler_orbit.hpp"
#include "physics/solar_system.hpp"
@@ -35,7 +34,7 @@ using integrators::methods::DormandالمكاوىPrince1986RKN434FM;
using integrators::methods::QuinlanTremaine1990Order12;
using physics::BodyCentredNonRotatingDynamicFrame;
using physics::DegreesOfFreedom;
using physics::DiscreteTraject0ry;
using physics::DiscreteTrajectory;
using physics::Ephemeris;
using physics::KeplerianElements;
using physics::KeplerOrbit;
@@ -77,7 +76,7 @@ class OrbitalElementsTest : public ::testing::Test {

// Completes |initial_osculating_elements| and returns a GCRS trajectory
// obtained by flowing the corresponding initial conditions in |ephemeris|.
static not_null<std::unique_ptr<DiscreteTraject0ry<GCRS>>>
static not_null<std::unique_ptr<DiscreteTrajectory<GCRS>>>
EarthCentredTrajectory(
KeplerianElements<GCRS>& initial_osculating_elements,
Instant const& initial_time,
@@ -86,7 +85,7 @@ class OrbitalElementsTest : public ::testing::Test {
MassiveBody const& earth = FindEarthOrDie(ephemeris);
ephemeris.Prolong(final_time);
BodyCentredNonRotatingDynamicFrame<ICRS, GCRS> gcrs{&ephemeris, &earth};
DiscreteTraject0ry<ICRS> icrs_trajectory;
DiscreteTrajectory<ICRS> icrs_trajectory;
KeplerOrbit<GCRS> initial_osculating_orbit{earth,
MasslessBody{},
initial_osculating_elements,
@@ -110,7 +109,7 @@ class OrbitalElementsTest : public ::testing::Test {
/*speed_integration_tolerance=*/1 * Milli(Metre) / Second
},
/*max_ephemeris_steps=*/std::numeric_limits<std::int64_t>::max());
auto result = make_not_null_unique<DiscreteTraject0ry<GCRS>>();
auto result = make_not_null_unique<DiscreteTrajectory<GCRS>>();
for (auto const& [time, degrees_of_freedom] : icrs_trajectory) {
result->Append(time, gcrs.ToThisFrameAtTime(time)(degrees_of_freedom));
}
38 changes: 18 additions & 20 deletions astronomy/standard_product_3.cpp
Original file line number Diff line number Diff line change
@@ -32,25 +32,24 @@ using quantities::si::Second;
// n-point finite difference formulæ on the positions to produce a trajectory
// with consistent velocities.
template<int n>
not_null<std::unique_ptr<DiscreteTraject0ry<ITRS>>> ComputeVelocities(
DiscreteTraject0ry<ITRS> const& arc) {
auto result = make_not_null_unique<DiscreteTraject0ry<ITRS>>();
CHECK_GE(arc.size(), n);
not_null<std::unique_ptr<DiscreteTrajectory<ITRS>>> ComputeVelocities(
DiscreteTrajectory<ITRS> const& arc) {
auto result = make_not_null_unique<DiscreteTrajectory<ITRS>>();
CHECK_GE(arc.Size(), n);
std::array<Instant, n> times;
std::array<Position<ITRS>, n> positions;
auto it = arc.begin();
for (int k = 0; k < n; ++k, ++it) {
// TODO(egg): we should check when reading the file that the times are
// equally spaced at the interval declared in columns 25-38 of SP3
// line two.
auto const& [time, degrees_of_freedom] = *it;
times[k] = time;
positions[k] = degrees_of_freedom.position();
times[k] = it->time;
positions[k] = it->degrees_of_freedom.position();
}
// We use a central difference formula wherever possible, so we keep
// |offset| at (n - 1) / 2 except at the beginning and end of the arc.
int offset = 0;
for (int i = 0; i < arc.size(); ++i) {
for (int i = 0; i < arc.Size(); ++i) {
result->Append(
times[offset],
{positions[offset],
@@ -65,15 +64,14 @@ not_null<std::unique_ptr<DiscreteTraject0ry<ITRS>>> ComputeVelocities(
} else {
std::move(positions.begin() + 1, positions.end(), positions.begin());
std::move(times.begin() + 1, times.end(), times.begin());
auto const& [time, degrees_of_freedom] = *it;
times.back() = time;
positions.back() = degrees_of_freedom.position();
times.back() = it->time;
positions.back() = it->degrees_of_freedom.position();
++it;
}
}
// Note that having the right number of calls to |Append| does not guarantee
// this, as appending at an existing time merely emits a warning.
CHECK_EQ(result->size(), arc.size());
CHECK_EQ(result->Size(), arc.Size());
return result;
}

@@ -195,7 +193,7 @@ StandardProduct3::StandardProduct3(
std::forward_as_tuple());
CHECK(inserted) << "Duplicate satellite identifier " << id << ": "
<< full_location;
it->second.push_back(make_not_null_unique<DiscreteTraject0ry<ITRS>>());
it->second.push_back(make_not_null_unique<DiscreteTrajectory<ITRS>>());
satellites_.push_back(id);
} else {
CHECK_EQ(columns(c, c + 2), " 0") << full_location;
@@ -332,9 +330,9 @@ StandardProduct3::StandardProduct3(
// from the check.
CHECK_EQ(id, satellites_[i]) << location;

std::vector<not_null<std::unique_ptr<DiscreteTraject0ry<ITRS>>>>& orbit =
std::vector<not_null<std::unique_ptr<DiscreteTrajectory<ITRS>>>>& orbit =
it->second;
DiscreteTraject0ry<ITRS>& arc = *orbit.back();
DiscreteTrajectory<ITRS>& arc = *orbit.back();

Position<ITRS> const position =
Displacement<ITRS>({float_columns(5, 18) * Kilo(Metre),
@@ -377,8 +375,8 @@ StandardProduct3::StandardProduct3(

// Bad or absent positional and velocity values are to be set to 0.000000.
if (position == ITRS::origin || velocity == ITRS::unmoving) {
if (!arc.empty()) {
orbit.push_back(make_not_null_unique<DiscreteTraject0ry<ITRS>>());
if (!arc.Empty()) {
orbit.push_back(make_not_null_unique<DiscreteTrajectory<ITRS>>());
}
} else {
arc.Append(epoch, {position, velocity});
@@ -392,7 +390,7 @@ StandardProduct3::StandardProduct3(
for (auto& [id, orbit] : orbits_) {
// Do not leave a final empty trajectory if the orbit ends with missing
// data.
if (orbit.back()->empty()) {
if (orbit.back()->Empty()) {
orbit.pop_back();
}
}
@@ -405,7 +403,7 @@ StandardProduct3::StandardProduct3(
arc = ComputeVelocities<n>(*arc); \
break

switch (arc->size()) {
switch (arc->Size()) {
COMPUTE_VELOCITIES_CASE(1);
COMPUTE_VELOCITIES_CASE(2);
COMPUTE_VELOCITIES_CASE(3);
@@ -441,7 +439,7 @@ StandardProduct3::satellites() const {
return satellites_;
}

std::vector<not_null<DiscreteTraject0ry<ITRS> const*>> const&
std::vector<not_null<DiscreteTrajectory<ITRS> const*>> const&
StandardProduct3::orbit(SatelliteIdentifier const& id) const {
return FindOrDie(const_orbits_, id);
}
10 changes: 5 additions & 5 deletions astronomy/standard_product_3.hpp
Original file line number Diff line number Diff line change
@@ -9,7 +9,7 @@
#include "astronomy/frames.hpp"
#include "base/not_null.hpp"
#include "geometry/named_quantities.hpp"
#include "physics/discrete_traject0ry.hpp"
#include "physics/discrete_trajectory.hpp"

namespace principia {
namespace astronomy {
@@ -19,7 +19,7 @@ using base::not_null;
using geometry::Instant;
using geometry::Position;
using geometry::Velocity;
using physics::DiscreteTraject0ry;
using physics::DiscreteTrajectory;

// A representation of data in the extended standard product 3 orbit format.
// Specification:
@@ -111,7 +111,7 @@ class StandardProduct3 {

// Each orbit may consist of several arcs, separated by missing data.
// The arcs are non-overlapping, and are ordered chronologically.
std::vector<not_null<DiscreteTraject0ry<ITRS> const*>> const& orbit(
std::vector<not_null<DiscreteTrajectory<ITRS> const*>> const& orbit(
SatelliteIdentifier const& id) const;

Version version() const;
@@ -124,13 +124,13 @@ class StandardProduct3 {
private:
std::vector<SatelliteIdentifier> satellites_;
std::map<SatelliteIdentifier,
std::vector<not_null<std::unique_ptr<DiscreteTraject0ry<ITRS>>>>>
std::vector<not_null<std::unique_ptr<DiscreteTrajectory<ITRS>>>>>
orbits_;
// |orbits_| is the same as |const_orbits_|, but with non-owning pointers to
// constant trajectories; this allows us to return references to these vectors
// from |StandardProduct3::orbit|.
std::map<SatelliteIdentifier,
std::vector<not_null<DiscreteTraject0ry<ITRS> const*>>>
std::vector<not_null<DiscreteTrajectory<ITRS> const*>>>
const_orbits_;
Version version_;

66 changes: 29 additions & 37 deletions astronomy/standard_product_3_test.cpp
Original file line number Diff line number Diff line change
@@ -27,7 +27,7 @@ using integrators::methods::QuinlanTremaine1990Order12;
using physics::BodySurfaceDynamicFrame;
using physics::ContinuousTrajectory;
using physics::DegreesOfFreedom;
using physics::DiscreteTraject0ry;
using physics::DiscreteTrajectory;
using physics::Ephemeris;
using physics::RotatingBody;
using physics::SolarSystem;
@@ -362,48 +362,40 @@ TEST_P(StandardProduct3DynamicsTest, PerturbedKeplerian) {
EXPECT_THAT(sp3.version(), Eq(GetParam().version));
EXPECT_THAT(sp3.file_has_velocities(), Eq(GetParam().file_has_velocities));
for (auto const& satellite : sp3.satellites()) {
for (not_null<DiscreteTraject0ry<ITRS> const*> const arc :
for (not_null<DiscreteTrajectory<ITRS> const*> const arc :
sp3.orbit(satellite)) {
auto it = arc->begin();
for (int i = 0;; ++i) {
DiscreteTraject0ry<ICRS> integrated_arc;
{
auto const& [time, degrees_of_freedom] = *it;
ephemeris_->Prolong(time);
integrated_arc.Append(
time,
itrs_.FromThisFrameAtTime(time)(degrees_of_freedom));
}
DiscreteTrajectory<ICRS> integrated_arc;
ephemeris_->Prolong(it->time);
integrated_arc.Append(
it->time,
itrs_.FromThisFrameAtTime(it->time)(it->degrees_of_freedom));
if (++it == arc->end()) {
break;
}
{
auto const& [time, degrees_of_freedom] = *it;
ephemeris_->FlowWithAdaptiveStep(
&integrated_arc,
Ephemeris<ICRS>::NoIntrinsicAcceleration,
time,
Ephemeris<ICRS>::AdaptiveStepParameters(
EmbeddedExplicitRungeKuttaNyströmIntegrator<
DormandالمكاوىPrince1986RKN434FM,
Position<ICRS>>(),
std::numeric_limits<std::int64_t>::max(),
/*length_integration_tolerance=*/1 * Milli(Metre),
/*speed_integration_tolerance=*/1 * Milli(Metre) / Second),
/*max_ephemeris_steps=*/std::numeric_limits<std::int64_t>::max());
DegreesOfFreedom<ICRS> const actual =
integrated_arc.rbegin()->second;
DegreesOfFreedom<ICRS> const expected =
itrs_.FromThisFrameAtTime(time)(degrees_of_freedom);
EXPECT_THAT(AbsoluteError(expected.position(), actual.position()),
Lt(25 * Metre))
<< "orbit of satellite " << satellite << " flowing from point "
<< i;
EXPECT_THAT(AbsoluteError(expected.velocity(), actual.velocity()),
Lt(1 * Deci(Metre) / Second))
<< "orbit of satellite " << satellite << " flowing from point "
<< i;
}
ephemeris_->FlowWithAdaptiveStep(
&integrated_arc,
Ephemeris<ICRS>::NoIntrinsicAcceleration,
it->time,
Ephemeris<ICRS>::AdaptiveStepParameters(
EmbeddedExplicitRungeKuttaNyströmIntegrator<
DormandالمكاوىPrince1986RKN434FM,
Position<ICRS>>(),
std::numeric_limits<std::int64_t>::max(),
/*length_integration_tolerance=*/1 * Milli(Metre),
/*speed_integration_tolerance=*/1 * Milli(Metre) / Second),
/*max_ephemeris_steps=*/std::numeric_limits<std::int64_t>::max());
DegreesOfFreedom<ICRS> actual =
integrated_arc.back().degrees_of_freedom;
DegreesOfFreedom<ICRS> expected =
itrs_.FromThisFrameAtTime(it->time)(it->degrees_of_freedom);
EXPECT_THAT(AbsoluteError(expected.position(), actual.position()),
Lt(25 * Metre))
<< "orbit of satellite " << satellite << " flowing from point " << i;
EXPECT_THAT(AbsoluteError(expected.velocity(), actual.velocity()),
Lt(1 * Deci(Metre) / Second))
<< "orbit of satellite " << satellite << " flowing from point " << i;
}
}
}
16 changes: 7 additions & 9 deletions astronomy/лидов_古在_test.cpp
Original file line number Diff line number Diff line change
@@ -7,7 +7,6 @@
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "physics/body_centred_non_rotating_dynamic_frame.hpp"
#include "physics/discrete_traject0ry.hpp"
#include "physics/solar_system.hpp"
#include "mathematica/mathematica.hpp"
#include "testing_utilities/matchers.hpp"
@@ -26,7 +25,7 @@ using integrators::SymmetricLinearMultistepIntegrator;
using integrators::methods::Quinlan1999Order8A;
using integrators::methods::QuinlanTremaine1990Order12;
using physics::BodyCentredNonRotatingDynamicFrame;
using physics::DiscreteTraject0ry;
using physics::DiscreteTrajectory;
using physics::Ephemeris;
using physics::MassiveBody;
using physics::MasslessBody;
@@ -83,12 +82,11 @@ class Лидов古在Test : public ::testing::Test {

#if !_DEBUG
TEST_F(Лидов古在Test, MercuryOrbiter) {
DiscreteTraject0ry<ICRS> icrs_trajectory;
DiscreteTrajectory<ICRS> icrs_trajectory;
icrs_trajectory.Append(MercuryOrbiterInitialTime,
MercuryOrbiterInitialDegreesOfFreedom<ICRS>);
auto const ircs_segment = icrs_trajectory.segments().begin();
ircs_segment->SetDownsampling({.max_dense_intervals = 10'000,
.tolerance = 10 * Metre});
icrs_trajectory.SetDownsampling({.max_dense_intervals = 10'000,
.tolerance = 10 * Metre});
auto const instance = ephemeris_->NewInstance(
{&icrs_trajectory},
Ephemeris<ICRS>::NoIntrinsicAccelerations,
@@ -101,7 +99,7 @@ TEST_F(Лидов古在Test, MercuryOrbiter) {
LOG(INFO) << "Flowing to " << t;
auto const status = ephemeris_->FlowWithFixedStep(t, *instance);
if (!status.ok()) {
LOG(INFO) << status << " at " << icrs_trajectory.rbegin()->first;
LOG(INFO) << status << " at " << icrs_trajectory.back().time;
break;
}
}
@@ -110,13 +108,13 @@ TEST_F(Лидов古在Test, MercuryOrbiter) {
PRINCIPIA_UNICODE_PATH("лидов_古在.generated.wl"),
/*make_unique=*/false);

DiscreteTraject0ry<MercuryCentredInertial> mercury_centred_trajectory;
DiscreteTrajectory<MercuryCentredInertial> mercury_centred_trajectory;
for (auto const& [t, dof] : icrs_trajectory) {
mercury_centred_trajectory.Append(t,
mercury_frame_.ToThisFrameAtTime(t)(dof));
logger.Append(
"q",
mercury_centred_trajectory.rbegin()->second.position(),
mercury_centred_trajectory.back().degrees_of_freedom.position(),
mathematica::ExpressIn(Metre));
}

6 changes: 3 additions & 3 deletions astronomy/молния_orbit_test.cpp
Original file line number Diff line number Diff line change
@@ -12,7 +12,7 @@
#include "integrators/methods.hpp"
#include "integrators/symmetric_linear_multistep_integrator.hpp"
#include "mathematica/mathematica.hpp"
#include "physics/discrete_traject0ry.hpp"
#include "physics/discrete_trajectory.hpp"
#include "physics/kepler_orbit.hpp"
#include "physics/massless_body.hpp"
#include "physics/oblate_body.hpp"
@@ -37,7 +37,7 @@ using geometry::Position;
using integrators::SymmetricLinearMultistepIntegrator;
using integrators::methods::Quinlan1999Order8A;
using integrators::methods::QuinlanTremaine1990Order12;
using physics::DiscreteTraject0ry;
using physics::DiscreteTrajectory;
using physics::Ephemeris;
using physics::KeplerOrbit;
using physics::KeplerianElements;
@@ -118,7 +118,7 @@ TEST_F(МолнияOrbitTest, DISABLED_Satellite) {
*earth_body, satellite, initial_elements, J2000);
auto const satellite_state_vectors = initial_orbit.StateVectors(J2000);

DiscreteTraject0ry<ICRS> trajectory;
DiscreteTrajectory<ICRS> trajectory;
trajectory.Append(J2000, earth_degrees_of_freedom + satellite_state_vectors);
auto const instance = ephemeris_->NewInstance(
{&trajectory},

0 comments on commit d3bc6b8

Please sign in to comment.