# Orbit Determination from Position and Velocity

I'm working on a game, somewhat similar to KSP (no building rockets, just beautiful orbital acrobatics). It's in a "small" star system, radius of 6 AU, and it features celestial bodies which orbit in deterministic keplerian orbits, and has a playerSpaceShip which can move freely between them.

It's very important to be able to perform orbital maneuvers in this game (prograde, retrograde, normal, antinormal, radial in/out, plane matching, rendezvous, etc) and to be able to see your playerShip orbit trajectory change as you maneuver. The playerShip can rotate and translate (using thrust) freely, and the orbit trajectory should change with the player's movements. KSP does this, for instance.

Time acceleration is also a feature, for obvious reasons considering the scale.

I'm working in Godot 4, and I made a copy of the Vector3 object that uses 64 bit floats to account for the size of my space ( i can't compile godot in 64 bit to get those sweet, sweet 64-bit Vector3 objects because i don't think the target platform will allow that build ). Tracking where all of these objects in the space are happens in a dedicated singleton separate from rendering, and data is pulled to be rendered as needed. This should allow for centimeter-level precision (for docking procedures), and be able to handle the vast differences in velocities and accelerations too...

So far, I've already put together things so that the planets and moons and whathaveyou all move correctly in keplerian orbits, and because you only need to input whatever Mean Anomaly is appropriate, they all move accurately even with any time acceleration. (I'm using my own Gametime singleton to manage this). I have a classical-orbital-elements-to-state-vectors mechanic already working just fine, for all kinds of elliptical orbits (including equatorial and circular). Hyperbolic trajectories seems way more difficult, so i'm leaving that out. I can display the orbits of these bodies by sampling Mean anomaly at equidistant points in time, and then displaying the coordinates as points, so the orbit paths are just dotted rings.

The playership can have this-or-that body as its central body, and the mass of that body impacts the physics calculations, naturally (always use the particular g*mass of body, never g directly). The player can change which body it orbits around by exiting or entering their spheres-of-influences (determined by calculating the hill sphere). i'm using a Runge-Kutta integrator for physics, and it seems to work just fine for integrating linear and rotational forces, but (as expected) it becomes less and less accurate the more you accelerate time (timestep increases). That's why the playerShip is intended to be in a deterministic orbit unless the player applies thrust, at which point the orbit will be recalculated. (if the proposed new orbit is not elliptical, then the playerShip "rejects" the proposed thrust and the playership wouldn't change trajectory). (if time is accelerated, then if the player applies thrust the time snaps to "realtime" so the physics would work correctly for applying thrust)

I'd planned on using a state-vector-to-classical-orbit function to be able to update the player's orbit after any change in thrust (then i could display the player's orbit like i do the bodies), but i just can't get it to work "well enough". I've cannibalized a bunch of different libraries out there (poliastro, juliaAstroBase, rene schwarz' memorandi, orbital mechanics for engineers pdf), but every time i wind up getting back an orbit that is (at best) very close to what the "correct" orbit should be after performing the maneuver, but then the player's position (as represented by the nu or mean anomaly) is incorrect. in the best version i have so far, the player's Mean anomaly is 180 degrees off, HALF the time (ie if i have find the orbit once, then it's somewhat accurate, but then i recalculate it again a split-second later and the playership "jumps" 180 degrees in orbit...or 45 degrees, it's not even just "always off by pi")

i've performed the necessary changes to make sure the data i get back fits with godot's frame (orbital calcs usually use Z-positive for up, and godot uses Y-positive, for instance). Again, this all works for classical-to-state (the deterministic orbits work fine). I'm using atan2 to be sure also. Some versions of state-to-classical equations don't work at all accurate, and others are better, but none are "doing it right"...

I'm at a wall here, and have been chewing on this problem for over a month. I know KSP got this kind of mechanic to work, and other similar-to-KSP programs as well (Space Simulator ipad app, SimpleRockets2 also I think).

### Does anyone know how this is done?

To sum up the problem as best i can:

• playerShip position and velocity vectors are relative to current central body.
• playerShip begins in orbit defined by classical orbital elements, and moves by incrementing mean anomaly
• player can apply thrust to their playerShip in any direction (to change orbit)
• ?????? how to get an accurate enough new orbit for playership?

i've had to claw my way up to figuring this all out so far over several months, but i'm worried my fundamental approach of "find orbit from newly changed state vectors" may be intractably wrong, and will never work out, and that there is some other clever way of doing it. can anyone share any insight? i'm happy to share my code if necessary (i'll post whatever snippets you may wonder about, because the entire thing is a bunch), i just didn't want to begin with that because i don't know which area a solution may be found in and this question-post is text-wally enough :|

...and thanks in advance for reading this far!

• I have done something similar in Unity (Gravity Engine asset). My "goto" for algorithms is the code that accompanies Vallado "Fundamentals of Astrodynamics and Applications". Specifically, I use the KEPLER algorithm for time based evolution of an orbit and RVtoCOE for orbit element determination. The code is available for various languages at celestrak.org/software/vallado-sw.php. I avoid evolving based on mean anomaly. Special care is required for near circular orbits, since argument of perigee is not well defined and true anomaly is defined with respect to it. Feb 19 at 15:33
• If you change the orbit, then you change where in the orbit you are. You should expect all its elements to change, including the mean anomaly. You are supposed to get discontinuities. The position doesn't change instantaneously, but the mean anomaly does, as do the other Keplerian elements. Feb 19 at 18:06
• @ryanC - of course the elements change. the problem is that my spaceship is teleporting to roughly the opposite side of the orbit (ie off by 180 degrees). not always the exact opposite, but teleporting nonetheless. this is because after the orbit is found, the ship then has its xyz position found deterministically according to its new orbit, and the anomaly must be off (hence the teleporting). @ nbodyphysics thanks for the tip! I'm going through his code right now, good stuff, fingers crossed. also, what do you evolve using instead of mean anomaly for even time steps? Feb 19 at 19:26
• @w94n9 So use a fictitious value of $\mu$ for these bodies. The problem is that $G$ is only known to less than five digits (the relative uncertainty in the most recently published CODATA value is $2.2\times10^{-5}$). You might as well be using single precision arithmetic if you use $G$. Feb 20 at 4:13
• Here's my standard rant on avoiding G & using mu: astronomy.stackexchange.com/a/48616/16685 Feb 20 at 19:39

## 1 Answer

First, some caveats:

• Always use the gravitational mass parameter $$\mu_\text{body}$$ rather than the product $$GM_\text{body}$$.
The gravitational mass parameter $$\mu_\text{body}$$ conceptually is the product $$GM_\text{body}$$. However, it is $$\mu_\text{body}$$ that is observable while mass is not observable. There are two problems with using $$GM_\text{body}$$: The universal gravitational constant $$G$$ and the mass of the central body. The problem with $$G$$ is that this is one of the least accurately known physical constants. The problem with $$M_\text{body}$$ is that we have yet to devise a scale that can directly measure the mass of an asteroid-sized body, let alone something larger.

• Be wary of the various gotchas using orbital elements.
Argument of periapsis is undefined for circular orbits and right ascension of ascending node is undefined for equatorial orbits. Things get a bit funky for parabolic orbits ($$1/a = 0$$), hyperbolic orbits ($$1/a < 0$$), and radial orbits ($$||\vec h|| = 0$$). There are easy fixes for circular orbits, equatorial orbits, and hyperbolic trajectories which I'll address. This answer will be skipping over the funkiness of parabolic and radial trajectories.

• Be very careful of how you use the two argument arctangent function.
In some programming languages, it's atan2(y,x) (think of it as numerator (y) over denominator (x)) but in others it's atan2(x,y) (because the x coordinate "always" comes before the y coordinate). Read the fine documentation on atan2 for your system.

• Classical orbital elements are only useful in the two body problem.
(Rather, they are only useful for propagation in the two body problem. They are useful for other purposes.) In cases where there are multiple gravitational bodies, or bodies with a non-spherical mass distribution, one has to either resort to using numerical integration, or to using the less accurate patched conic approach.

With those caveats, off we go.

Preliminaries
I'll be using italics for scalars, bold for vectors. For example, $$r = ||\mathbf r|| = \sqrt{\mathbf r \cdot \mathbf r}$$. I'll be using

• $$\mathbf r$$ as the equatorial inertial frame position vector with respect to the central body. The magnitude of this vector is denoted as $$r$$
• $$\mathbf v$$ as the equatorial inertial frame velocity vector with respect to the central body ($$\mathbf v = \frac {d \mathbf r}{dt} = \dot{\mathbf r}$$. The magnitude of this vector is denoted as $$v$$
• $$\mathbf h$$ as the specific angular momentum vector: $$\mathbf h \equiv \mathbf r \times \mathbf v$$, where $$\times$$ denotes the three dimensional cross product (make sure to use the right hand rule). The magnitude of this vector is denoted as $$h$$. Note that the specific angular momentum (and hence its magnitude) is a constant of motion for Keplerian orbits
• Unit vectors constructed from a non-unit vector wear a hat. For example, $$\hat{\mathbf h} = \frac{\mathbf h} h$$. Note however that for historical reasons (hysterical reasons?), this particular unit vector is also denoted as $$\hat{\mathbf w}$$
• $$\mu$$ is the gravitational mass parameter for the central body in question.

Remarks
A radial orbit or trajectory results when $$\mathbf h = \boldsymbol 0$$. Radial trajectories can be bound (total mechanical energy < 0) or unbound (total mechanical energy ≥ 0). These are a royal pain. (Some call these rectilinear orbits or rectilinear trajectories.) Since you are using a digital computer with a standard floating point representation, you probably won't get exactly the zero vector when you compute $$\mathbf r \times \mathbf v$$ for a radial trajectory, but you will get something close: $$h \lll rv$$. Treat such cases as a radial trajectory.

Semi-major axis length $$a$$
This is fairly easy via the vis-viva equation: $$v^2 = \mu\left(\frac2r - \frac1a\right)$$, or $$\frac1a = \frac2r - \frac{v^2}\mu \tag{1}$$ where

• $$a$$ is the semi-major axis length of the current orbit (trajectory)
• $$\mu$$ is the central body's gravitational mass parameter.

Remarks
Equation (1) computes the multiplicative inverse of the semi-major axis length. This inverse can conceptually range from -infinity to +infinity. This creates a problem when $$\frac1a = 0$$: $$a$$ is infinity (or -infinity). This represents a parabolic trajectory. As equation (1) involves a subtraction, on a computer using floating point arithmetic, a value if $$1/a$$ that is close to zero should be treated as zero. Negative values of $$\frac1a$$ (and hence of $$a$$) can also arise. These represent hyperbolic trajectories.

Eccentricity
Here it will be useful to first calculate the eccentricity vector, $$\mathbf e = \frac{\mathbf v \times \mathbf h}{\mu} - \hat{\mathbf r}\tag{2}$$ Another way to write this is $$\mathbf e = \left(\frac1r-\frac1a\right)\mathbf r - \frac{\mathbf r \cdot \mathbf v}{\mu} \mathbf v \tag{2a}$$

The eccentricity vector, like the specific angular momentum vector, is a constant of motion for Keplerian orbits. It is either the zero vector or always points from the central body to periapsis. The eccentricity itself is simply the magnitude of the eccentricity vector: $$e = || \mathbf e || \tag{3}$$ The unit vector defined by $$\mathbf e$$ is oftentimes denoted as $$\hat {\mathbf p}$$ because it points to periapsis. This will be used later.

Remarks
Equation (2) makes it clear what happens in the case of a radial trajectory ($$\mathbf h \approx 0$$). In that case $$\mathbf e = -\hat{\mathbf r}$$, so $$e=1$$. This is one of the two cases where $$e=1$$. The other case is a parabolic trajectory. Equation (2a) makes it clear what happens in a circular orbit as $$1/r - 1/a$$ is zero, as is $$\mathbf r \cdot \mathbf v$$: The eccentricity is zero. This makes $$\hat {\mathbf p}$$ be undefined in the case of circular orbits.

Anomalies
It first helps to define define the perifocal frame as the right-handed triple $$\hat{\mathbf p}, \hat{\mathbf q}, \hat{\mathbf w}$$, where $$\hat{\mathbf q} = \hat{\mathbf w} \times \hat{\mathbf p}$$. Note that the vector $$\hat{\mathbf q}$$ lies on the orbital plane and points from the central body toward 90° true anomaly. Note also that this perifocal frame (aka the PQW frame) is undefined in the case of circular orbits and radial orbits, as $$e = 0$$ and hence $$\hat{\mathbf p}$$ is undefined in the case of circular orbits, and $$\hat{\mathbf w}$$ is undefined in the case of radial orbits. Note also that the PQW frame is well-defined in the case of non-radial parabolic and hyperbolic trajectories.

With this, the true anomaly is easy to find (assuming a non-circular, non-radial orbit): $$\theta = \tan^{-1}\left(\frac {\mathbf r \cdot \hat{\mathbf q}} {\mathbf r \cdot \hat{\mathbf p}} \right)\tag{4}$$ Use the two-argument arctangent function to calculate equation (4). Some like to add $$2\pi$$ to the result of equation (4) when the result is negative. (The result of equation (4) is in the range $$[-\pi,\pi]$$ with most versions of atan2.) I prefer to leave it as-is.

The eccentric anomaly and mean anomaly are given by $$\tan{\frac E2} = \sqrt{\frac{1-e}{1+e}}\tan{\frac\theta2}\tag{5}$$ $$M = E - e\sin E\tag{6}$$ The standard one-argument arctangent function suffices for solving equation (5). Equation (6) is Kepler's equation. You'll need to calculate $$E$$ from $$M$$ via equation (6) when you use Keplerian elements for propagation. Here there is no need for that. Simply apply equation (6) as-is.

Inclination and right ascension of ascending node

The $$\hat{\mathbf w}$$ unit vector in equatorial coordinates is $$\hat{\mathbf w} = \left[\sin i \sin\Omega, -\sin i \cos\Omega,\cos i\right]^T$$ This unit vector has already been calculated, and it provides the information to compute both the inclination and the right ascension of ascending node. Use $$\hat{\mathbf w}_z = \cos i$$ (the $$z$$ equatorial component of $$\hat{\mathbf w}$$, or $$i = \cos^{-1}(\hat{\mathbf w}_z) \tag{8}$$ The $$x$$ and $$y$$ components provide the information needed to calculate the right ascension of ascending node $$\Omega$$ (assuming it is defined): $$\Omega = \tan^{-1}\left(\frac {\phantom{-}\hat{\mathbf w}_x} {-\hat{\mathbf w}_y} \right)\tag{9}$$ Use the two argument inverse tangent, and be sure to negate the $$y$$ component argument.

Remarks
The right ascension of ascending node is undefined in the case of equatorial orbits as $$\sin i = 0$$ when the inclination is 0° or 180°. A widely used approach is to set the right ascension of ascending node to 0 in such cases.

Argument of periapsis
The equatorial $$\hat{\mathbf z}$$ axis expressed in the PQW frame is given by $$\hat{\mathbf z} = \left[\sin i \sin\omega, \sin i \cos\omega,\cos i\right]^T$$ Thus $$\hat{\mathbf z} \cdot \hat{\mathbf p} = \hat{\mathbf p}_z = \sin i \sin\omega$$ and $$\hat{\mathbf z} \cdot \hat{\mathbf q} = \hat{\mathbf q}_z = \sin i \cos\omega$$. With this, the argument of periapsis can be computed via $$\omega = \tan^{-1} \left(\frac {\hat{\mathbf p}_z} {\hat{\mathbf q}_z} \right) \tag{10}$$

Timestamp
One last thing that is needed for propagation purposes is a timestamp, the time at which the position and velocity vectors were calculated (or given). This is typically denoted by $$t_0$$.

Propagation
The only orbital elements that change over time for a Keplerian orbit are the anomalies. The mean anomaly varies linearly with time: $$M(t) = M(t_0) + n(t-t_0)\tag{11}$$ where $$n$$, the mean motion, is given by $$n = \sqrt{\frac {\mu} {a^3}}\tag{12}$$

It is useful to take the result of equation (11) modulo $$2\pi$$. A mean anomaly of $$M+2k\pi$$ where $$k$$ is some integer is the same mean anomaly as $$M$$. Bring $$M$$ into the range $$[0,2\pi)$$ will be beneficial here.

The eccentric anomaly $$E$$ is the solution for $$E$$ for equation (6) given $$e$$ and $$M$$. One way to solve this is via Newton-Raphson iteration. The Newton-Raphson iterator for solving equation (6) is $$E_{n+1} = E_n - \frac{E_n-e\sin E_n - M}{1-e\cos E_n} = \frac{M+e(\sin E_n - E\cos E_n)}{1-e\cos E_n}$$ An initial guess is needed. Many people use $$E_0 = M$$, which is a good guess for low eccentricity orbits. It is not a good guess for moderate or high eccentricity orbits. An eccentricity of 0.6 is enough to send this Newton-Raphson iterator into orbit. A better (more stable) initial guess is $$E_0 = \pi$$ -- assuming $$M$$ has been placed in the range $$[0, 2\pi)$$. This works for all but the highest eccentricity orbits. If you're going to be having orbits within a small epsilon of 1 you will need to use something more stable but slower than Newton-Raphson.

Once the eccentric anomaly has been found, equation (5) can be used to find the true anomaly. No iteration is needed for this final part of the problem.

Hyperbolic trajectories
For a non-radial hyperbolic trajectory, the semi-major axis length by equation (1) is negative and the eccentricity via equation (3) is greater than one. This would lead to imaginary numbers for the eccentric anomaly and for the mean motion if used as-is. The workaround is fairly simple: Use hyperbolic functions where needed. The hyperbolic equivalent of equations (5), (6), and (12) are $$\tanh\frac H2 = \sqrt{\frac{e-1}{e+1}} \tan\frac\theta2 \tag{5H}$$ $$M = e\sinh H-H\tag{6H}$$ $$n = \sqrt{\frac{\mu}{-a^3}} \tag{12H}$$

• Thanks for such a great explanation! it's very helpful for checking the logic throughout. vallado's code also helped considerably. i think i MAY have found the problem too...namely, for undefined orbits (as has been pointed out) returned values for whatever equations involved must be kept in mind, and then later on appropriate checks must be inserted. i think what was happening was that an edge-case orbit was returning a wonky value, and then this was throwing off the m or raan. i've cleaned things up CONSIDERABLY, and it seems good now...MANY MANY MANY THANKS! :D i was going nuts lol Feb 20 at 17:59
• also, i use REALLY Big timesteps in time acceleration, and the Kepler algorithm from Vallado will eventually break (at steps like delta=1 month, crazy fast) and throw the satellite waaaay out of orbit (or worse). using mean anomaly to advance bodies in orbit keeps things perfectly tight, hence why it was so important to get that rv2coe algo working for any orbital changes (and it SEEMS like it is...i'll try to break it of course :P but g*ddam what a catharsis that its spitting out accurate orbits now!! again, much obliged for your insight, gents, you've really helped lower my stress levels:) Feb 20 at 18:04