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: 44d75edf04a4
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 72923f553e45
Choose a head ref
  • 14 commits
  • 10 files changed
  • 2 contributors

Commits on Apr 7, 2019

  1. Copy the full SHA
    f80f7e0 View commit details
  2. Copy the full SHA
    364d281 View commit details
  3. the test passes

    eggrobin committed Apr 7, 2019
    Copy the full SHA
    1e713df View commit details
  4. Copy the full SHA
    a7e38b6 View commit details
  5. add the new files

    eggrobin committed Apr 7, 2019
    Copy the full SHA
    69aa0f8 View commit details
  6. document

    eggrobin committed Apr 7, 2019
    Copy the full SHA
    582d7eb View commit details
  7. GRGS dialect

    eggrobin committed Apr 7, 2019
    Copy the full SHA
    36b196a View commit details
  8. has velocities

    eggrobin committed Apr 7, 2019
    Copy the full SHA
    103a2bf View commit details
  9. lint

    eggrobin committed Apr 7, 2019
    Copy the full SHA
    42f7f58 View commit details

Commits on Apr 8, 2019

  1. after phl's review

    eggrobin committed Apr 8, 2019
    Copy the full SHA
    6df4b69 View commit details
  2. const orbits

    eggrobin committed Apr 8, 2019
    Copy the full SHA
    1da06b5 View commit details
  3. lint

    eggrobin committed Apr 8, 2019
    Copy the full SHA
    13014a2 View commit details
  4. Copy the full SHA
    55e9de4 View commit details
  5. Merge pull request #2125 from eggrobin/grgs

    Computed velocities and GRGS dialect
    pleroy authored Apr 8, 2019
    Copy the full SHA
    72923f5 View commit details
3 changes: 3 additions & 0 deletions astronomy/astronomy.vcxproj
Original file line number Diff line number Diff line change
@@ -20,12 +20,15 @@
<None Include="standard_product_3\COD0MGXFIN_20181260000_01D_05M_ORB.SP3" />
<None Include="standard_product_3\COD0MGXFIN_20183640000_01D_05M_ORB.SP3" />
<None Include="standard_product_3\esa11802.eph" />
<None Include="standard_product_3\grgja203.b08243.e08247.D_S.sp3" />
<None Include="standard_product_3\grgtop03.b97344.e97348.D_S.sp3" />
<None Include="standard_product_3\ilrsa.orb.lageos2.160319.v35.sp3" />
<None Include="standard_product_3\ilrsa.orb.lageos2.180804.v70.sp3" />
<None Include="standard_product_3\ilrsb.orb.lageos2.160319.v35.sp3" />
<None Include="standard_product_3\mcc14000.sp3" />
<None Include="standard_product_3\nga20342.eph" />
<None Include="standard_product_3\README.md" />
<None Include="standard_product_3\ssaja102.b03007.e03017.DGS.sp3" />
<None Include="trappist_gravity_model.cfg" />
<None Include="trappist_gravity_model_slippist1.cfg" />
<None Include="trappist_initial_state_jd_2457000_000000000.cfg" />
9 changes: 9 additions & 0 deletions astronomy/astronomy.vcxproj.filters
Original file line number Diff line number Diff line change
@@ -92,6 +92,15 @@
<None Include="standard_product_3\README.md">
<Filter>Resource Files\Standard Product 3</Filter>
</None>
<None Include="standard_product_3\grgja203.b08243.e08247.D_S.sp3">
<Filter>Resource Files\Standard Product 3</Filter>
</None>
<None Include="standard_product_3\grgtop03.b97344.e97348.D_S.sp3">
<Filter>Resource Files\Standard Product 3</Filter>
</None>
<None Include="standard_product_3\ssaja102.b03007.e03017.DGS.sp3">
<Filter>Resource Files\Standard Product 3</Filter>
</None>
</ItemGroup>
<ItemGroup>
<Text Include="sol_initial_state_jd_2433282_500000000.proto.txt">
18 changes: 10 additions & 8 deletions astronomy/geodesy_test.cpp
Original file line number Diff line number Diff line change
@@ -73,7 +73,7 @@ class GeodesyTest : public ::testing::Test {
earth_trajectory_(*ephemeris_->trajectory(earth_)),
itrs_(ephemeris_.get(), earth_) {}

SolarSystem<ICRS> solar_system_2010_;
SolarSystem<ICRS> const solar_system_2010_;
not_null<std::unique_ptr<Ephemeris<ICRS>>> const ephemeris_;
not_null<OblateBody<ICRS> const*> const earth_;
ContinuousTrajectory<ICRS> const& earth_trajectory_;
@@ -106,19 +106,21 @@ TEST_F(GeodesyTest, LAGEOS2) {
StandardProduct3::SatelliteIdentifier const lageos2_id{
StandardProduct3::SatelliteGroup::General, 52};

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

Instant const initial_time = initial_ilrsa.orbit(lageos2_id).Begin().time();
Instant const initial_time =
initial_ilrsa.orbit(lageos2_id).front()->Begin().time();
DegreesOfFreedom<ITRS> const initial_dof_ilrsa =
initial_ilrsa.orbit(lageos2_id).Begin().degrees_of_freedom();
initial_ilrsa.orbit(lageos2_id).front()->Begin().degrees_of_freedom();

DegreesOfFreedom<ITRS> const initial_dof_ilrsb =
initial_ilrsb.orbit(lageos2_id).Begin().degrees_of_freedom();
initial_ilrsb.orbit(lageos2_id).front()->Begin().degrees_of_freedom();

Instant const final_time = final_ilrsa.orbit(lageos2_id).Begin().time();
Instant const final_time =
final_ilrsa.orbit(lageos2_id).front()->Begin().time();
DegreesOfFreedom<ITRS> const expected_final_dof =
final_ilrsa.orbit(lageos2_id).Begin().degrees_of_freedom();
final_ilrsa.orbit(lageos2_id).front()->Begin().degrees_of_freedom();

ephemeris_->Prolong(final_time);

156 changes: 142 additions & 14 deletions astronomy/standard_product_3.cpp
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
#include "astronomy/standard_product_3.hpp"

#include <algorithm>
#include <fstream>
#include <optional>
#include <string>
@@ -10,20 +11,70 @@
#include "astronomy/time_scales.hpp"
#include "base/map_util.hpp"
#include "glog/logging.h"
#include "numerics/finite_difference.hpp"

namespace principia {
namespace astronomy {
namespace internal_standard_product_3 {

using base::FindOrDie;
using base::make_not_null_unique;
using geometry::Displacement;
using numerics::FiniteDifference;
using quantities::NaN;
using quantities::Speed;
using quantities::si::Deci;
using quantities::si::Kilo;
using quantities::si::Metre;
using quantities::si::Second;

// Given a trajectory whose velocities are bad or absent (e.g., NaN), uses
// n-point finite difference formulæ on the positions to produce a trajectory
// with consistent velocities.
template<int 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.
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) {
result->Append(
times[offset],
{positions[offset],
FiniteDifference(
/*values=*/positions,
/*step=*/(times[n - 1] - times[0]) / (n - 1),
offset)});
// At every iteration, either |offset| advances, or the |positions|
// window shifts and |it| advances.
if (offset < (n - 1) / 2 || it == arc.End()) {
++offset;
} else {
std::move(positions.begin() + 1, positions.end(), positions.begin());
std::move(times.begin() + 1, times.end(), times.begin());
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());
return result;
}

StandardProduct3::StandardProduct3(
std::filesystem::path const& filename,
StandardProduct3::Dialect const dialect) {
@@ -136,11 +187,13 @@ StandardProduct3::StandardProduct3(
}
id.index = integer_columns(c + 1, c + 2);
CHECK_GT(id.index, 0) << full_location;
auto const [it, inserted] = orbits_.emplace(std::piecewise_construct, // NOLINT
std::forward_as_tuple(id),
std::forward_as_tuple());
auto const [it, inserted] = // NOLINT(whitespace/braces)
orbits_.emplace(std::piecewise_construct,
std::forward_as_tuple(id),
std::forward_as_tuple());
CHECK(inserted) << "Duplicate satellite identifier " << id << ": "
<< full_location;
it->second.push_back(make_not_null_unique<DiscreteTrajectory<ITRS>>());
satellites_.push_back(id);
} else {
CHECK_EQ(columns(c, c + 2), " 0") << full_location;
@@ -270,15 +323,18 @@ StandardProduct3::StandardProduct3(
// from the check.
CHECK_EQ(id, satellites_[i]) << location;

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

Position<ITRS> const position =
Displacement<ITRS>({float_columns(5, 18) * Kilo(Metre),
float_columns(19, 32) * Kilo(Metre),
float_columns(33, 46) * Kilo(Metre)}) +
ITRS::origin;
// TODO(egg): use a difference formula to compute the velocities if they
// are not provided.
// If the file does not provide velocities, fill the trajectory with NaN
// velocities; we then replace it with another trajectory whose velocities
// are computed using a finite difference formula.
Velocity<ITRS> velocity({NaN<Speed>(), NaN<Speed>(), NaN<Speed>()});

read_line();
@@ -295,10 +351,11 @@ StandardProduct3::StandardProduct3(
CHECK_EQ(SatelliteGroup{column(2)}, id.group) << location;
}
CHECK_EQ(integer_columns(3, 4), id.index) << location;
velocity =
Velocity<ITRS>({float_columns(5, 18) * (Deci(Metre) / Second),
float_columns(19, 32) * (Deci(Metre) / Second),
float_columns(33, 46) * (Deci(Metre) / Second)});
Speed const speed_unit =
dialect == Dialect::GRGS ? Metre / Second : Deci(Metre) / Second;
velocity = Velocity<ITRS>({float_columns(5, 18) * speed_unit,
float_columns(19, 32) * speed_unit,
float_columns(33, 46) * speed_unit});

read_line();
if (version_ >= Version::C && line.has_value() &&
@@ -309,30 +366,85 @@ StandardProduct3::StandardProduct3(
}
}

orbit.Append(epoch, {position, velocity});
// Bad or absent positional and velocity values are to be set to 0.000000.
if (position == ITRS::origin || velocity == Velocity<ITRS>()) {
if (!arc.Empty()) {
orbit.push_back(make_not_null_unique<DiscreteTrajectory<ITRS>>());
}
} else {
arc.Append(epoch, {position, velocity});
}
}
}
if (dialect != Dialect::ILRSA) {
CHECK_EQ(columns(1, 3), "EOF") << location;
read_line();
}
for (auto& [id, orbit] : orbits_) {
// Do not leave a final empty trajectory if the orbit ends with missing
// data.
if (orbit.back()->Empty()) {
orbit.pop_back();
}
}
CHECK(!line.has_value()) << location;
if (!has_velocities_) {
for (auto& [id, orbit] : orbits_) {
for (auto& arc : orbit) {
#define COMPUTE_VELOCITIES_CASE(n) \
case n: \
arc = ComputeVelocities<n>(*arc); \
break

switch (arc->Size()) {
COMPUTE_VELOCITIES_CASE(1);
COMPUTE_VELOCITIES_CASE(2);
COMPUTE_VELOCITIES_CASE(3);
COMPUTE_VELOCITIES_CASE(4);
COMPUTE_VELOCITIES_CASE(5);
COMPUTE_VELOCITIES_CASE(6);
COMPUTE_VELOCITIES_CASE(7);
COMPUTE_VELOCITIES_CASE(8);
default:
arc = ComputeVelocities<9>(*arc);
break;
}

#undef COMPUTE_VELOCITIES_CASE
}
}
}
for (auto& [id, orbit] : orbits_) {
auto const [it, inserted] = // NOLINT(whitespace/braces)
const_orbits_.emplace(std::piecewise_construct,
std::forward_as_tuple(id),
std::forward_as_tuple());
CHECK(inserted) << id;
auto& const_orbit = it->second;
for (auto& arc : orbit) {
const_orbit.push_back(arc.get());
}
}
}

std::vector<StandardProduct3::SatelliteIdentifier> const&
StandardProduct3::satellites() const {
return satellites_;
}

DiscreteTrajectory<ITRS> const& StandardProduct3::orbit(
SatelliteIdentifier const& id) const {
return FindOrDie(orbits_, id);
std::vector<not_null<DiscreteTrajectory<ITRS> const*>> const&
StandardProduct3::orbit(SatelliteIdentifier const& id) const {
return FindOrDie(const_orbits_, id);
}

StandardProduct3::Version StandardProduct3::version() const {
return version_;
}

bool StandardProduct3::file_has_velocities() const {
return has_velocities_;
}

bool operator==(StandardProduct3::SatelliteIdentifier const& left,
StandardProduct3::SatelliteIdentifier const& right) {
return left.group == right.group && left.index == right.index;
@@ -349,6 +461,22 @@ std::ostream& operator<<(std::ostream& out,
return out << std::string(1, static_cast<char>(version));
}

std::ostream& operator<<(std::ostream& out,
StandardProduct3::Dialect const& dialect) {
switch (dialect) {
case StandardProduct3::Dialect::Standard:
return out << "Standard SP3";
case StandardProduct3::Dialect::ILRSA:
return out << "ILRSA SP3 dialect";
case StandardProduct3::Dialect::ILRSB:
return out << "ILRSB SP3 dialect";
case StandardProduct3::Dialect::GRGS:
return out << "GRGS SP3 dialect";
default:
return out << "Unknown SP3 dialect (" << static_cast<int>(dialect) << ")";
}
}

std::ostream& operator<<(std::ostream& out,
StandardProduct3::SatelliteGroup const& group) {
return out << std::string(1, static_cast<char>(group));
36 changes: 34 additions & 2 deletions astronomy/standard_product_3.hpp
Original file line number Diff line number Diff line change
@@ -7,13 +7,15 @@
#include <vector>

#include "astronomy/frames.hpp"
#include "base/not_null.hpp"
#include "geometry/named_quantities.hpp"
#include "physics/discrete_trajectory.hpp"

namespace principia {
namespace astronomy {
namespace internal_standard_product_3 {

using base::not_null;
using geometry::Instant;
using geometry::Position;
using geometry::Velocity;
@@ -51,6 +53,17 @@ class StandardProduct3 {
// - the minute field of the epoch record takes the value 60 at the end of
// the hour.
ILRSB,
// Files produced by the Groupe de Recherche de Géodesie Spatiale for the
// International DORIS Service.
// Note that these files, while having the letters grg in the file name,
// have the legacy agency field LCA, for Laboratoire d’Études en Géophysique
// et Océanographie Spatiales / Collecte Localisation Satellites (LEGOS/CLS)
// Analysis center.
// The files are syntactically correct, since the nonconformance only
// affects the numerical value of the velocity fields.
// Divergence from the specification:
// - the velocities are in m/s instead of the standard dm/s.
GRGS,
};

enum class SatelliteGroup : char {
@@ -80,13 +93,29 @@ class StandardProduct3 {
// (that order is the same in the satellite ID records and within each epoch).
std::vector<SatelliteIdentifier> const& satellites() const;

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

Version version() const;

// Whether velocities were given in the SP3 file. If this is false, the
// velocities provided by this object are computed from the positions by a
// finite difference formula.
bool file_has_velocities() const;

private:
std::vector<SatelliteIdentifier> satellites_;
std::map<SatelliteIdentifier, DiscreteTrajectory<ITRS>> orbits_;
std::map<SatelliteIdentifier,
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 references to these vectors from
// |StandardProduct3::orbit|.
std::map<SatelliteIdentifier,
std::vector<not_null<DiscreteTrajectory<ITRS> const*>>>
const_orbits_;
Version version_;

bool has_velocities_;
@@ -101,6 +130,9 @@ bool operator<(StandardProduct3::SatelliteIdentifier const& left,
std::ostream& operator<<(std::ostream& out,
StandardProduct3::Version const& version);

std::ostream& operator<<(std::ostream& out,
StandardProduct3::Dialect const& dialect);

std::ostream& operator<<(std::ostream& out,
StandardProduct3::SatelliteGroup const& group);

Loading