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

Commits on Apr 20, 2019

  1. Copy the full SHA
    15560f4 View commit details
  2. Copy the full SHA
    2309c2e View commit details

Commits on Apr 22, 2019

  1. Merge pull request #2140 from pleroy/GeopotentialTest

     Randomized comparison of the Fortran and C++ geopotentials
    pleroy authored Apr 22, 2019
    Copy the full SHA
    97a9607 View commit details
Showing with 74 additions and 37 deletions.
  1. +73 −36 physics/geopotential_test.cpp
  2. +1 −1 physics/solar_system_body.hpp
109 changes: 73 additions & 36 deletions physics/geopotential_test.cpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@

#include "physics/geopotential.hpp"

#include <random>
#include <vector>

#include "astronomy/fortran_astrodynamics_toolkit.hpp"
@@ -34,6 +35,7 @@ using quantities::Angle;
using quantities::AngularFrequency;
using quantities::Degree2SphericalHarmonicCoefficient;
using quantities::Degree3SphericalHarmonicCoefficient;
using quantities::GravitationalParameter;
using quantities::ParseQuantity;
using quantities::Pow;
using quantities::SIUnit;
@@ -71,7 +73,7 @@ class GeopotentialTest : public ::testing::Test {
declination_of_pole_) {}

template<typename Frame>
Vector<Quotient<Acceleration, GravitationalParameter>, Frame>
static Vector<Quotient<Acceleration, GravitationalParameter>, Frame>
SphericalHarmonicsAcceleration(Geopotential<Frame> const& geopotential,
Instant const& t,
Displacement<Frame> const& r) {
@@ -81,7 +83,7 @@ class GeopotentialTest : public ::testing::Test {
}

template<typename Frame>
Vector<Quotient<Acceleration, GravitationalParameter>, Frame>
static Vector<Quotient<Acceleration, GravitationalParameter>, Frame>
GeneralSphericalHarmonicsAcceleration(Geopotential<Frame> const& geopotential,
Instant const& t,
Displacement<Frame> const& r) {
@@ -92,6 +94,41 @@ class GeopotentialTest : public ::testing::Test {
t, r, r_norm, r², one_over_r³);
}

static Vector<Acceleration, ITRS> AccelerationCpp(
Displacement<ITRS> const& displacement,
Geopotential<ICRS> const& geopotential,
OblateBody<ICRS> const& earth) {
Displacement<ICRS> const icrs_displacement =
earth.FromSurfaceFrame<ITRS>(Instant())(displacement);
auto const icrs_acceleration =
earth.gravitational_parameter() *
(GeneralSphericalHarmonicsAcceleration(
geopotential, Instant(), icrs_displacement) -
icrs_displacement / Pow<3>(icrs_displacement.Norm()));
return earth.ToSurfaceFrame<ITRS>(Instant())(icrs_acceleration);
}

static Vector<Acceleration, ITRS> AccelerationF90(
Displacement<ITRS> const& displacement,
OblateBody<ICRS> const& earth) {
double const mu =
earth.gravitational_parameter() / SIUnit<GravitationalParameter>();
double const rbar = earth.reference_radius() / Metre;
numerics::FixedMatrix<double, 10, 10> cnm;
numerics::FixedMatrix<double, 10, 10> snm;
for (int n = 0; n <= 9; ++n) {
for (int m = 0; m <= n; ++m) {
cnm[n][m] = earth.cos()[n][m] * LegendreNormalizationFactor[n][m];
snm[n][m] = earth.sin()[n][m] * LegendreNormalizationFactor[n][m];
}
}
return Vector<Acceleration, ITRS>(
SIUnit<Acceleration>() *
astronomy::fortran_astrodynamics_toolkit::
ComputeGravityAccelerationLear<9, 9>(
displacement.coordinates() / Metre, mu, rbar, cnm, snm));
}

// The axis of rotation is along the z axis for ease of testing.
AngularFrequency const angular_frequency_ = -1.5 * Radian / Second;
Angle const right_ascension_of_pole_ = 0 * Degree;
@@ -300,45 +337,45 @@ TEST_F(GeopotentialTest, TestVector) {
2.3585194144421 * Metre / Second / Second,
-2.9531441870790 * Metre / Second / Second});

Vector<Acceleration, ITRS> actual_acceleration_cpp;
Vector<Acceleration, ITRS> actual_acceleration_f90;
Vector<Acceleration, ITRS> const actual_acceleration_cpp =
AccelerationCpp(displacement, geopotential, earth);
Vector<Acceleration, ITRS> const actual_acceleration_f90 =
AccelerationF90(displacement, earth);

{
Displacement<ICRS> const icrs_displacement =
earth.FromSurfaceFrame<ITRS>(Instant())(
Displacement<ITRS>(displacement));
auto const icrs_acceleration =
earth_μ * (GeneralSphericalHarmonicsAcceleration(
geopotential, Instant(), icrs_displacement) -
icrs_displacement / Pow<3>(icrs_displacement.Norm()));
actual_acceleration_cpp =
earth.ToSurfaceFrame<ITRS>(Instant())(icrs_acceleration);
EXPECT_THAT(RelativeError(actual_acceleration_cpp, expected_acceleration),
IsNear(1.3e-8));
}
{
double mu = earth_μ / SIUnit<GravitationalParameter>();
double rbar = earth_reference_radius / Metre;
numerics::FixedMatrix<double, 10, 10> cnm;
numerics::FixedMatrix<double, 10, 10> snm;
for (int n = 0; n <= 9; ++n) {
for (int m = 0; m <= n; ++m) {
cnm[n][m] = earth.cos()[n][m] * LegendreNormalizationFactor[n][m];
snm[n][m] = earth.sin()[n][m] * LegendreNormalizationFactor[n][m];
}
}
actual_acceleration_f90 = Vector<Acceleration, ITRS>(
SIUnit<Acceleration>() *
astronomy::fortran_astrodynamics_toolkit::
ComputeGravityAccelerationLear<9, 9>(
displacement.coordinates() / Metre, mu, rbar, cnm, snm));
EXPECT_THAT(RelativeError(actual_acceleration_f90, expected_acceleration),
IsNear(1.3e-8));
}
EXPECT_THAT(RelativeError(actual_acceleration_cpp, expected_acceleration),
IsNear(1.3e-8));
EXPECT_THAT(RelativeError(actual_acceleration_f90, expected_acceleration),
IsNear(1.3e-8));
EXPECT_THAT(actual_acceleration_cpp,
AlmostEquals(actual_acceleration_f90, 4, 8));
}

TEST_F(GeopotentialTest, CppF90Comparison) {
SolarSystem<ICRS> solar_system_2000(
SOLUTION_DIR / "astronomy" / "sol_gravity_model.proto.txt",
SOLUTION_DIR / "astronomy" /
"sol_initial_state_jd_2451545_000000000.proto.txt");
solar_system_2000.LimitOblatenessToDegree("Earth", /*max_degree=*/9);
auto earth_message = solar_system_2000.gravity_model_message("Earth");
auto const earth = solar_system_2000.MakeOblateBody(earth_message);
Geopotential<ICRS> const geopotential(earth.get(), /*tolerance=*/0);

std::mt19937_64 random(42);
std::uniform_real_distribution<double> length_distribution(-1e7, 1e7);
for (int i = 0; i < 1000; ++i) {
Displacement<ITRS> const displacement(
{length_distribution(random) * Metre,
length_distribution(random) * Metre,
length_distribution(random) * Metre});
Vector<Acceleration, ITRS> const actual_acceleration_cpp =
AccelerationCpp(displacement, geopotential, *earth);
Vector<Acceleration, ITRS> const actual_acceleration_f90 =
AccelerationF90(displacement, *earth);
EXPECT_THAT(actual_acceleration_cpp,
AlmostEquals(actual_acceleration_f90, 0, 584));
}
}

TEST_F(GeopotentialTest, HarmonicDamping) {
HarmonicDamping σ(1 * Metre);
EXPECT_THAT(σ.inner_threshold(), Eq(1 * Metre));
2 changes: 1 addition & 1 deletion physics/solar_system_body.hpp
Original file line number Diff line number Diff line change
@@ -399,7 +399,7 @@ SolarSystem<Frame>::MakeOblateBody(
serialization::GravityModel::Body const& body) {
Check(body);
CHECK(body.has_mean_radius());
CHECK(body.has_j2());
CHECK(body.has_j2() || body.has_geopotential());
auto const massive_body_parameters = MakeMassiveBodyParameters(body);
auto const rotating_body_parameters = MakeRotatingBodyParameters(body);
auto const oblate_body_parameters = MakeOblateBodyParameters(body);