Coding Orbital Mechanics

I forget what inspired me, probably discussions about propellant depots, but after an hour of two playing around in the great spaceflight simulator Orbiter I decided to give the old orbital mechanics another go.

I've written dynamic physics simulators before which have really nice accuracy and can handle the big numbers required for orbital simulations, including three body problems, but I never managed to tackle the classical orbital mechanics - too much Greek terminology I guess - until today.

I set myself a somewhat difficult task:

Given a position vector r and a velocity vector v, somewhere above a reference (say, a planetary body) calculate all the required classical orbital elements:

a - semi-major axis
e - eccentricity
i - inclination
l - longitude of the ascending node
w - argument of periapsis
t - true anomaly

There's others, but they can be derived from these. You'll also note I haven't used any Greek letters that my keyboard doesn't have. I can also do the reverse: convert classical orbital elements back into state vector elements.

So what is this good for? The true anomaly tells us where the object is in the orbit relative to the lowest point in the orbit (the periapsis) and has an interesting operation: you can add (or subtract) time to the orbit to get a new position. To do this, one calculates the mean anomaly and adds to it a multiple of the mean motion. From the mean anomaly you can calculate the eccentric anomaly, but this has to be done numerically. Then the eccentric anomaly can be converted into a true anomaly and we're done - but don't forget to correct for range and remember what quadrant of the orbit you're in.

    void step(double time)
    {
        double M = meanAnomaly();
        M += meanMotion() * time;
        while (M < -2*M_PI)
            M = M + 2*M_PI;
        if (M < 0)
            M = 2*M_PI + M;
        while (M > 2*M_PI)
            M = M - 2*M_PI;
        double E = calcEccentricAnomaly(M);
        calcTrueAnomaly(E);
    }

This means we can add any duration of time, with any desired accuracy, and expect to get an accurate result. In comparison, numerical simulations of newton's laws often give more inaccurate results the longer you run them. On the other hand, doing accurate orbital perturbation and rocket burn simulation with classical orbital elements is very difficult, and accurate three body solutions are still considered unsolved problems. As such, having the ability to accurately convert from these two very different systems gives you the best of both worlds.

Here is the code, in all it's C++ glory. Feel free to use it for whatever you like. In writing it I found the descriptions on orbital state vectors from Wikipedia particularly helpful, along with appendix C of the Orbiter manual and this Numerit application note.

How about that! I combined the two topics of this blog in a single post.