Converting State Vectors to Keplerian Orbital Elements for Binary Objects

[SethFusion]

Homework Statement:: I'm working on a personal project to convert objects from a simulation using state vectors for position and velocity to Keplerian orbital elements (semimajor axis, eccentricity, argument of periapsis, etc.). However, the equations I am using do not calculate the eccentricity correctly for binary objects.
Relevant Equations:: mu = gravitation parameter * mass of parent body
momentum = position x velocity
eccentricity = ((velocity x momentum) / mu) - (position / positionMagnitude)
where velocity, mommentum, and position are vectors.

Specifically I have been using this short paper:
https://downloads.rene-schwarz.com/...State_Vectors_to_Keplerian_Orbit_Elements.pdf

Now, I've been working on this on and off for a while, so thank you for taking the time to read through this. Hopefully someone can give me some insight as to what the problem is, because I'm a computer science major, not a physics major.

The program I am working on converts objects from a program called Universe Sandbox (US) into objects for another program called Space Engine (SE). More specifically, this is for solar systems, real or imagined, and deals with stars, planets, moons, and any celestial objects. US is a simulator and stores information about objects using cartesian state vectors, which can be exported to represent a snapshot of the simulation. SE is an emulator and only uses simple Keplerian elements.

The paper I linked above has worked near-perfect so far and I have had few problems. But eccentricity is not calculated correctly for any binary relationships, and I'm not sure why. Here is the relevant code:

void CalcOrbit(Object& obj)
{
    for (int i = 0; i < obj.child.size(); i++)
        CalcOrbit(*obj.child.at(i));

    if (&obj == root)
        return;

    double mu = G * obj.parent->mass;

    obj.position = Subtract(obj.position, obj.parent->position);
    obj.velocity = Subtract(obj.velocity, obj.parent->velocity);

    StateVect momentVect;
    // calculate the momentum vector h
    momentVect = CrossProduct(obj.position, obj.velocity);

    .
    .
    .

    // calculate eccentricity vector and eccentricity
    StateVect eccentVect, VcrossH = CrossProduct(obj.velocity, momentVect);
    eccentVect.x = ((VcrossH.x / mu) - (obj.position.x / obj.position.magnitude()));
    eccentVect.y = ((VcrossH.y / mu) - (obj.position.y / obj.position.magnitude()));
    eccentVect.z = ((VcrossH.z / mu) - (obj.position.z / obj.position.magnitude()));
    obj.eccentricity = eccentVect.magnitude();
    
    .
    .
    .
}

Let me give a little context for this code: What happens here is that every object in a solar system is connected through vectors of "children" and every object has a parent. For instance, a star would be at the top of the hierarchy and a moon would be at the bottom. Objects are sent through this CalcOrbit function from the bottom of the hierarchy (this is why the parent's position/velocity are subtracted from the object's), and their Keplerian elements are calculated as necessary.

Just to be clear, this equation works fine for most objects. But when an object passes through that is in a binary relationship, its "parent" is the barycenter between it and its partner. The barycenter has a total mass of the combined mass of the two objects. When eccentricity is calculated, it often comes out wildly incorrect (often around 0.9).

So I have tried to use other formulas, but had no success (I may have implemented them wrong, though). I have also tried changing the value of mu to its partner's mass, the barycenter's parent's mass, and other values, but none have worked. Unfortunately, I think I don't have the background in physics needed to know what the actual problem is, which is why I'm posting here. Is this problem even possible? Or are binary relationships too complicated to derive an eccentricity from a single snapshot of position/velocity?

If anyone has any ideas about why this doesn't work, or if I haven't provided enough context, please let me know.

The full code can be found here: https://github.com/SethFusion/US2-to-SE-Converter
But be warned that most of the program is dedicated to identifying objects and their relationships, and much of it I didn't write myself. The relevant part will be the full CalcOrbit function around line 353.

Also, I apologize if this is not the correct forum to post this in. I was afraid there was too much code to post in the regular physics forum. Mods, please feel free to move this to the correct forum or to contact me about where would be the correct place for this post.

 

[pbuk]

Thank you for posting this, it has solved a difficult Christmas present problem!

I don't have time to look at this fully at the moment, but briefly:

• This is not the best topic for this here, I'll get it moved.
• If mu in your code is the barycentric mass (i.e. the sum of the masses of the two elements) then it is not the same as $\mu$ in the paper you linked which is given by [Edit: $\mu = GM_b$ - originally posted incorrectly as $\mu = \frac{M_b}{G}$]. You need to divide the sum of the masses by the gravitational constant $G \approx 6.67428 \times 10^{-11} \mathrm {\frac{m^3}{kg \cdot s^2}}$, being careful to get the units right.

• I have also tried changing the value of mu to its partner's mass, the barycenter's parent's mass, and other values, but none have worked.

  Orbital mechanics is much more complicated than it looks; stabbing around in the dark like this is not a good idea. There are some good primers around suitable for a computational background, I'll dig out some links.

 

[SethFusion]

Orbital mechanics is much more complicated than it looks; stabbing around in the dark like this is not a good idea. There are some good primers around suitable for a computational background, I'll dig out some links.

In my heart, I knew none of these would work. But I also wasn't sure where to start looking for a solution, so I was a bit desperate. Thank you for your advice and, if you could find those links that might help me, I would greatly appreciate it.

 

[SethFusion]

• If mu in your code is the barycentric mass (i.e. the sum of the masses of the two elements) then it is not the same as $\mu$ in the paper you linked which is given by $\mu = \frac{M_b}{G}$. You need to divide the sum of the masses by the gravitational constant $G \approx 6.67428 \times 10^{-11} \mathrm {\frac{m^3}{kg \cdot s^2}}$, being careful to get the units right.

I want to make sure I'm understanding you right, so please bear with me if I ask a dumb question. It seems to me that in the paper, mu is defined as mass multiplied by G, not divided by G. This would cancel out the kg in the denominator of the units, and I believe the eccentricity equation itself uses units m³ / s²

Does this make sense? That is currently how I calculate mu for all objects since every object's parent body could be different. And like I said, this works fine for simple star-planet or planet-moon relationships but doesn't work for binary relationships.

If I am overlooking something important, please let me know.

 

[Mark44]

Line 6 of your code doesn't look right to me. obj is a reference parameter in CalcOrbit(), but you have this in line 6:

    if (&obj == root)

It could be that the compiler does the right thing anyway, but the comparison above should be obj == root, I believe, not &obj == root.

 

[SethFusion]

It could be that the compiler does the right thing anyway, but the comparison above should be obj == root, I believe, not &obj == root.

Right, this is just a check to see if the current object is at the top of the hierarchy, which means it wouldn't have an orbit. Root is a pointer to the top-most object. Lines 3-4 recursively call the function so that orbits are calculated from the bottom-up.

I probably could have left that code out, since it isn't directly related to the issue at hand. Sorry for the confusion.

 

[pbuk]

I want to make sure I'm understanding you right, so please bear with me if I ask a dumb question. It seems to me that in the paper, mu is defined as mass multiplied by G, not divided by G. This would cancel out the kg in the denominator of the units, and I believe the eccentricity equation itself uses units m³ / s²

Does this make sense? That is currently how I calculate mu for all objects since every object's parent body could be different. And like I said, this works fine for simple star-planet or planet-moon relationships but doesn't work for binary relationships.

If I am overlooking something important, please let me know.

No, of course you are right: $\mu = GM$, and now I am not rushing I can see that you are calculating this correctly on line 9.

There used to be a good primer on the JPL/NASA site but I can't find it now. For a book that contains exactly as much maths as you need and no more, Howard Curtis's Orbital Mechanics for Engineering Students is great. There are copies of this posted on the internet in breach of copyright, please do not post links here.

You also need to understand the concept of roundoff error in order to know for instance whether the results of these calculations

    obj.position = Subtract(obj.position, obj.parent->position);
    obj.velocity = Subtract(obj.velocity, obj.parent->velocity);

are meaningful (my apologies if this is all basic to you, but I am trying not to make any assumptions).

Now the next important point is that all orbits are calculated by considering the combined mass and the barycentre. It's just that when we have a planet orbiting its star, or a moon orbiting its planet, the relative size of the masses and the distance between them is such that the difference between the barycentre and the position of the more massive object is negligible.

So exactly the same calculations should apply for binary-barycentre orbits as for moon-planet orbits (bearing in mind that no orbit is exactly Keplerian where the total number of bodies in the system is greater than 2).

How do you know that the calculated eccentricities are wrong?

 

[SethFusion]

Now the next important point is that all orbits are calculated by considering the combined mass and the barycentre. It's just that when we have a planet orbiting its star, or a moon orbiting its planet, the relative size of the masses and the distance between them is such that the difference between the barycentre and the position of the more massive object is negligible.

So exactly the same calculations should apply for binary-barycentre orbits as for moon-planet orbits (bearing in mind that no orbit is exactly Keplerian where the total number of bodies in the system is greater than 2).

How do you know that the calculated eccentricities are wrong?

I certainly don't mind going over the basics. The subtraction of the parent's position and velocity is meaningful because the position of the parent object in the simulation is not always at the center of the simulation grid. For example, a star could have position (100,000, 0, 0) and the orbit of the planets are relative to that. The subtraction centers the current relationship on (0, 0, 0).

And yes, I would think the exact same calculations would work for all of the reasons you gave. I know the calculations are wrong because the program outputs the final values of the objects being converted. And Visual Studio allows you to watch how variables change during runtime step-by-step, so I know it is the formula that outputs these numbers. Here is an example output for the Earth-Moon system:

Barycenter                    "Earth-Moon"
{
    ParentBody            "Sun"
    Orbit
    {
        RefPlane        "Equator"
        SemiMajorAxis    0.999347
        Period            0.998924
        Eccentricity    0.0165413 // this value is roughly correct for the Eath's orbit around the Sun
        Inclination        0.00224701
        AscendingNode    12.6919
        ArgOfPericenter    271.956
        MeanAnomaly        -84.1143
    }
}

Planet                    "Earth"
{
    ParentBody            "Earth-Moon"
    Orbit
    {
        RefPlane        "Equator"
        SemiMajorAxis    1.54412e-05
        Period            0.071725
        Eccentricity    0.999998
        Inclination        4.91858
        AscendingNode    301.961   
        ArgOfPericenter    17.8443
        MeanAnomaly        179.996
    }
}

Moon                    "Moon"
{
    ParentBody            "Earth-Moon"
    Orbit
    {
        RefPlane        "Equator"
        SemiMajorAxis    0.00248531
        Period            0.071725
        Eccentricity    0.0250595 // this is closer, but still not accurate
        Inclination        4.91858
        AscendingNode    301.961
        ArgOfPericenter    224.996
        MeanAnomaly        151.516
    }
}

Now here is some context. There are three objects output here: The Earth-Moon barycenter, the Earth, and the Moon. You'll see that the barycenter was calculated correctly (I have had problems with ascending node and arg of pericenter, but that is for another topic), roundoff error is responsible for the difference between the real values.

Now when you look at the orbits of the Earth and Moon (measured relative to the barycenter), their eccentricities are not correct. The Moon is close, but I think this has something to do with distance from the barycenter. Other tests where the two objects are similar in mass have shown both values are wildly wrong.

This, I suppose, is why I'm posting here. I'm clearly missing something, but I'm not sure how to edit the formulas to account for whatever is missing.

 

[pbuk]

Barycenter                    "Earth-Moon"
{
    ParentBody            "Sun"
    Orbit
    {
        RefPlane        "Equator"
        SemiMajorAxis    0.999347
        Period            0.998924
        Eccentricity    0.0165413 // this value is roughly correct for the Eath's orbit around the Sun
        Inclination        0.00224701
        AscendingNode    12.6919
        ArgOfPericenter    271.956
        MeanAnomaly        -84.1143
    }
}

It looks like you are using AU for distance and some form of years (sideral?) for time, and heaven(!) knows what for mass - Earth masses? Are you sure you have got the correct value for $G$ in ${\mathrm{\dfrac{AU^3}{M_E \cdot {y_{sid}}^2}}}$? Personally I like to use metres, seconds and kilograms everywhere, but km and days (86400 seconds) are OK too.

 

[pbuk]

I don't like any of what you (or rather Roshex) are doing with binaries, I've added a comment on GitHub but I think the problem is more widespread than that.

For these Keplerian calculations to work you need to completely ignore the orbital partner - the only thing that matters is the body under consideration and the barycentre.

 

[SethFusion]

It looks like you are using AU for distance and some form of years (sideral?) for time, and heaven(!) knows what for mass - Earth masses? Are you sure you have got the correct value for $G$ in ${\mathrm{\dfrac{AU^3}{M_E \cdot {y_{sid}}^2}}}$? Personally I like to use metres, seconds and kilograms everywhere, but km and days (86400 seconds) are OK too.

I know the units look messy, but this is because the output file here is used as an input for the Space Engine program, which takes a very specific format to render objects. The real calculations are simpler than this, and then converted into a form Space Engine will accept. The original units are meters, kg, and seconds.

I don't like any of what you (or rather Roshex) are doing with binaries, I've added a comment on GitHub but I think the problem is more widespread than that.

For these Keplerian calculations to work you need to completely ignore the orbital partner - the only thing that matters is the body under consideration and the barycentre.

Quite frankly, I don't know exactly why he wrote that bit of code in that way. But eccentricity was broken long before that was written, and his code seemed to work for its purpose so I left it. But the only place the partner is considered is when calculating the semimajor axis. At the time an object's eccentricity is calculated, only the object and its parent are considered. In this case, an object and its barycenter.

I agree with you that we are doing something wrong with binaries. But I don't think that is the issue, unless I'm overlooking something simple.

 

[pbuk]

The original units are meters, kg, and seconds.

- and you are clearly handling the conversions correctly as most of the results are right.

I agree with you that we are doing something wrong with binaries.

Or you are doing something wrong with every body but it is only showing up with (what you are calling) binaries.

Have you worked through how the calculations should work, given some actual ephemerides?Try the following (in AU/kg/d units just for fun) from Horizons.

2020-Nov-18 00:00:00.0000
Earth GM, km^3/s^2             = 398600.435436
X = 5.489949218109829E-01 Y = 8.240267486214050E-01 Z = 5.688932676260653E-05
VX=-1.452755470046848E-02 VY= 9.595265033031247E-03 VZ= 4.803523735379163E-07

Luna GM, km^3/s^2          = 4902.800066
X = 5.491722934830898E-01 Y = 8.215687915965815E-01 Z = 8.241227889171440E-07
VX=-1.391405630742886E-02 VY= 9.603350712330360E-03 VZ=-5.486534523053657E-05

Earth-Moon Barycenter
 X = 5.489970769804319E-01 Y = 8.239968830074429E-01 Z = 5.620810177708234E-05
 VX=-1.452010033654438E-02 VY= 9.595363278758949E-03 VZ=-1.921301891792545E-07

 

[pbuk]

I've had a bit more of a look at this,

• I can't see anything wrong with your code, it seems to be implementing the calculations in the paper you linked.
• I can't see anything wrong with the calculations in that paper (although the referencing is terrible - Wikipedia!), for instance the eccentricity equation $\mathbf {e} = \dfrac { \mathbf {\dot r} \times \mathbf {h}}{\mu} - \dfrac{\mathbf {r}}{r}$ is Equation 2.30 (rearranged) in Curtis's book.
• I checked a few calculations using Horizons data of current state vectors. I get Jupiter's eccentricity accurate (in the sense that it matches Horizon's figure) to 2%, Earth's to about 16% but the Moon's is way out and apparently Titan is on its way out of the Solar System!
• I took a step back and realized that this is a pointless task - in any gravity field with more than 2 bodies, instantaneous state vectors cannot be translated into meaningful Keplerian orbits. To do this you need to fit the orbital parameters to the state vectors over at least one orbit, either by some averaging calculation or by optimising a best fit (which is I believe how JPL does it).

So I don't think there is anything wrong with your code, just what you are trying to achieve with it. Good luck!

 

[SethFusion]

• I can't see anything wrong with the calculations in that paper (although the referencing is terrible - Wikipedia!), for instance the eccentricity equation $\mathbf {e} = \dfrac { \mathbf {\dot r} \times \mathbf {h}}{\mu} - \dfrac{\mathbf {r}}{r}$ is Equation 2.30 (rearranged) in Curtis's book.

I've been taking a look through the Curtis book you mentioned. It's a good resource and much more useful than the paper I started out with. Thank you for that.

• I took a step back and realized that this is a pointless task - in any gravity field with more than 2 bodies, instantaneous state vectors cannot be translated into meaningful Keplerian orbits. To do this you need to fit the orbital parameters to the state vectors over at least one orbit, either by some averaging calculation or by optimising a best fit (which is I believe how JPL does it).

I was afraid this was the case, but I didn't have the background or mathematical skill to prove it. I'm glad you told me this. Even if it can't be done in this way, my mind is at ease now. I just want to thank you for taking time out of your day to look through my amateur code and for explaining the reasoning behind the problem.

Unfortunately, it isn't possible in this situation to follow objects along their orbits to obtain that data, so for this program, it is something that will just have to be fixed manually.