Keplerian orbits #279

dylan-mcleod · 2024-05-10T22:57:35Z

Added a Keplerian orbit type in scientific/kepler.h, to later build patched conics on top of

src/osp/scientific/kepler.cpp

jonesmz · 2024-05-11T05:33:26Z

src/osp/scientific/kepler.cpp

+    /**
+     * @brief Wrap an angle to the range [0, 2π), useful for keplerian orbits
+     */
+    inline constexpr double wrap_angle(double angle)


when practical use const for function parameters or local variables that aren't modified.

src/osp/scientific/kepler.h

src/osp/scientific/kepler.cpp

src/testapp/sessions/orbits.cpp

jonesmz · 2024-05-11T05:52:18Z

please also make sure the continuous integration passes before you finish your changes for this PR.

jonesmz · 2024-05-11T05:57:41Z

src/testapp/sessions/orbits.cpp

+        Session &uniCore,
+        Session &uniScnFrame)
+    {
+        feenableexcept(FE_INVALID | FE_OVERFLOW);


this function appears to only exist on linux?

I'm not sure enabling floating point exception trapping in a platform specific manner is the right approach.

if you're doing math operations that are prone to encountering floating point exceptions, i'd rather see standardized / standard-library functionality used only, unless there's a substantial performance improvement to be had from calling out to platform specific functions.

seems like this header has the standardized stuff in it? https://en.cppreference.com/w/cpp/numeric/fenv

src/osp/scientific/kepler.cpp

src/osp/scientific/kepler.h

Capital-Asterisk · 2024-05-11T06:56:26Z

src/osp/scientific/kepler.cpp

+    {
+        if (is_elliptic())
+        {
+#if 1 // Use standard kepler solver


whatever is going to be here in the final version, if constexpr can likely do the job here

src/osp/scientific/kepler.cpp

jonesmz · 2024-05-11T07:07:43Z

src/osp/scientific/kepler.cpp

+    return E;
+}
+
+


can you remove the extra newlines?

src/osp/scientific/kepler.cpp

…tform-specific M_PI constant

src/osp/scientific/kepler.cpp

jonesmz · 2024-05-11T16:41:19Z

src/osp/scientific/kepler.cpp

@@ -133,63 +143,63 @@ static double mean_motion(double const mu, double const a)

 static double get_orbiting_radius(double const ecc, double const semiMajorAxis, double const trueAnomaly)


all of these math functions look like they are noexcept right?

i am a noob at c++ but in (76-77) of kerpler.cpp double df_E = 1.0 - eccentricity * cos(E); // if eccentricity happens to be 1, a div by 0 , since c++ do-while loops execute at least once, idk if that causes issues in c++ / this specific implementation (ie it may check if eccentricity is 1 and use some other code) but every language i actually have skills in throws a exception

my understanding is that floating point exceptions are not the same as C++ exceptions. I think if it does a divide by zero, the program will simply crash

regardless, if crashing isn't the right answer, then what should the function do?

More accurately, on posix platforms it'll receive a signal that interrupts program flow and calls a signal handler (which we have no handler registered, so it'll default to crashing)

on windows it'll likely result in a structured exception handler exception being thrown (which is not the same thing as a C++ exception, and is handled completely differently, just like signal handling is done totally differently)

Regardless, throwing a c++ exception isnt the right answer, we need to have a solid determination of what we want the program to do if div by zero happens in this function.

Crashing is one choice, or we can return some kind of error code, but we want to document what the behaivor is going to be and then stick to it.

A floating point math error isn't necessarily an unrecoverable error for the end application, so I don't know if crashing is a good default choice choice for handling them.

I'm not well versed in c++ error handling though, so I don't know what the right way to handle floating point errors would be.

My recommendation is to add noexcept unless someone provides a signal-handler -> exception translator, which very well might happen.

Until then, noexcept can result in better code-gen in many situations, helps analysis tools, and it can always be removed later.

…ectories

jonesmz · 2024-06-23T12:33:14Z

src/osp/scientific/kepler.cpp

+ */
+static double solve_kepler_hyperbolic(double const meanAnomaly, double const eccentricity, int steps)
+{
+    double const tolerance = 1.0E-10;


you can make tolerance into a constexpr variable here can't you?

static constexpr tolerance = 1.0E-10;

jonesmz · 2024-06-23T12:34:52Z

src/osp/scientific/kepler.cpp

+    do
+    {
+        double f_E = eccentricity * sinh(E) - E - meanAnomaly;
+        double df_E = eccentricity * cosh(E) - 1.0;


f_E and df_E can both be const.

I know this might seem very silly to point out, but I've watched so many people make really brain dead mistakes because they've modified a variable that they intended to never be modified.

I think it's a serious mis-design of C++ that variables aren't immutable by default, and require a keyword to make mutable. Instead of the situation we have now.

jonesmz · 2024-06-23T12:35:49Z

src/osp/scientific/kepler.cpp

+        delta = f_E / df_E;
+        E = E - delta;
+        --steps;
+    } while (std::abs(delta) > tolerance && steps > 1);


minor: I recommend wrapping the individual portions of your conditional expression in parens.

Watched a lot of people mistake the order of evaluation here too... Had a few bugs in prod because of it at work.

jonesmz · 2024-06-23T12:36:11Z

src/osp/scientific/kepler.cpp

+    double const beta = ecc / (1 + std::sqrt(1 - ecc * ecc));
+    return E + 2.0 * atan2(beta * sin(E), 1 - beta * cos(E));
+
+    // return 2.0 * atan2(std::sqrt(1 + ecc) * sin(E / 2.0), std::sqrt(1 - ecc) * cos(E / 2.0));


minor: Best to remove unused/commented-out code

jonesmz · 2024-06-23T12:39:13Z

src/osp/scientific/kepler.cpp

+/**
+ * @brief Convert true anomaly (ν) to hyperbolic eccentric anomaly (E)
+ */ 
+static double true_anomaly_to_hyperbolic_eccentric_anomaly(double const &ecc, double const &v)


Is there any value in taking ecc and v as references instead of by value?

At least on x86_64 linux, with clang, sizeof(double) is 8bytes, which is the same as a pointer.

So on most platforms, when passing by value, the double should be passed via register (assuming few enough parameters). On platforms where there isnt enough space to pass by register, it'll implicitly pass by reference / pointer regardless.

jonesmz · 2024-06-23T12:42:09Z

src/osp/scientific/kepler.cpp

+
+    u = Vector3d(cos(a) * cos(b) - sin(a) * sin(b) * cos(c),
+                 cos(a) * sin(b) + sin(a) * cos(b) * cos(c),
+                 sin(a) * sin(c));


is it fair to assume that C++ standard library doesn't provide any special math functions that do much / any of this for us? i'm guessing that cos and sin are all we get?

jonesmz · 2024-06-23T12:42:58Z

src/osp/scientific/kepler.cpp

+        delta = f_E / df_E;
+        E = E - delta;
+
+        E = wrap_angle(E); // Sometimes this sometimes fails to converge without this line


sometimes this sometimes -> the algorithm sometimes

jonesmz · 2024-06-23T12:43:31Z

src/osp/scientific/kepler.cpp

+    // newton-rhapson
+    do
+    {
+        double f_E = eccentricity * sinh(E) - E - meanAnomaly;


minor: not a show-stopper, but i personally prefer to see functions from the standard library explicitly called out as coming from the std:: namespace.

jonesmz · 2024-06-23T12:45:54Z

src/osp/scientific/kepler.cpp

+                              Vector3d const v, 
+                              Vector3d const w,
+                              Vector3d      &position, 
+                              Vector3d      &velocity)


it doesn't appear that position and velocity are used as input parameters, but are just output parameters.

Consider returning std::pair<Vector3d, Vector3d> instead, and then unpacking the tuple at the place where the function is called, either using std::tie or just structured bindings, depending on how the calling code wants to.

This is less about codegen and more about just how using the function looks, but there could be advantages for named return value optimization.

jonesmz · 2024-06-23T12:47:44Z

src/osp/scientific/kepler.cpp

+    }
+}
+
+void KeplerOrbit::get_uvw_vectors(Vector3d &u, Vector3d &v, Vector3d &w) const


same as my comment for get_state_vectors

jonesmz · 2024-06-23T12:49:43Z

src/osp/scientific/kepler.cpp

+using namespace osp;
+
+// M_PI isn't standard, so define it here.
+static constexpr double PI = Magnum::Math::Constants<double>::pi();


you can use pi_v from https://en.cppreference.com/w/cpp/numeric/constants, if you prefer.

jonesmz · 2024-06-23T12:51:45Z

src/osp/scientific/kepler.cpp

+}
+
+std::optional<double> KeplerOrbit::get_next_time_radius_equals(const double r, const double currentTime) const {
+    std::optional<double> apoapsis = get_apoapsis();


i like std::optional a lot from a programmer ease of use perspective, but just wanted to double check that you don't consider it a barrier to good codegen compared to using explicit checks for NaN ?

jonesmz · 2024-06-23T12:52:31Z

src/osp/scientific/kepler.h

+
+    KeplerOrbitParams m_params;
+
+    KeplerOrbit(KeplerOrbitParams const params) : m_params(params) {}


KeplerOrbitParams is large enough that it should probably always be passed by const& to the function.

jonesmz · 2024-06-23T12:55:21Z

src/testapp/sessions/orbits.cpp

+    auto const [vx, vy, vz] = sat_views(rMainSpaceCommon.m_satVelocities, rMainSpaceCommon.m_data, planetCount);
+    auto const [qx, qy, qz, qw] = sat_views(rMainSpaceCommon.m_satRotations, rMainSpaceCommon.m_data, planetCount);
+
+    std::mt19937 gen(seed);


you might get slightly better random number generation performance with the _64 version of the mt19937 because it works off of a bigger block size of random data

dylan-mcleod added 2 commits May 10, 2024 14:48

Add keplerian orbits, as well as an orbit testing scenario

ba1a82b

Move some functions around in kepler.cpp

23d6a79

dylan-mcleod marked this pull request as draft May 11, 2024 00:20

jonesmz reviewed May 11, 2024

View reviewed changes

src/osp/scientific/kepler.h Outdated Show resolved Hide resolved

jonesmz reviewed May 11, 2024

View reviewed changes

src/osp/scientific/kepler.h Outdated Show resolved Hide resolved

jonesmz reviewed May 11, 2024

View reviewed changes

src/osp/scientific/kepler.cpp Outdated Show resolved Hide resolved

jonesmz reviewed May 11, 2024

View reviewed changes

src/osp/scientific/kepler.cpp Outdated Show resolved Hide resolved

jonesmz reviewed May 11, 2024

View reviewed changes

src/osp/scientific/kepler.cpp Outdated Show resolved Hide resolved

jonesmz reviewed May 11, 2024

View reviewed changes

src/osp/scientific/kepler.cpp Outdated Show resolved Hide resolved

jonesmz reviewed May 11, 2024

View reviewed changes

src/testapp/sessions/orbits.cpp Outdated Show resolved Hide resolved

jonesmz reviewed May 11, 2024

View reviewed changes

src/testapp/sessions/orbits.cpp Show resolved Hide resolved

jonesmz approved these changes May 11, 2024

View reviewed changes

jonesmz reviewed May 11, 2024

View reviewed changes

Capital-Asterisk reviewed May 11, 2024

View reviewed changes

src/osp/scientific/kepler.cpp Show resolved Hide resolved

Capital-Asterisk reviewed May 11, 2024

View reviewed changes

src/osp/scientific/kepler.cpp Outdated Show resolved Hide resolved

Incorporate reviewer feedback

899bea8

Capital-Asterisk reviewed May 11, 2024

View reviewed changes

src/osp/scientific/kepler.h Outdated Show resolved Hide resolved

Capital-Asterisk reviewed May 11, 2024

View reviewed changes

src/osp/scientific/kepler.cpp Outdated Show resolved Hide resolved

jonesmz reviewed May 11, 2024

View reviewed changes

src/osp/scientific/kepler.cpp Outdated

return E;

}

Copy link

Member

jonesmz May 11, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can you remove the extra newlines?

jonesmz reviewed May 11, 2024

View reviewed changes

src/osp/scientific/kepler.cpp Show resolved Hide resolved

dylan-mcleod added 4 commits May 11, 2024 00:12

Remove std prefix from trigonometric functions, and remove use of pla…

1720587

…tform-specific M_PI constant

Change orbital elements to live inside internal PoD struct

43f468b

Resolve warnings and MSVC build errors

b32de43

Add math for determining an orbit's time until escape

77af67f

jonesmz reviewed May 11, 2024

View reviewed changes

src/osp/scientific/kepler.cpp Show resolved Hide resolved

jonesmz reviewed May 11, 2024

View reviewed changes

Add doxygen documentation, make orbit params struct non-private

53146a0

dylan-mcleod added 2 commits May 13, 2024 16:28

Move a few things around in kepler.cpp/.h and fix some edge case traj…

1df4986

…ectories

Improve legibility of kepler.cpp

5170598

jonesmz reviewed Jun 23, 2024

View reviewed changes

jonesmz requested a review from EhWhoAmI June 23, 2024 12:41

jonesmz reviewed Jun 23, 2024

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Keplerian orbits #279

Keplerian orbits #279

dylan-mcleod commented May 10, 2024

jonesmz May 11, 2024

jonesmz commented May 11, 2024

jonesmz May 11, 2024

Capital-Asterisk May 11, 2024

jonesmz May 11, 2024

jonesmz May 11, 2024

tigercoding56 May 11, 2024 •

edited

Loading

jonesmz May 11, 2024

jonesmz May 11, 2024 •

edited

Loading

dylan-mcleod May 11, 2024

jonesmz Jun 23, 2024 •

edited

Loading

jonesmz Jun 23, 2024

jonesmz Jun 23, 2024

jonesmz Jun 23, 2024

jonesmz Jun 23, 2024

jonesmz Jun 23, 2024

jonesmz Jun 23, 2024

jonesmz Jun 23, 2024

jonesmz Jun 23, 2024

jonesmz Jun 23, 2024 •

edited

Loading

jonesmz Jun 23, 2024

jonesmz Jun 23, 2024

jonesmz Jun 23, 2024

jonesmz Jun 23, 2024

jonesmz Jun 23, 2024 •

edited

Loading

		@@ -133,63 +143,63 @@ static double mean_motion(double const mu, double const a)

		static double get_orbiting_radius(double const ecc, double const semiMajorAxis, double const trueAnomaly)


		KeplerOrbitParams m_params;

		KeplerOrbit(KeplerOrbitParams const params) : m_params(params) {}

Keplerian orbits #279

Are you sure you want to change the base?

Keplerian orbits #279

Conversation

dylan-mcleod commented May 10, 2024

Choose a reason for hiding this comment

jonesmz commented May 11, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tigercoding56 May 11, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jonesmz May 11, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jonesmz Jun 23, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jonesmz Jun 23, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jonesmz Jun 23, 2024 • edited Loading

Choose a reason for hiding this comment

tigercoding56 May 11, 2024 •

edited

Loading

jonesmz May 11, 2024 •

edited

Loading

jonesmz Jun 23, 2024 •

edited

Loading

jonesmz Jun 23, 2024 •

edited

Loading

jonesmz Jun 23, 2024 •

edited

Loading