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: 8e0d73ec77a4
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: f48ea15eeecf
Choose a head ref

Commits on Jul 30, 2019

  1. OrbitalElements

    eggrobin committed Jul 30, 2019
    Copy the full SHA
    12fa64b View commit details
  2. reorder

    eggrobin committed Jul 30, 2019
    Copy the full SHA
    7aab57e View commit details

Commits on Aug 1, 2019

  1. typos

    eggrobin committed Aug 1, 2019
    Copy the full SHA
    3db5139 View commit details
  2. test

    eggrobin committed Aug 1, 2019
    Copy the full SHA
    a273300 View commit details
  3. a test that passes

    eggrobin committed Aug 1, 2019
    Copy the full SHA
    03a8753 View commit details
  4. strange numbers

    eggrobin committed Aug 1, 2019
    Copy the full SHA
    d1a0fb8 View commit details
  5. log

    eggrobin committed Aug 1, 2019
    Copy the full SHA
    bdd7de5 View commit details
  6. ???

    eggrobin committed Aug 1, 2019
    Copy the full SHA
    04c001b View commit details
  7. 64

    eggrobin committed Aug 1, 2019
    Copy the full SHA
    2e3a6b3 View commit details

Commits on Aug 2, 2019

  1. two tests

    eggrobin committed Aug 2, 2019
    Copy the full SHA
    5fb4ae7 View commit details

Commits on Aug 4, 2019

  1. more testing

    eggrobin committed Aug 4, 2019
    Copy the full SHA
    962544b View commit details
  2. renaming

    eggrobin committed Aug 4, 2019
    Copy the full SHA
    6ddefe4 View commit details
  3. remove

    eggrobin committed Aug 4, 2019
    Copy the full SHA
    5a683b3 View commit details

Commits on Aug 5, 2019

  1. unused unwind, line length

    eggrobin committed Aug 5, 2019
    Copy the full SHA
    cfed8ea View commit details
  2. names

    eggrobin committed Aug 5, 2019
    Copy the full SHA
    117d6b9 View commit details

Commits on Aug 6, 2019

  1. iwyu

    eggrobin committed Aug 6, 2019
    Copy the full SHA
    f5b8550 View commit details
  2. ndebug

    eggrobin committed Aug 6, 2019
    Copy the full SHA
    5d523e8 View commit details

Commits on Aug 7, 2019

  1. address some comments

    eggrobin committed Aug 7, 2019
    Copy the full SHA
    fbc55f1 View commit details
  2. NDEBUG

    eggrobin committed Aug 7, 2019
    Copy the full SHA
    4917bbb View commit details

Commits on Aug 8, 2019

  1. comments

    eggrobin committed Aug 8, 2019
    Copy the full SHA
    7a758f2 View commit details
  2. simplify the test

    eggrobin committed Aug 8, 2019
    Copy the full SHA
    59bfb62 View commit details
  3. compactify

    eggrobin committed Aug 8, 2019
    Copy the full SHA
    31b1631 View commit details
  4. more comments

    eggrobin committed Aug 8, 2019
    Copy the full SHA
    e8c8cd1 View commit details
  5. fits on previous line

    eggrobin committed Aug 8, 2019
    Copy the full SHA
    bb8f27e View commit details
  6. line length

    eggrobin committed Aug 8, 2019
    Copy the full SHA
    da46e7c View commit details
  7. blank lines

    eggrobin committed Aug 8, 2019
    Copy the full SHA
    41511ea View commit details

Commits on Aug 9, 2019

  1. ʃ

    eggrobin committed Aug 9, 2019
    Copy the full SHA
    080809c View commit details
  2. Merge pull request #2264 from eggrobin/elements

    Elements
    eggrobin authored Aug 9, 2019
    Copy the full SHA
    f48ea15 View commit details
3 changes: 3 additions & 0 deletions astronomy/astronomy.vcxproj
Original file line number Diff line number Diff line change
@@ -62,6 +62,8 @@
<ClInclude Include="date_time.hpp" />
<ClInclude Include="date_time_body.hpp" />
<ClInclude Include="fortran_astrodynamics_toolkit.hpp" />
<ClInclude Include="orbital_elements.hpp" />
<ClInclude Include="orbital_elements_body.hpp" />
<ClInclude Include="orbit_recurrence.hpp" />
<ClInclude Include="orbit_recurrence_body.hpp" />
<ClInclude Include="solar_system_fingerprints.hpp" />
@@ -87,6 +89,7 @@
<ClCompile Include="ksp_system_test.cpp" />
<ClCompile Include="geodesy_test.cpp" />
<ClCompile Include="lunar_orbit_test.cpp" />
<ClCompile Include="orbital_elements_test.cpp" />
<ClCompile Include="orbit_recurrence_test.cpp" />
<ClCompile Include="solar_system_dynamics_test.cpp" />
<ClCompile Include="standard_product_3.cpp" />
9 changes: 9 additions & 0 deletions astronomy/astronomy.vcxproj.filters
Original file line number Diff line number Diff line change
@@ -219,6 +219,12 @@
<ClInclude Include="orbit_recurrence_body.hpp">
<Filter>Source Files</Filter>
</ClInclude>
<ClInclude Include="orbital_elements_body.hpp">
<Filter>Source Files</Filter>
</ClInclude>
<ClInclude Include="orbital_elements.hpp">
<Filter>Header Files</Filter>
</ClInclude>
</ItemGroup>
<ItemGroup>
<ClCompile Include="lunar_eclipse_test.cpp">
@@ -278,5 +284,8 @@
<ClCompile Include="orbit_recurrence_test.cpp">
<Filter>Test Files</Filter>
</ClCompile>
<ClCompile Include="orbital_elements_test.cpp">
<Filter>Test Files</Filter>
</ClCompile>
</ItemGroup>
</Project>
2 changes: 1 addition & 1 deletion astronomy/frames.hpp
Original file line number Diff line number Diff line change
@@ -60,7 +60,7 @@ using ICRS = Frame<serialization::Frame::SolarSystemTag,
// the BCRS; in particular, the equations of motion in the GCRS include de
// Sitter and Lense–Thirring precession.
// - The time coordinate of the BCRS is TCB, and the time coordinate of the GCRS
// is TCG; TT is a linear scaling of TCG, and TDB is a linear scaling of TDB;
// is TCG; TT is a linear scaling of TCG, and TDB is a linear scaling of TCB;
// for practical purposes TT and TDB are within 2 ms of each other;
// Principia's |Instant| is TT.
using GCRS = Frame<serialization::Frame::SolarSystemTag,
2 changes: 1 addition & 1 deletion astronomy/orbit_recurrence.hpp
Original file line number Diff line number Diff line change
@@ -32,7 +32,7 @@ using quantities::Time;
// definitions of the day:
// — for a sun-synchronous orbit (see sections 7.5.4, 7.5.5), this day is the
// mean solar day;
// — for an orbit with no nodal precession (, i.e., a strictly polar orbit, see
// — for an orbit with no nodal precession (i.e., a strictly polar orbit, see
// section 7.1.3), this day is the stellar day.
// These days are counted negatively for a body with retrograde rotation.
class OrbitRecurrence final {
202 changes: 202 additions & 0 deletions astronomy/orbital_elements.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,202 @@
#pragma once

#include <vector>

#include "base/status_or.hpp"
#include "geometry/named_quantities.hpp"
#include "physics/body.hpp"
#include "physics/discrete_trajectory.hpp"
#include "physics/massive_body.hpp"
#include "quantities/named_quantities.hpp"
#include "quantities/quantities.hpp"

namespace principia {
namespace astronomy {
namespace internal_orbital_elements {

using base::StatusOr;
using geometry::Instant;
using physics::Body;
using physics::DiscreteTrajectory;
using physics::MassiveBody;
using quantities::Angle;
using quantities::AngularFrequency;
using quantities::Difference;
using quantities::Infinity;
using quantities::Length;
using quantities::Time;

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

// The classical Keplerian elements (a, e, i, Ω, ω, M),
// together with an epoch.
struct ClassicalElements {
Instant time;
Length semimajor_axis;
double eccentricity;
Angle inclination;
Angle longitude_of_ascending_node;
Angle argument_of_periapsis;
Angle mean_anomaly;
};

// Mean element time series. These elements are free of short-period
// variations, i.e., variations whose period is the orbital period.
std::vector<ClassicalElements> const& mean_elements() const;

// The period of the (osculating) mean longitude λ = Ω + ω + M.
// Note that since our mean elements are filtered by integration over this
// period, it does not make much sense to recompute it based on our mean
// elements.
Time sidereal_period() const;
// The period of the (mean) mean argument of latitude u = ω + M.
Time nodal_period() const;
// The period of the (mean) mean anomaly M.
Time anomalistic_period() const;

// The rate of precession of Ω.
AngularFrequency nodal_precession() const;

// NOTE(egg): The argument of periapsis ω typically precesses as well.
// However, long-period variations tend to be comparatively large, so that a
// precession rate computed over a few orbits would be highly inaccurate.
// More importantly, whereas the actual value of Ω′ is relevant to, e.g.,
// orbit recurrence computation or sun-synchronicity, one typically cares
// about ω′ only when requiring that ω′ be 0 (in a frozen orbit), in which
// case the more relevant requirement is that ω stays close to some reference
// value.

// Of the mean classical elements (a, e, i, Ω, ω, M), under the influence of
// gravitational forces,
// — M always exhibits a fast secular variation (anomalistic mean motion);
// — Ω often exhibits a secular variation (nodal precession); there are
// however rare cases where it is kept constant (so-called inertial orbits
// that achieve Ω′ = 0 by being polar, e.g., CoRoT or Gravity Probe B); in
// that case, the frozen value may occasionally be relevant: for CoRoT, it
// determines the region of the sky that may be observed.
// — ω exhibits a secular variation, except for frozen orbits or orbits at the
// critical inclination; For frozen orbits (type II frozen orbits in the
// terminology of Ulrich Walter (2018), Astronautics: The Physics of Space
// Flight), its constant value must be either 90° or 270°; for orbits at the
// critical inclination (type I frozen orbits), ω is arbitrary; in highly
// eccentric cases, it is often chosen to be 270° so that the apogee is at
// high latitudes (Молния, みちびき, etc.).
// — a, e, i exhibit no secular variation.
// However, the elements that exhibit no secular variation still have
// long-period variations; instead of trying to characterize these complex
// effects, we provide the interval of values taken by these elements over the
// trajectory being analysed.

// Represents the interval [min, max].
// TODO(egg): This makes sense for T = instant, but InfinitePast and
// InfiniteFuture work differently from Infinity.
template<typename T>
struct Interval {
T min = +Infinity<T>();
T max = -Infinity<T>();

// The Lebesgue measure of this interval.
Difference<T> measure() const;
// The midpoint of this interval; NaN if the interval is empty (min > max).
T midpoint() const;

// Extends this interval so that it contains x.
void Include(T const& x);
};

Interval<Length> mean_semimajor_axis_interval() const;
Interval<double> mean_eccentricity_interval() const;
Interval<Angle> mean_inclination_interval() const;
Interval<Angle> mean_longitude_of_ascending_node_interval() const;
Interval<Angle> mean_argument_of_periapsis_interval() const;

// The equinoctial elements, and in particular the osculating equinoctial
// elements, are not directly interesting; anything that could be derived from
// them should be directly computed by this class instead. They are however
// useful for experimentation in Mathematica, to see whether the
// transformation from osculating to mean elements is well-behaved, whether
// the mean elements are stable, and what useful quantities can be derived
// from the mean elements.

// The equinoctial elements, together with an epoch.
// See Broucke and Cefola (1972), On the equinoctial orbit elements.
struct EquinoctialElements {
Instant t; // The epoch of the elements.
Length a; // The semimajor axis.
double h; // e sin ϖ = e sin (Ω + ω).
double k; // e cos ϖ = e cos (Ω + ω).
Angle λ; // The mean longitude ϖ + M = Ω + ω + M.
double p; // tg i/2 sin Ω.
double q; // tg i/2 cos Ω.
// pʹ and qʹ use the cotangent of the half-inclination instead of its
// tangent; they are better suited to retrograde orbits.
double pʹ; // cotg i/2 sin Ω.
double qʹ; // cotg i/2 cos Ω.
};

std::vector<EquinoctialElements> const& osculating_equinoctial_elements()
const;
std::vector<EquinoctialElements> const& mean_equinoctial_elements() const;

private:
OrbitalElements() = default;

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

// |equinoctial_elements| must contain at least 2 elements.
static 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(
std::vector<EquinoctialElements> const& osculating,
Time const& period);

static 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();

// |mean_classical_elements_| must have been computed; sets
// |mean_*_interval_| accordingly.
void ComputeMeanElementIntervals();

std::vector<EquinoctialElements> osculating_equinoctial_elements_;
Time sidereal_period_;
std::vector<EquinoctialElements> mean_equinoctial_elements_;
std::vector<ClassicalElements> mean_classical_elements_;
Time anomalistic_period_;
Time nodal_period_;
AngularFrequency nodal_precession_;

Interval<Length> mean_semimajor_axis_interval_;
Interval<double> mean_eccentricity_interval_;
Interval<Angle> mean_inclination_interval_;
Interval<Angle> mean_longitude_of_ascending_node_interval_;
Interval<Angle> mean_argument_of_periapsis_interval_;
};

} // namespace internal_orbital_elements

using internal_orbital_elements::OrbitalElements;

} // namespace astronomy
} // namespace principia

#include "astronomy/orbital_elements_body.hpp"
410 changes: 410 additions & 0 deletions astronomy/orbital_elements_body.hpp

Large diffs are not rendered by default.

425 changes: 425 additions & 0 deletions astronomy/orbital_elements_test.cpp

Large diffs are not rendered by default.

6 changes: 5 additions & 1 deletion physics/kepler_orbit.hpp
Original file line number Diff line number Diff line change
@@ -89,7 +89,11 @@ std::ostream& operator<<(std::ostream& out,

template<typename Frame>
class KeplerOrbit final {
static_assert(Frame::is_inertial, "Frame must be inertial");
// TODO(egg): This class used to require:
// static_assert(Frame::is_inertial, "Frame must be inertial");
// However, it only requires a non-rotating frame; for instance, it is often
// used with a body-centred non-rotating frame. Perhaps we should have a
// concept of non-rotating frame.

public:
// Exactly one of the |optional|s must be filled in the given