diff --git a/ksp_plugin/interface_external.cpp b/ksp_plugin/interface_external.cpp index 81e97ecbec..6cfadcae0e 100644 --- a/ksp_plugin/interface_external.cpp +++ b/ksp_plugin/interface_external.cpp @@ -3,6 +3,7 @@ #include +#include "absl/strings/str_cat.h" #include "base/array.hpp" #include "base/status.hpp" #include "base/status_or.hpp" @@ -23,6 +24,7 @@ using ksp_plugin::Vessel; using physics::BodyCentredNonRotatingDynamicFrame; using physics::ComputeApsides; using physics::DiscreteTrajectory; +using physics::OblateBody; using physics::RigidMotion; using physics::RigidTransformation; @@ -71,6 +73,79 @@ Status principia__ExternalFlowFreefall( "|ExternalFlowFreefall| is not yet implemented")); } +Status principia__ExternalGeopotentialGetCoefficient( + Plugin const* const plugin, + int const body_index, + int const degree, + int const order, + XY* const coefficient) { + journal::Method m{ + {plugin, + body_index, + degree, + order}, + {coefficient}}; + if (plugin == nullptr) { + return m.Return( + MakeStatus(Error::INVALID_ARGUMENT, "|plugin| must not be null")); + } + if (!plugin->HasCelestial(body_index)) { + return m.Return(MakeStatus( + Error::NOT_FOUND, + absl::StrCat("No celestial with index ", body_index))); + } + if (order < 0 || order > degree) { + return m.Return(MakeStatus( + Error::INVALID_ARGUMENT, + absl::StrCat(u8"Expected 0 ≤ order ≤ degree; got degree = ", + degree, ", order = ", order))); + } + if (degree == 0) { + *coefficient = {1, 0}; + return m.Return(OK()); + } + auto const& body = *plugin->GetCelestial(body_index).body(); + if (!body.is_oblate()) { + *coefficient = {0, 0}; + return m.Return(OK()); + } + auto const& oblate_body = dynamic_cast const&>(body); + if (degree > oblate_body.geopotential_degree()) { + *coefficient = {0, 0}; + return m.Return(OK()); + } + *coefficient = {oblate_body.cos()[degree][order], + oblate_body.sin()[degree][order]}; + return m.Return(OK()); +} + +Status principia__ExternalGeopotentialGetReferenceRadius( + Plugin const* const plugin, + int const body_index, + double* const reference_radius) { + journal::Method m{ + {plugin, + body_index}, + {reference_radius}}; + if (plugin == nullptr) { + return m.Return( + MakeStatus(Error::INVALID_ARGUMENT, "|plugin| must not be null")); + } + if (!plugin->HasCelestial(body_index)) { + return m.Return(MakeStatus( + Error::NOT_FOUND, + absl::StrCat("No celestial with index ", body_index))); + } + auto const& body = *plugin->GetCelestial(body_index).body(); + if (!body.is_oblate()) { + *reference_radius = body.mean_radius() / Metre; + return m.Return(OK()); + } + auto const& oblate_body = dynamic_cast const&>(body); + *reference_radius = oblate_body.reference_radius() / Metre; + return m.Return(OK()); +} + Status principia__ExternalGetNearestPlannedCoastDegreesOfFreedom( Plugin const* const plugin, int const central_body_index, diff --git a/ksp_plugin_test/interface_external_test.cpp b/ksp_plugin_test/interface_external_test.cpp index ba0c8ccbe5..0e0e77376e 100644 --- a/ksp_plugin_test/interface_external_test.cpp +++ b/ksp_plugin_test/interface_external_test.cpp @@ -20,6 +20,7 @@ using ksp_plugin::PartId; using ksp_plugin::FakePlugin; using ksp_plugin::Vessel; using physics::SolarSystem; +using quantities::Sqrt; using quantities::si::Centi; using quantities::si::Hour; using quantities::si::Kilo; @@ -114,5 +115,74 @@ TEST_F(InterfaceExternalTest, GetNearestPlannedCoastDegreesOfFreedom) { Lt(1 * Centi(Metre) / Second))))); } +TEST_F(InterfaceExternalTest, Geopotential) { + XY coefficient; + double radius; + auto status = principia__ExternalGeopotentialGetCoefficient( + &plugin_, + SolarSystemFactory::Earth, + /*degree=*/2, + /*order=*/0, + &coefficient); + EXPECT_THAT(base::Status(static_cast(status.error), ""), IsOk()); + EXPECT_THAT(-coefficient.x * Sqrt(5), IsNear(1.08e-3)); + EXPECT_THAT(coefficient.y, Eq(0)); + + status = principia__ExternalGeopotentialGetCoefficient( + &plugin_, + SolarSystemFactory::Earth, + /*degree=*/3, + /*order=*/1, + &coefficient); + EXPECT_THAT(base::Status(static_cast(status.error), ""), IsOk()); + EXPECT_THAT(coefficient.x, IsNear(2.03e-6)); + EXPECT_THAT(coefficient.y, IsNear(0.248e-6)); + + status = principia__ExternalGeopotentialGetCoefficient( + &plugin_, + SolarSystemFactory::Earth, + /*degree=*/1729, + /*order=*/163, + &coefficient); + EXPECT_THAT(base::Status(static_cast(status.error), ""), IsOk()); + EXPECT_THAT(coefficient.x, Eq(0)); + EXPECT_THAT(coefficient.y, Eq(0)); + + status = principia__ExternalGeopotentialGetReferenceRadius( + &plugin_, + SolarSystemFactory::Earth, + &radius); + EXPECT_THAT(base::Status(static_cast(status.error), ""), IsOk()); + EXPECT_THAT(radius, Eq(6'378'136.3)); + + status = principia__ExternalGeopotentialGetCoefficient( + &plugin_, + SolarSystemFactory::Ariel, + /*degree=*/2, + /*order=*/2, + &coefficient); + EXPECT_THAT(base::Status(static_cast(status.error), ""), IsOk()); + EXPECT_THAT(coefficient.x, Eq(0)); + EXPECT_THAT(coefficient.y, Eq(0)); + + status = principia__ExternalGeopotentialGetCoefficient( + &plugin_, + SolarSystemFactory::Ariel, + /*degree=*/0, + /*order=*/0, + &coefficient); + EXPECT_THAT(base::Status(static_cast(status.error), ""), IsOk()); + EXPECT_THAT(coefficient.x, Eq(1)); + EXPECT_THAT(coefficient.y, Eq(0)); + + status = principia__ExternalGeopotentialGetReferenceRadius( + &plugin_, + SolarSystemFactory::Ariel, + &radius); + EXPECT_THAT(base::Status(static_cast(status.error), ""), IsOk()); + EXPECT_THAT(radius, Eq(578'900)); + +} + } // namespace interface } // namespace principia diff --git a/serialization/journal.proto b/serialization/journal.proto index b10f95182a..d8ed382158 100644 --- a/serialization/journal.proto +++ b/serialization/journal.proto @@ -151,7 +151,7 @@ message XY { } message Method { - extensions 5000 to 5999; // Last used: 5155. + extensions 5000 to 5999; // Last used: 5157. } message AdvanceTime { @@ -423,6 +423,55 @@ message ExternalFlowFreefall { optional Return return = 3; } +// Sets |coefficient| to the normalized geopotential coefficient of the given +// |degree| and |order| of the body with index |body_index|. +// |coefficient.x| is set to Cnm, |coefficient.y| is set to Snm. +message ExternalGeopotentialGetCoefficient { + extend Method { + optional ExternalGeopotentialGetCoefficient extension = 5156; + } + message In { + required fixed64 plugin = 1 [(pointer_to) = "Plugin const", + (is_subject) = true]; + required int32 body_index = 2; + required int32 degree = 3; + required int32 order = 4; + } + message Out { + required XY coefficient = 1; + } + message Return { + required Status result = 1; + } + optional In in = 1; + optional Out out = 2; + optional Return return = 3; +} + +// Sets |reference_radius| to the value in metres of the reference radius of +// the geopotential model for the body with index |body_index|. +// If the body does not have a geopotential model, i.e., if it is modeled as an +// isotropic sphere, the mean radius is used. +message ExternalGeopotentialGetReferenceRadius { + extend Method { + optional ExternalGeopotentialGetReferenceRadius extension = 5157; + } + message In { + required fixed64 plugin = 1 [(pointer_to) = "Plugin const", + (is_subject) = true]; + required int32 body_index = 2; + } + message Out { + required double reference_radius = 1; + } + message Return { + required Status result = 1; + } + optional In in = 1; + optional Out out = 2; + optional Return return = 3; +} + // Returns the first point of the coast phase following the given manoeuvre // which is locally nearest to the reference position, or the nearest endpoint // if there is no local minimum.