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: b8d380d8273c
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: 932348f634ed
Choose a head ref
  • 12 commits
  • 16 files changed
  • 1 contributor

Commits on Nov 27, 2020

  1. Copy the full SHA
    3aa692e View commit details
  2. Copy the full SHA
    5208b1d View commit details
  3. flow once

    eggrobin committed Nov 27, 2020
    Copy the full SHA
    3bac1c4 View commit details
  4. more stopping

    eggrobin committed Nov 27, 2020
    Copy the full SHA
    160990c View commit details
  5. more stopping, more moving

    eggrobin committed Nov 27, 2020
    Copy the full SHA
    181db36 View commit details
  6. more moving

    eggrobin committed Nov 27, 2020
    Copy the full SHA
    7b79830 View commit details
  7. Copy the full SHA
    85bdcdf View commit details
  8. move it, move it

    eggrobin committed Nov 27, 2020
    Copy the full SHA
    21f5563 View commit details

Commits on Nov 28, 2020

  1. try to make it build

    eggrobin committed Nov 28, 2020
    Copy the full SHA
    eb21510 View commit details
  2. lint

    eggrobin committed Nov 28, 2020
    Copy the full SHA
    66fff58 View commit details
  3. fancier RETURN_IF_ERROR

    eggrobin committed Nov 28, 2020
    Copy the full SHA
    aba47b9 View commit details
  4. Merge pull request #2799 from eggrobin/orbit-analyser

    Interrupt the orbit analysis if the mission duration changes
    eggrobin authored Nov 28, 2020

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    932348f View commit details
14 changes: 7 additions & 7 deletions astronomy/orbit_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -172,10 +172,10 @@ class OrbitAnalysisTest : public ::testing::Test {
std::tuple<OrbitalElements, OrbitRecurrence, OrbitGroundTrack>
ElementsAndRecurrence(SP3Orbit const& orbit) {
auto const earth_centred_trajectory = EarthCentredTrajectory(orbit);
auto const elements = OrbitalElements::ForTrajectory(
*earth_centred_trajectory,
earth_,
MasslessBody{}).ValueOrDie();
auto elements = OrbitalElements::ForTrajectory(
*earth_centred_trajectory,
earth_,
MasslessBody{}).ValueOrDie();

{
auto const identifier = (std::stringstream() << orbit.satellite).str();
@@ -221,14 +221,14 @@ class OrbitAnalysisTest : public ::testing::Test {
129'602'768.13 * ArcSecond / (100 * JulianYear),
1.089 * ArcSecond / Pow<2>(100 * JulianYear)},
"1899-12-31T12:00:00"_TT);
auto const ground_track = OrbitGroundTrack::ForTrajectory(
auto ground_track = OrbitGroundTrack::ForTrajectory(
*earth_centred_trajectory,
earth_,
{{/*epoch=*/J2000,
/*mean_longitude_at_epoch=*/newcomb_mean_longitude(J2000),
/*year*/ 2 * π * Radian /
newcomb_mean_longitude.EvaluateDerivative(J2000)}});
return {elements, recurrence, ground_track};
newcomb_mean_longitude.EvaluateDerivative(J2000)}}).ValueOrDie();
return {std::move(elements), recurrence, std::move(ground_track)};
}

SolarSystem<ICRS> earth_1957_;
10 changes: 9 additions & 1 deletion astronomy/orbit_ground_track.hpp
Original file line number Diff line number Diff line change
@@ -5,6 +5,7 @@

#include "astronomy/orbit_recurrence.hpp"
#include "astronomy/orbital_elements.hpp"
#include "base/status_or.hpp"
#include "geometry/interval.hpp"
#include "physics/discrete_trajectory.hpp"
#include "physics/rotating_body.hpp"
@@ -13,6 +14,7 @@ namespace principia {
namespace astronomy {
namespace internal_orbit_ground_track {

using base::StatusOr;
using geometry::Instant;
using geometry::Interval;
using physics::DiscreteTrajectory;
@@ -55,11 +57,17 @@ class OrbitGroundTrack {
friend class OrbitGroundTrack;
};

OrbitGroundTrack(OrbitGroundTrack const&) = delete;
OrbitGroundTrack(OrbitGroundTrack&&) = default;

OrbitGroundTrack& operator=(OrbitGroundTrack const&) = delete;
OrbitGroundTrack& operator=(OrbitGroundTrack&&) = default;

// Returns an object that describes the properties of the ground track of
// |trajectory| as an orbit around |primary|; if |mean_sun| is provided,
// sun-synchronicity is analysed.
template<typename PrimaryCentred, typename Inertial>
static OrbitGroundTrack ForTrajectory(
static StatusOr<OrbitGroundTrack> ForTrajectory(
DiscreteTrajectory<PrimaryCentred> const& trajectory,
RotatingBody<Inertial> const& primary,
std::optional<MeanSun> const& mean_sun);
14 changes: 7 additions & 7 deletions astronomy/orbit_ground_track_body.hpp
Original file line number Diff line number Diff line change
@@ -118,19 +118,19 @@ inline OrbitGroundTrack::EquatorCrossingLongitudes::EquatorCrossingLongitudes(
: nominal_recurrence_(nominal_recurrence) {}

template<typename PrimaryCentred, typename Inertial>
OrbitGroundTrack OrbitGroundTrack::ForTrajectory(
StatusOr<OrbitGroundTrack> OrbitGroundTrack::ForTrajectory(
DiscreteTrajectory<PrimaryCentred> const& trajectory,
RotatingBody<Inertial> const& primary,
std::optional<MeanSun> const& mean_sun) {
DiscreteTrajectory<PrimaryCentred> ascending_nodes;
DiscreteTrajectory<PrimaryCentred> descending_nodes;
OrbitGroundTrack ground_track;
ComputeNodes(trajectory.begin(),
trajectory.end(),
Vector<double, PrimaryCentred>({0, 0, 1}),
/*max_points=*/std::numeric_limits<int>::max(),
ascending_nodes,
descending_nodes);
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()) {
ground_track.mean_solar_times_of_ascending_nodes_ =
17 changes: 12 additions & 5 deletions astronomy/orbital_elements.hpp
Original file line number Diff line number Diff line change
@@ -15,6 +15,7 @@ namespace principia {
namespace astronomy {
namespace internal_orbital_elements {

using base::Status;
using base::StatusOr;
using geometry::Instant;
using geometry::Interval;
@@ -30,6 +31,12 @@ using quantities::Time;

class OrbitalElements {
public:
OrbitalElements(OrbitalElements const&) = delete;
OrbitalElements(OrbitalElements&&) = default;

OrbitalElements& operator=(OrbitalElements const&) = delete;
OrbitalElements& operator=(OrbitalElements&&) = default;

template<typename PrimaryCentred>
static StatusOr<OrbitalElements> ForTrajectory(
DiscreteTrajectory<PrimaryCentred> const& trajectory,
@@ -152,28 +159,28 @@ class OrbitalElements {
DiscreteTrajectory<PrimaryCentred> const& trajectory);

// |equinoctial_elements| must contain at least 2 elements.
static Time SiderealPeriod(
static StatusOr<Time> SiderealPeriod(
std::vector<EquinoctialElements> const& equinoctial_elements);

// |osculating| must contain at least 2 elements.
// The resulting elements are averaged over one period, centred on
// their |EquinoctialElements::t|.
static std::vector<EquinoctialElements> MeanEquinoctialElements(
static StatusOr<std::vector<EquinoctialElements>> MeanEquinoctialElements(
std::vector<EquinoctialElements> const& osculating,
Time const& period);

static std::vector<ClassicalElements> ToClassicalElements(
static StatusOr<std::vector<ClassicalElements>> ToClassicalElements(
std::vector<EquinoctialElements> const& equinoctial_elements);

// |mean_classical_elements_| must have been computed; sets
// |anomalistic_period_|, |nodal_period_|, and |nodal_precession_|
// accordingly. Note that this does not compute |sidereal_period_| (our mean
// element computation is based on it, so it gets computed earlier).
void ComputePeriodsAndPrecession();
Status ComputePeriodsAndPrecession();

// |radial_distances_| and |mean_classical_elements_| must have been computed;
// sets |radial_distance_interval_| and |mean_*_interval_| accordingly.
void ComputeIntervals();
Status ComputeIntervals();

std::vector<EquinoctialElements> osculating_equinoctial_elements_;
std::vector<Length> radial_distances_;
39 changes: 29 additions & 10 deletions astronomy/orbital_elements_body.hpp
Original file line number Diff line number Diff line change
@@ -5,6 +5,7 @@
#include <algorithm>
#include <vector>

#include "base/jthread.hpp"
#include "physics/kepler_orbit.hpp"
#include "quantities/elementary_functions.hpp"

@@ -14,6 +15,7 @@ namespace internal_orbital_elements {

using base::Error;
using base::Status;
using base::this_stoppable_thread;
using geometry::Velocity;
using physics::DegreesOfFreedom;
using physics::KeplerOrbit;
@@ -42,17 +44,22 @@ StatusOr<OrbitalElements> OrbitalElements::ForTrajectory(
orbital_elements.osculating_equinoctial_elements_ =
OsculatingEquinoctialElements(trajectory, primary, secondary);
orbital_elements.radial_distances_ = RadialDistances(trajectory);
orbital_elements.sidereal_period_ =
auto const sidereal_period =
SiderealPeriod(orbital_elements.osculating_equinoctial_elements_);
RETURN_IF_ERROR(sidereal_period);
orbital_elements.sidereal_period_ = sidereal_period.ValueOrDie();
if (!IsFinite(orbital_elements.sidereal_period_)) {
// Guard against NaN sidereal periods (from hyperbolic orbits).
return Status(
Error::OUT_OF_RANGE,
"sidereal period is " + DebugString(orbital_elements.sidereal_period_));
}
orbital_elements.mean_equinoctial_elements_ =
auto mean_equinoctial_elements =
MeanEquinoctialElements(orbital_elements.osculating_equinoctial_elements_,
orbital_elements.sidereal_period_);
RETURN_IF_ERROR(mean_equinoctial_elements);
orbital_elements.mean_equinoctial_elements_ =
std::move(mean_equinoctial_elements).ValueOrDie();
if (orbital_elements.mean_equinoctial_elements_.size() < 2) {
return Status(
Error::OUT_OF_RANGE,
@@ -62,10 +69,13 @@ StatusOr<OrbitalElements> OrbitalElements::ForTrajectory(
DebugString(trajectory.back().time -
trajectory.front().time));
}
orbital_elements.mean_classical_elements_ =
auto mean_classical_elements =
ToClassicalElements(orbital_elements.mean_equinoctial_elements_);
orbital_elements.ComputePeriodsAndPrecession();
orbital_elements.ComputeIntervals();
RETURN_IF_ERROR(mean_classical_elements);
orbital_elements.mean_classical_elements_ =
std::move(mean_classical_elements).ValueOrDie();
RETURN_IF_ERROR(orbital_elements.ComputePeriodsAndPrecession());
RETURN_IF_ERROR(orbital_elements.ComputeIntervals());
return orbital_elements;
}

@@ -187,7 +197,7 @@ OrbitalElements::mean_equinoctial_elements() const {
return mean_equinoctial_elements_;
}

inline Time OrbitalElements::SiderealPeriod(
inline StatusOr<Time> OrbitalElements::SiderealPeriod(
std::vector<EquinoctialElements> const& equinoctial_elements) {
Time const Δt =
equinoctial_elements.back().t - equinoctial_elements.front().t;
@@ -200,6 +210,7 @@ inline Time OrbitalElements::SiderealPeriod(
for (auto previous = first, it = first + 1;
it != equinoctial_elements.end();
previous = it, ++it) {
RETURN_IF_STOPPED;
Time const t = it->t - t0;
auto const λt = it->λ * t;
Time const dt = it->t - previous->t;
@@ -209,7 +220,7 @@ inline Time OrbitalElements::SiderealPeriod(
return 2 * π * Radian * Pow<3>(Δt) / (12 * ʃ_λt_dt);
}

inline std::vector<OrbitalElements::EquinoctialElements>
inline StatusOr<std::vector<OrbitalElements::EquinoctialElements>>
OrbitalElements::MeanEquinoctialElements(
std::vector<EquinoctialElements> const& osculating,
Time const& period) {
@@ -249,6 +260,7 @@ OrbitalElements::MeanEquinoctialElements(
for (auto previous = osculating.begin(), it = osculating.begin() + 1;
it != osculating.end();
previous = it, ++it) {
RETURN_IF_STOPPED;
integrals.push_back(integrals.back());
integrals.back().t_max = it->t;
Time const dt = it->t - previous->t;
@@ -266,6 +278,7 @@ OrbitalElements::MeanEquinoctialElements(
std::vector<EquinoctialElements> mean_elements;
int j = 0;
for (auto const& up_to_tᵢ : integrals) {
RETURN_IF_STOPPED;
// We are averaging the elements over the interval [tᵢ, tᵢ + period].
Instant const tᵢ = up_to_tᵢ.t_max;
while (integrals[j].t_max < tᵢ + period) {
@@ -318,11 +331,12 @@ OrbitalElements::MeanEquinoctialElements(
return mean_elements;
}

inline std::vector<OrbitalElements::ClassicalElements>
inline StatusOr<std::vector<OrbitalElements::ClassicalElements>>
OrbitalElements::ToClassicalElements(
std::vector<EquinoctialElements> const& equinoctial_elements) {
std::vector<ClassicalElements> classical_elements;
for (auto const& equinoctial : equinoctial_elements) {
RETURN_IF_STOPPED;
double const tg_½i = Sqrt(Pow<2>(equinoctial.p) + Pow<2>(equinoctial.q));
double const cotg_½i =
Sqrt(Pow<2>(equinoctial.pʹ) + Pow<2>(equinoctial.qʹ));
@@ -355,7 +369,7 @@ OrbitalElements::ToClassicalElements(
return classical_elements;
}

inline void OrbitalElements::ComputePeriodsAndPrecession() {
inline Status OrbitalElements::ComputePeriodsAndPrecession() {
Time const Δt = mean_classical_elements_.back().time -
mean_classical_elements_.front().time;
auto const Δt³ = Pow<3>(Δt);
@@ -384,6 +398,7 @@ inline void OrbitalElements::ComputePeriodsAndPrecession() {
for (auto previous = first, it = first + 1;
it != mean_classical_elements_.end();
previous = it, ++it) {
RETURN_IF_STOPPED;
Angle const u = it->argument_of_periapsis + it->mean_anomaly;
Angle const& M = it->mean_anomaly;
Angle const& Ω = it->longitude_of_ascending_node;
@@ -406,13 +421,16 @@ inline void OrbitalElements::ComputePeriodsAndPrecession() {
anomalistic_period_ = 2 * π * Radian * Δt³ / (12 * ʃ_Mt_dt);
nodal_period_ = 2 * π * Radian * Δt³ / (12 * ʃ_ut_dt);
nodal_precession_ = 12 * ʃ_Ωt_dt / Δt³;
return Status::OK;
}

inline void OrbitalElements::ComputeIntervals() {
inline Status OrbitalElements::ComputeIntervals() {
for (auto const& r : radial_distances_) {
RETURN_IF_STOPPED;
radial_distance_interval_.Include(r);
}
for (auto const& elements : mean_classical_elements_) {
RETURN_IF_STOPPED;
mean_semimajor_axis_interval_.Include(elements.semimajor_axis);
mean_eccentricity_interval_.Include(elements.eccentricity);
mean_inclination_interval_.Include(elements.inclination);
@@ -423,6 +441,7 @@ inline void OrbitalElements::ComputeIntervals() {
mean_periapsis_distance_interval_.Include(elements.periapsis_distance);
mean_apoapsis_distance_interval_.Include(elements.apoapsis_distance);
}
return Status::OK;
}


10 changes: 10 additions & 0 deletions base/jthread.hpp
Original file line number Diff line number Diff line change
@@ -7,6 +7,7 @@

#include "absl/synchronization/mutex.h"
#include "base/not_null.hpp"
#include "base/status.hpp"

namespace principia {
namespace base {
@@ -123,6 +124,15 @@ class this_stoppable_thread {
friend jthread MakeStoppableThread(Function&& f, Args&&... args);
};

#define RETURN_IF_STOPPED \
do { \
if (::principia::base::this_stoppable_thread::get_stop_token() \
.stop_requested()) { \
return ::principia::base::Status(::principia::base::Error::CANCELLED, \
"Cancelled by stop token"); \
} \
} while (false)

} // namespace internal_jthread

using internal_jthread::MakeStoppableThread;
6 changes: 6 additions & 0 deletions base/jthread_body.hpp
Original file line number Diff line number Diff line change
@@ -118,6 +118,12 @@ inline jthread::jthread(jthread&& other)
thread_(std::move(other.thread_)) {}

inline jthread& jthread::operator=(jthread&& other) {
if (stop_state_ != nullptr) {
stop_state_->request_stop();
}
if (thread_.joinable()) {
thread_.join();
}
stop_state_ = std::move(other.stop_state_);
thread_ = std::move(other.thread_);
return *this;
7 changes: 6 additions & 1 deletion base/status.hpp
Original file line number Diff line number Diff line change
@@ -103,6 +103,10 @@ class Status final {
std::string message_;
};

inline Status const& GetStatus(Status const& s) {
return s;
}

// Prints a human-readable representation of |s| to |os|.
std::ostream& operator<<(std::ostream& os, Status const& s);

@@ -113,7 +117,8 @@ std::ostream& operator<<(std::ostream& os, Status const& s);
#define RETURN_IF_ERROR(expr) \
do { \
/* Using _status below to avoid capture problems if expr is "status". */ \
::principia::base::Status const _status = (expr); \
::principia::base::Status const _status = \
(::principia::base::GetStatus(expr)); \
if (!_status.ok()) \
return _status; \
} while (false)
Loading