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

Commits on Aug 22, 2021

  1. logger

    eggrobin committed Aug 22, 2021
    Copy the full SHA
    66b10c4 View commit details

Commits on Aug 23, 2021

  1. add the files…

    eggrobin committed Aug 23, 2021
    Copy the full SHA
    cc352ad View commit details
  2. citations

    eggrobin committed Aug 23, 2021
    Copy the full SHA
    14daba5 View commit details
  3. not in debug

    eggrobin committed Aug 23, 2021
    Copy the full SHA
    dc920ee View commit details
  4. grammar

    eggrobin committed Aug 23, 2021
    Copy the full SHA
    a1de6f7 View commit details
  5. Copy the full SHA
    a891c26 View commit details

Commits on Aug 24, 2021

  1. after pleroy’s review

    eggrobin committed Aug 24, 2021
    Copy the full SHA
    97fe163 View commit details
  2. add the renamed file

    eggrobin committed Aug 24, 2021
    Copy the full SHA
    f7bc98a View commit details
  3. Merge pull request #3113 from eggrobin/Лидов古在

    Лидов–古在
    eggrobin authored Aug 24, 2021
    Copy the full SHA
    07e8049 View commit details
2 changes: 2 additions & 0 deletions astronomy/astronomy.vcxproj
Original file line number Diff line number Diff line change
@@ -100,6 +100,7 @@
<ClInclude Include="experimental_eop_c02.generated.h" />
<ClInclude Include="frames.hpp" />
<ClInclude Include="fortran_astrodynamics_toolkit_body.hpp" />
<ClInclude Include="mercury_orbiter.hpp" />
</ItemGroup>
<ItemGroup>
<ClCompile Include="..\base\bundle.cpp" />
@@ -123,6 +124,7 @@
<ClCompile Include="lunar_eclipse_test.cpp" />
<ClCompile Include="mercury_perihelion_test.cpp" />
<ClCompile Include="trappist_dynamics_test.cpp" />
<ClCompile Include="лидов_古在_test.cpp" />
<ClCompile Include="молния_orbit_test.cpp" />
</ItemGroup>
</Project>
6 changes: 6 additions & 0 deletions astronomy/astronomy.vcxproj.filters
Original file line number Diff line number Diff line change
@@ -267,6 +267,9 @@
<ClInclude Include="orbit_ground_track_body.hpp">
<Filter>Source Files</Filter>
</ClInclude>
<ClInclude Include="mercury_orbiter.hpp">
<Filter>Header Files</Filter>
</ClInclude>
</ItemGroup>
<ItemGroup>
<ClCompile Include="lunar_eclipse_test.cpp">
@@ -335,5 +338,8 @@
<ClCompile Include="..\geometry\instant_output.cpp">
<Filter>Source Files</Filter>
</ClCompile>
<ClCompile Include="лидов_古在_test.cpp">
<Filter>Test Files</Filter>
</ClCompile>
</ItemGroup>
</Project>
44 changes: 44 additions & 0 deletions astronomy/mercury_orbiter.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@
#pragma once

#include "astronomy/time_scales.hpp"
#include "geometry/named_quantities.hpp"
#include "physics/degrees_of_freedom.hpp"

namespace principia {
namespace astronomy {
namespace internal_mercury_orbiter {

using astronomy::operator""_TT;
using geometry::Displacement;
using geometry::Velocity;
using geometry::Instant;
using physics::DegreesOfFreedom;
using quantities::si::Metre;
using quantities::si::Second;

// State of the spacecraft Mercury Orbiter 1 at the start of September, from a
// save by Butcher given in issue #1119. This is used both from a plugin
// compatibility test which reads that save, and from an astronomical test which
// looks for the Лидов–古在 mechanism in the orbit.

constexpr Instant MercuryOrbiterInitialTime =
"1966-09-01T00:16:55"_TT + 0.2571494579315186 * Second;

template<typename Barycentric>
physics::DegreesOfFreedom<Barycentric> const
MercuryOrbiterInitialDegreesOfFreedom = {
Barycentric::origin +
Displacement<Barycentric>({-2.40627773705000000e+10 * Metre,
+3.52445087251250000e+10 * Metre,
+2.13640458684375000e+10 * Metre}),
Velocity<Barycentric>({-5.19594188203811646e+04 * (Metre / Second),
-2.23741500134468079e+04 * (Metre / Second),
-7.15344990825653076e+03 * (Metre / Second)})};

} // namespace internal_mercury_orbiter

using internal_mercury_orbiter::MercuryOrbiterInitialDegreesOfFreedom;
using internal_mercury_orbiter::MercuryOrbiterInitialTime;

} // namespace astronomy
} // namespace principia
178 changes: 178 additions & 0 deletions astronomy/лидов_古在_test.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,178 @@

#include <memory>

#include "astronomy/mercury_orbiter.hpp"
#include "astronomy/orbital_elements.hpp"
#include "astronomy/frames.hpp"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "physics/body_centred_non_rotating_dynamic_frame.hpp"
#include "physics/solar_system.hpp"
#include "mathematica/mathematica.hpp"
#include "testing_utilities/matchers.hpp"
#include "testing_utilities/is_near.hpp"

namespace principia {
namespace astronomy {

using base::not_null;
using geometry::Frame;
using geometry::Instant;
using geometry::Interval;
using geometry::NonRotating;
using geometry::Position;
using integrators::SymmetricLinearMultistepIntegrator;
using integrators::methods::Quinlan1999Order8A;
using integrators::methods::QuinlanTremaine1990Order12;
using physics::BodyCentredNonRotatingDynamicFrame;
using physics::DiscreteTrajectory;
using physics::Ephemeris;
using physics::MassiveBody;
using physics::MasslessBody;
using physics::SolarSystem;
using physics::Trajectory;
using quantities::Cos;
using quantities::Length;
using quantities::Pow;
using quantities::Sin;
using quantities::si::Degree;
using quantities::si::Kilo;
using quantities::si::Metre;
using quantities::si::Milli;
using quantities::si::Minute;
using quantities::si::Radian;
using quantities::si::Second;
using testing_utilities::operator""_⑴;
using testing_utilities::IsNear;

constexpr std::int64_t MaxDenseIntervals = 10'000;
constexpr Length DownsamplingTolerance = 10 * Metre;

// A test that showcases the eccentricity-inclination exchange mechanism
// described in [Лид61] and [Koz62]. We follow the treatment in [Лид61].
class Лидов古在Test : public ::testing::Test {
protected:
using MercuryCentredInertial =
Frame<enum class MercuryCentredInertialTag, NonRotating>;

Лидов古在Test()
: solar_system_1950_(
SOLUTION_DIR / "astronomy" / "sol_gravity_model.proto.txt",
SOLUTION_DIR / "astronomy" /
"sol_initial_state_jd_2433282_500000000.proto.txt"),
ephemeris_(solar_system_1950_.MakeEphemeris(
/*accuracy_parameters=*/{/*fitting_tolerance=*/1 * Milli(Metre),
/*geopotential_tolerance=*/0x1p-24},
Ephemeris<ICRS>::FixedStepParameters(
SymmetricLinearMultistepIntegrator<QuinlanTremaine1990Order12,
Position<ICRS>>(),
/*step=*/10 * Minute))),
mercury_(*solar_system_1950_.massive_body(*ephemeris_, "Mercury")),
mercury_frame_(ephemeris_.get(), &mercury_) {
google::SetStderrLogging(google::INFO);
}

~Лидов古在Test() override {
google::SetStderrLogging(FLAGS_stderrthreshold);
}

SolarSystem<ICRS> solar_system_1950_;
not_null<std::unique_ptr<Ephemeris<ICRS>>> ephemeris_;
MassiveBody const& mercury_;
BodyCentredNonRotatingDynamicFrame<ICRS, MercuryCentredInertial>
mercury_frame_;
};

#if !_DEBUG
TEST_F(Лидов古在Test, MercuryOrbiter) {
DiscreteTrajectory<ICRS> icrs_trajectory;
icrs_trajectory.Append(MercuryOrbiterInitialTime,
MercuryOrbiterInitialDegreesOfFreedom<ICRS>);
icrs_trajectory.SetDownsampling(MaxDenseIntervals, DownsamplingTolerance);
auto const instance = ephemeris_->NewInstance(
{&icrs_trajectory},
Ephemeris<ICRS>::NoIntrinsicAccelerations,
{SymmetricLinearMultistepIntegrator<Quinlan1999Order8A, Position<ICRS>>(),
/*step=*/10 * Second});

for (int year = 1967; year < 1980; ++year) {
Instant const t = DateTimeAsTT(date_time::DateTime::BeginningOfDay(
date_time::Date::Ordinal(year, 1)));
LOG(INFO) << "Flowing to " << t;
auto const status = ephemeris_->FlowWithFixedStep(t, *instance);
if (!status.ok()) {
LOG(INFO) << status << " at " << icrs_trajectory.back().time;
break;
}
}
mathematica::Logger logger(
SOLUTION_DIR / "mathematica" / u"лидов_古在.generated.wl",
/*make_unique=*/false);

DiscreteTrajectory<MercuryCentredInertial> mercury_centred_trajectory;
for (auto const& [t, dof] : icrs_trajectory) {
mercury_centred_trajectory.Append(t,
mercury_frame_.ToThisFrameAtTime(t)(dof));
logger.Append(
"q",
mercury_centred_trajectory.back().degrees_of_freedom.position(),
mathematica::ExpressIn(Metre));
}

OrbitalElements const elements = OrbitalElements::ForTrajectory(
mercury_centred_trajectory, mercury_, MasslessBody{}).value();
// The constants c₁ and c₂ are defined in [Лид61], equations (58) and (59)
// respectively.
Interval<double> c₁;
Interval<double> c₂;
for (auto const elements : elements.mean_elements()) {
double const ε = 1 - Pow<2>(elements.eccentricity);
double const cos²_i = Pow<2>(Cos(elements.inclination));
double const sin²_i = Pow<2>(Sin(elements.inclination));
double const sin²_ω = Pow<2>(Sin(elements.argument_of_periapsis));
c₁.Include(ε * cos²_i);
c₂.Include((1 - ε) * (2.0 / 5.0 - sin²_i * sin²_ω));
logger.Append("t", elements.time, mathematica::ExpressIn(Second));
logger.Append("a", elements.semimajor_axis, mathematica::ExpressIn(Metre));
logger.Append("e", elements.eccentricity);
logger.Append("i", elements.inclination, mathematica::ExpressIn(Radian));
logger.Append(R"(\[Omega])",
elements.argument_of_periapsis,
mathematica::ExpressIn(Radian));
}
// The elements e, i, and ω all vary quite a lot.
EXPECT_THAT(elements.mean_eccentricity_interval().min, IsNear(0.40_⑴));
EXPECT_THAT(elements.mean_eccentricity_interval().max, IsNear(0.88_⑴));
EXPECT_THAT(elements.mean_inclination_interval().min, IsNear(62_⑴ * Degree));
EXPECT_THAT(elements.mean_inclination_interval().max, IsNear(77_⑴ * Degree));
EXPECT_THAT(elements.mean_argument_of_periapsis_interval().min,
IsNear(51_⑴ * Degree));
EXPECT_THAT(elements.mean_argument_of_periapsis_interval().max,
IsNear(129_⑴ * Degree));

// The conservation of the “тривиального интеграла a = const” [Лид61, p. 25]
// is excellent: while the sun is nudging and deforming the orbit, it is not
// pumping energy into nor out of it.
EXPECT_THAT(elements.mean_semimajor_axis_interval().min,
IsNear(14'910.0_⑴ * Kilo(Metre)));
EXPECT_THAT(elements.mean_semimajor_axis_interval().max,
IsNear(14'910.3_⑴ * Kilo(Metre)));

// The integral c₁ is preserved quite well: we have an exchange between
// inclination and eccentricity.
EXPECT_THAT(c₁.min, IsNear(0.042_⑴));
EXPECT_THAT(c₁.max, IsNear(0.050_⑴));

// The integral c₂ is also conserved: the long-term evolution of the orbital
// elements is as described in [Лид61].
EXPECT_THAT(c₂.min, IsNear(-0.091_⑴));
EXPECT_THAT(c₂.max, IsNear(-0.083_⑴));

// TODO(egg): The above are integrals of motion only when averaging over an
// orbit of the perturbing body (so here, over the orbit of Mercury); see what
// things look like under a moving average.
}
#endif

} // namespace astronomy
} // namespace principia
34 changes: 34 additions & 0 deletions documentation/bibliography.bib
Original file line number Diff line number Diff line change
@@ -816,6 +816,17 @@ @article{IAUWGCC2009
volume = {109},
}

@article{Kozai1962,
author = {Kozai, Yoshihide},
date = {1962-11},
doi = {10.1086/108790},
journaltitle = {The Astronomical Journal},
number = {9},
pages = {591--598},
title = {Secular Perturbations of Asteroids with High Inclination and Eccentricity},
volume = {67},
}

@article{Kudryavtsev2007,
author = {Kudryavtsev, S. M.},
url = {http://dx.doi.org/10.1051/0004-6361:20077568},
@@ -861,6 +872,18 @@ @article{Laplace1772
volume = {1792},
}

@article{LidovCleaves1962,
author = {Lidov., M. L. and Cleaves, H. F.},
date = {1961-10},
doi = {10.1016/0032-0633(62)90129-0},
journaltitle = {Planetary and Space Science},
note = {Translated from the Russian \cite{Лидов1961}},
number = {10},
pages = {719--759},
title = {The evolution of orbits of artificial satellites of planets under the action of gravitational perturbations of external bodies},
volume = {9},
}

@article{LǐLǐLǐ2012,
author = {\chinese{李松明} and \chinese{李岩} and \chinese{李劲东}},
shortauthor = {Lǐ, Sōng Míng and Lǐ, Yán and Lǐ, Jìn Dōng},
@@ -1355,6 +1378,17 @@ @article{ZhōuLǐXìnWéiZhāng2012
volume = {29},
}

@article{Лидов1961,
author = {Лидов, М. Л.},
publisher = {Академия Наук СССР},
addendum = {English translation: \cite{LidovCleaves1962}},
date = {1961},
journaltitle = {Искусственные Спутники Земли},
number = {8},
pages = {5--45},
title = {Эволюция орбит искусственных спутников планет под действием гравитационных возмущений внешних тел},
}

@book{Babbage1827,
author = {Babbage, Charles},
publisher = {J. Mawman},
Binary file modified documentation/bibliography.pdf
Binary file not shown.
21 changes: 6 additions & 15 deletions ksp_plugin_test/plugin_compatibility_test.cpp
Original file line number Diff line number Diff line change
@@ -5,6 +5,7 @@
#include <vector>

#include "astronomy/time_scales.hpp"
#include "astronomy/mercury_orbiter.hpp"
#include "base/file.hpp"
#include "base/not_null.hpp"
#include "base/pull_serializer.hpp"
@@ -22,6 +23,8 @@ namespace principia {
namespace interface {

using astronomy::operator""_TT;
using astronomy::MercuryOrbiterInitialDegreesOfFreedom;
using astronomy::MercuryOrbiterInitialTime;
using astronomy::TTSecond;
using astronomy::date_time::DateTime;
using astronomy::date_time::operator""_DateTime;
@@ -333,24 +336,12 @@ TEST_F(PluginCompatibilityTest, DISABLED_Butcher) {
-7.76112133789062500e+03 * (Metre / Second)}))));

// We arrive in late August. Check the state in the beginning of September.
// TODO(egg): Move the expected values for this initial state to a header in
// astronomy, and use that to look for the Лидов–古在 mechanism.
auto const it = orbiter.psychohistory().LowerBound("1966-09-01T00:00:00"_TT);
EXPECT_THAT(it->time,
Eq("1966-09-01T00:16:55"_TT + 0.2571494579315186 * Second));
EXPECT_THAT(it->time, Eq(MercuryOrbiterInitialTime));
EXPECT_THAT(it->degrees_of_freedom,
Eq(DegreesOfFreedom<Barycentric>(
Barycentric::origin + Displacement<Barycentric>(
{-2.40627773705000000e+10 * Metre,
+3.52445087251250000e+10 * Metre,
+2.13640458684375000e+10 * Metre}),
Velocity<Barycentric>(
{-5.19594188203811646e+04 * (Metre / Second),
-2.23741500134468079e+04 * (Metre / Second),
-7.15344990825653076e+03 * (Metre / Second)}))));
Eq(MercuryOrbiterInitialDegreesOfFreedom<Barycentric>));
EXPECT_THAT((it->degrees_of_freedom.position() -
mercury.trajectory().EvaluatePosition(it->time))
.Norm(),
mercury.trajectory().EvaluatePosition(it->time)).Norm(),
IsNear(19'163_⑴ * Kilo(Metre)));

// Make sure that we can upgrade, save, and reload.