Find Orbital Elements given Orbital State Vectors

[FireStorm000]

Hi everyone,
I've been banging my head against this problem for a few days now so I'm ready to find some help with it.
In short, I have a set of orbital state vectors(a position vector and velocity vector) describing an object's starting position, as well as the mass of the focus, and I need to figure out the math of how to convert to conventional orbital elements.
The biggest hitch seems to be determining the length of the orbit's semi-major axis. I've found one formula which describes near-circular orbits very well, but it seems to fail on highly eccentric orbits:

a = (2 / r - v^2 / ( m * G))^-1

In the case of a circular orbit, the equation is perfect. For example, plugging in Earth's orbit about the sun produces an answer very close to 1 AU, which is satisfactory to me.
Now, however, consider the case of a highly eccentric, or even radial orbit. Let's say that the orbit has zero initial velocity, and starts at a distance of 1 AU. The orbiting body will fall in towards the sun and (considering an ideal case with point masses) continue through the other side, and come to rest 1 AU on the other side. That yields an ellipse with a semi-major axis equal to r, or 1 AU, however simplifying this equation will give you one half r as the semi-major axis.
So I guess that comes down to if anyone knows or can find a better equation for me to use; one that will hopefully work in the case of elliptic, hyperbolic, parabolic, and radial trajectories.
Thanks in advance for any help
--Firestorm--

 

[Filip Larsen]

Welcome to PF!

If the initial velocity is zero at some finite distance r from the primary body then you have a case of a so-called elliptic rectilinear orbit which is a degenerated ellipse where the orbital eccentricity is 1 (or very close to 1) while perihelion (closest distance to the Sun in orbit) is 0 (or very close to 0).

Getting a semi-major axis of a = 1/2 r is correct since the orbiting object will be "wipped around" the primary body and head out the same way it came in. If you draw the orbit it will be a line (flat ellipse) from the distance r down to the primary, with the center of the ellipse being half-way between r and the primary body. The reason this is so is that the orbital model assumes that the masses are point-like, that is, the model does not include any effects due to impact or due to the orbiting body "entering" or "passing through" the primary. In cases where an orbital body comes so close to its primary that non-gravitational forces becomes significant you need to model this phase of the orbit/impact in a different way if you want to know what physically would happen.

That aside, I'm sure you are aware that no natural body will be able to pass through the sun, so if you are not interested in the impact but merely made a theoretical "simplification" to understand purely orbital dynamics, you should model the sun as a point-like mass that you cannot collide with.

 

[FireStorm000]

Thanks! Looks like that formula works after all.

I'm using this to write a simulation engine of orbital dynamics; major bodies such as the planets, and even things like binary star systems, can be estimated to sufficient accuracy for my purposes with less processor time using Kepler orbits so that's what I'm using where possible. I think I can figure out everything else knowing the semi-major axis.

 

[BobG]

a = (2 / r - v^2 / ( m * G))^-1

This works because of relationship between specific energy and the semi-major axis:

ε = GM/(-2a) or a = GM/(-2ε)

and specific energy is also equal to the sum of kinetic energy and potential energy:

ε = v^2/2 - GM/r

In other words, get your specific energy and your angular momentum vector and everything else flows from those two parameters.

Technically, your equation just bypasses calculating the specific energy, but it's always nice to know why you're using that particular equation. That same principle of conservation of energy/conservation of angular momentum comes into play in determining your eccentricity and it also often comes into play in determining your state vector in the first place, since few antennas can calculate satellite velocity.

 

[FireStorm000]

So, I was at last able to derive that equation from conservation of energy and conservation of angular momentum; I'll attach the derivation for anyone who's interested.

I do have another issue that I'd like some advice on: Now that I have everything nicely in terms of orbital elements, is there an easy way to derive a position with respect to time function? For the sake of simplicity, I'm working in 2D. This is for a computer program, so ideally I want it to be executed very fast, as there will be many objects to deal with. I've already found Kepler's equation, though an s(t) function is anything but simple with that. I want to stay away from using Newtons force law as that has problems with taking loads of processor power and accuracy over multiple iterations(orbits calculated with the force law are unstable). Does anyone know how to get position from orbital elements without evaluating Kepler's equation again for every different time?

 

[BobG]

Does anyone know how to get position from orbital elements without evaluating Kepler's equation again for every different time?

Generally, solving Kepler's equation is the trade-off for using mean anomaly instead of true anomaly. For tracking a satellite using orbital elements that may be days old, it's a good trade-off.

If you're trying to do orbit simulations with short time steps, you pretty much have to accept that it's going to take a lot of calculations with a fast computer. The most aesthetic way to do it is to calculate your position & velocity for each time step in your time span and put them in an ephemeris file. Your simulation can use the ephemeris file to create the orbit, display its position, etc. At least the delay is all on the front end before your simulation starts instead of slowing down the simulation itself.

The drawback is that your ephemeris file can be huge if you're simply plotting the positions. So your ephemeris file will probably have larger time steps than your simulation and your simulation will interpolate between each time step (which is why you need the velocity as well as the position).

(Nothing original in this idea - that's how your high end commercial simulators, such as STK, do it when using SGP4 elsets.)

 

[gneill]

Depending upon the accuracy required you may be able to get away with a simple integration method for your orbit simulation. Look up terms like "leapfrog integrator" and "verlet integration". These are not particularly CPU-intensive.

 

[gneill]

I can demonstrate a fairly quick way to get at the parameters of your orbit given a velocity and position vector. I won't derive the formulae used here, but they can generally be found in any good book on Astrophysics.

Given: \vec{r}, \vec{v} are the position and velocity vectors respectively, then the specific angular momentum (pseudo)vector is given by:

\vec{h} = \vec{r} \times \vec{v}

The eccentricity vector, \vec{e} has a magnitude equal to the eccentricity constant for the orbit, and has a direction that points toward periapsis (this can be convenient to know!). This vector is given by:

\vec{e} = \frac{\vec{v} \times \vec{h}}{\mu} - \frac{\vec{r}}{r}

where \mu is the gravitational parameter of the primary (G*M).
The semi latus-rectum is given by:

p = \frac{h^2}{\mu}

Now you can find periapsis distance, apoapsis distance, semimajor axis,...:

r_p = \frac{p}{1 + e}

r_a = \frac{p}{1 - e}

a = \frac{p}{1 - e^2}