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

Commits on Oct 23, 2021

  1. Copy the full SHA
    d3bc6b8 View commit details
  2. Merge pull request #3166 from mockingbirdnest/revert-3163-Astronomy

    Revert "Convert astronomy to DiscreteTraject0ry"
    eggrobin authored Oct 23, 2021
    Copy the full SHA
    601f720 View commit details
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(
Loading