Three-Body orbital simulation reference plane

[AussieDave]

Hello,

I have a little three-body (well n-body but let's just worry about 3 now) simulator. I input six Keplerian orbital elements (period P, eccentricity E, inclination I, longitude of the ascending node W, argument of periapse w and initial true anomaly f). They then get converted to Cartesian coordinates and integrated using a simple RK4.

A problem, however, is the plane with which the angles are measured with respect to. I also have code that converts from Cartesian back to orbital elements, so I can see how they progress over time. This seems to be with respect to the invariable plane though, which is the plane of the total angular momentum of the system (and due to conservation of angular momentum it will not change).

When comparing my results to papers, I see I am using the invariable plane as the reference. I am able to successfully go from orbital elements -> Cartesian -> orbital elements and recover them, which makes me believe I am both inputting and outputting the orbital elements with respect to the invariable plane.

What I want though is to use the plane of the sky, which is how orbital elements are defined in exoplanet papers since we measure the orbits...in the sky. I think it would be just a rotation between the two but I am not sure what it is exactly and my tests so far have failed. I have also looked a fair bit online but not seen what the transformation is. Quite possibly it's simple and I just made an error, but I am not sure.

Kind regards in advance,

Dave

 

[BobG]

Which plane? The equatorial plane or the ecliptic plane? In other words, is it a simulator for planetary orbits in the Solar System or satellite orbits around the Earth? (It can't be generic, since you have to have the mass of the central body established somehow - either by default or manually)

I have no idea what you mean by "plane of the sky". Do you mean in topocentric-horizon coordinates?

Regardless, your coordinate transformations require both rotation and translation. The translation accounts for the difference in origins of the two coordinate systems. You get a break in going from your Keplerian orbital elements to geocentric Cartesian (or heliocentric Cartesian for planets) in that the origins are the same.

But when the origins of the two coordinate systems aren't the same, you have to account for the difference between the two by addition or subtraction. And, obviously, the two vectors have to be rotated into the same coordinate system (or parallel coordinate systems, since you haven't done the translation yet) in order to add or subtract them.

 

[lpetrich]

A good reference plane would be the one perpendicular to the total angular-momentum vector, the "invariable plane". One can then select a direction in that plane as a reference in it.

Here's an alternative approach: calculate the angular-momentum and the Laplace-Runge-Lenz (eccentricity) vectors for each orbit as you go. You can then plot them in rectangular-coordinate fashion.

 

[AussieDave]

I am looking at the orbits of exoplanets.

I understand Ipetrich the invariable plane, as that is what the orbital elements are currently defined with respect to in my program, but I want to be able to input the orbital elements as seen in discovery papers.

For example, by taking the X-Y plane as that of the sky, a transiting exoplanet will have have an inclination of ~90, but with respect to the invariable plane it will have an inclination of ~0

I want to know how to convert between the two. I have attached a screenshot from the Kepler 16 transiting circumbinary paper, and its supplementary online material. In the discovery paper (on the right) you see the sky-plane coordinates, with an inclination of ~90. On the left you see the variation in orbital elements over time with respect to the invariable plane, where the inclination is now ~0. I just don't know what the conversion is :(

 

[lpetrich]

AussieDave, I think that you ought to look for stuff on binary-star coordinate systems. Details should be easy to find, and exoplanet orbits use the same conventions as binary-star orbits. The inclinations are relative to the plane of the sky (our sky, of course), but I'd have to look up further details.

Direct detection (Fomalhaut B) and astrometric detection (Barnard's Star's now-discredited planets) are like visual binaries. They can get all six elements, though the inclination has a sign ambiguity, and the semimajor axis is angular and may be multiplied by (planet mass) / (star mass).

Radial-velocity and transit detections are like spectroscopic and eclipsing binaries. Radial velocity cannot get the orbit inclination, though transit detection can. Nether method can get how the orbit is rotated along the line of sight to the system.


Also, could you please shrink your inlined image a bit? Half its current width would be OK. It's currently at 2298*930 pixels. I ask that because your image messes up the formatting and makes this page hard to read.