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.

Sign up for GitHub

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 astronomy to DiscreteTraject0ry #3163

Merged
merged 3 commits into from
Oct 23, 2021
Merged

Convert astronomy to DiscreteTraject0ry #3163

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
527dd21
Astronomy compiles.
pleroy Oct 22, 2021
3e4bf92
Merge remote-tracking branch 'mockingbirdnest/Ondelette' into Astronomy
pleroy Oct 22, 2021
3f86696
Code cleanup.
pleroy Oct 22, 2021
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
37 changes: 17 additions & 20 deletions astronomy/geodesy_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
#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"
Expand All @@ -32,7 +33,7 @@ using integrators::methods::QuinlanTremaine1990Order12;
using physics::BodySurfaceDynamicFrame;
using physics::ContinuousTrajectory;
using physics::DegreesOfFreedom;
using physics::DiscreteTrajectory;
using physics::DiscreteTraject0ry;
using physics::Ephemeris;
using physics::KeplerianElements;
using physics::KeplerOrbit;
Expand Down Expand Up @@ -107,33 +108,29 @@ TEST_F(GeodesyTest, DISABLED_LAGEOS2) {
StandardProduct3::SatelliteIdentifier const lageos2_id{
StandardProduct3::SatelliteGroup::General, 52};

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

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_time, initial_dof_ilrsa] =
*initial_ilrsa.orbit(lageos2_id).front()->begin();

DegreesOfFreedom<ITRS> const initial_dof_ilrsb =
initial_ilrsb.orbit(lageos2_id).front()->front().degrees_of_freedom;
auto const& [_, initial_dof_ilrsb] =
*initial_ilrsb.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;
auto const& [final_time, expected_final_dof] =
*final_ilrsa.orbit(lageos2_id).front()->begin();

ephemeris_->Prolong(final_time);

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

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

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

DiscreteTrajectory<ICRS> trajectory;
DiscreteTraject0ry<ICRS> trajectory;
trajectory.Append(J2000, initial_state);
auto const instance = ephemeris_->NewInstance(
{&trajectory},
Expand All @@ -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.
DiscreteTrajectory<LunarSurface> surface_trajectory;
DiscreteTraject0ry<LunarSurface> surface_trajectory;
for (auto const& [time, degrees_of_freedom] : trajectory) {
surface_trajectory.Append(
time, lunar_frame_.ToThisFrameAtTime(time)(degrees_of_freedom));
Expand Down Expand Up @@ -413,18 +413,20 @@ TEST_P(LunarOrbitTest, NearCircularRepeatGroundTrackOrbit) {
mathematica::ExpressIn(Metre));
}

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

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

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

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

std::vector<double> descending_node_eccentricities;
Expand Down Expand Up @@ -530,7 +532,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];
Expand Down
2 changes: 0 additions & 2 deletions astronomy/mercury_perihelion_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@
#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"
Expand All @@ -33,7 +32,6 @@ 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;
Expand Down
9 changes: 5 additions & 4 deletions astronomy/orbit_analysis_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#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"
Expand All @@ -35,7 +36,7 @@ using numerics::EstrinEvaluator;
using numerics::PolynomialInMonomialBasis;
using physics::BodyCentredNonRotatingDynamicFrame;
using physics::BodySurfaceDynamicFrame;
using physics::DiscreteTrajectory;
using physics::DiscreteTraject0ry;
using physics::Ephemeris;
using physics::MasslessBody;
using physics::RotatingBody;
Expand Down Expand Up @@ -143,18 +144,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<DiscreteTrajectory<GCRS>>> EarthCentredTrajectory(
not_null<std::unique_ptr<DiscreteTraject0ry<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<DiscreteTrajectory<GCRS>>();
auto result = make_not_null_unique<DiscreteTraject0ry<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<DiscreteTrajectory<ITRS> const*>> const& orbit =
std::vector<not_null<DiscreteTraject0ry<ITRS> const*>> const& orbit =
sp3.orbit(sp3_orbit.satellite);
CHECK_EQ(orbit.size(), 1);
auto const& arc = *orbit.front();
Expand Down
6 changes: 3 additions & 3 deletions astronomy/orbit_ground_track.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#include "astronomy/orbit_recurrence.hpp"
#include "astronomy/orbital_elements.hpp"
#include "geometry/interval.hpp"
#include "physics/discrete_trajectory.hpp"
#include "physics/discrete_traject0ry.hpp"
#include "physics/rotating_body.hpp"

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

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

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

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

template<typename PrimaryCentred, typename Inertial>
absl::StatusOr<OrbitGroundTrack> OrbitGroundTrack::ForTrajectory(
DiscreteTrajectory<PrimaryCentred> const& trajectory,
DiscreteTraject0ry<PrimaryCentred> const& trajectory,
RotatingBody<Inertial> const& primary,
std::optional<MeanSun> const& mean_sun) {
DiscreteTrajectory<PrimaryCentred> ascending_nodes;
DiscreteTrajectory<PrimaryCentred> descending_nodes;
DiscreteTraject0ry<PrimaryCentred> ascending_nodes;
DiscreteTraject0ry<PrimaryCentred> descending_nodes;
OrbitGroundTrack ground_track;
RETURN_IF_ERROR(ComputeNodes(trajectory.begin(),
RETURN_IF_ERROR(ComputeNodes(trajectory,
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);
}
Expand All @@ -147,8 +148,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.front().time < ascending_nodes.front().time;
!ascending_nodes.empty() && !descending_nodes.empty() &&
descending_nodes.begin()->first < ascending_nodes.begin()->first;
return ground_track;
}

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

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

Expand Down Expand Up @@ -146,13 +146,13 @@ class OrbitalElements {

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

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

// |equinoctial_elements| must contain at least 2 elements.
static absl::StatusOr<Time> SiderealPeriod(
Expand Down
Loading