How can I predict what an object's orbital state vectors will be in the future?

Question

I have a position and velocity vector for an orbiting object and can calculate the orbital elements. From a simple two-body situation, how can I predict where in the orbit the object will be at some future time? What are the mathematical steps involved?

Answer 1

There are two ways. One is to numerically integrate from the current time to the future time. That is conceptually the easiest, as well as the simplest to implement (given a good integration routine). On the other hand, it is the most compute intensive, most prone to accumulated numerical error depending on how far you have to go, and provides no insight into the orbit.

For situations that are not your two-body problem, that is really the only option.

The second is to derive the six orbital elements from the six components of position and velocity (assuming some convention for $t=0$). That tells you a great deal about the future of the object, such as the size, shape, and orientation of the orbit, and the orbit period. From the elements you can compute the position and velocity at any future time with a few simple formulas, without having to integrate through the times in between. Again, assuming your two-body problem.

To answer the question in the comments, the position and velocity can be computed using the orbital elements thusly. In the plane of the orbit, where $\mu$ is the $GM$, $a$ is the semi-major axis, $e$ is the eccentricity, and $\tau$ is the eccentric anomaly (more on that one in a bit):

$x=a\left(\cos\tau-e\right)$

$y=a\sqrt{1-e^2}\sin\tau$

$z=0$

$v_x=-\sqrt{\mu\over a}{\sin\tau\over 1-e\cos\tau}$

$v_y=\sqrt{\mu\over a}{\sqrt{1-e^2}\cos\tau\over 1-e\cos\tau}$

$v_z=0$

You can then rotate those coordinates to the actual plane of the orbit by applying an Euler rotation with the angles $\Omega$, $i$, and $\omega$, where those are the longitude of the ascending node, the inclination, and the argument of periapsis respectively.

The eccentric anomaly goes from $0$ to $2\pi$ over one orbit, and is a convenient substitute for time. It is related to time by:

$t=\sqrt{a^3\over\mu}\left(\tau-e\sin\tau\right)$

There is no closed form solution to get $\tau$ from $t$, so you need to solve that numerically. Or if you are, for example, plotting, you can calculate the coordinates and $t$ all from $\tau$ and plot parametrically.

This all assumes the convention that $t=0$ and $\tau=0$ is at periapsis. You will need to offset $t$ for your solution, since I expect that you will want $t=0$ to be at your initial condition, which is probably not at periapsis. That is where the sixth orbital element for the anomaly comes into play.

This will give you the answer in a fake coordinate system for the two-body problem that is pretty darned close to the real coordinate system if the body being orbited is much, much heavier than the object orbiting. However if that's not the case, then you need to convert the coordinates back to the real situation for both bodies, since they are both orbiting their combined center of mass. An example of that case is the Moon orbiting the Earth.

Answer 2

@Stu, I'll add some information from the equations of motion perspective. These comments are more geared towards artificial satellite orbiting Earth problem--additional bodies and perturbations throws a lot of this out the window, but some fundamentals remain. Apologies in advance if you already know this.

We make some assumptions on the relative masses of the two-bodies involved, and assume the central body has a far greater mass. We also assume an inertial reference frame, and that the bodies are point masses. Using these we can develop the 2BP equation of motion where $r$ is the position vector from the central body center of mass, to the second body

$$ \ddot{\vec{r}}=-\frac{\mu}{r^2} \frac{\vec{r}}{r} $$

Now, the fact that we have an $r^3$ in the denominator shouldn't worry you - This is the inverse square gravity law, and the additional $r$ scales $\vec{r}$ to produce a unit vector. This looks a little nasty to solve in closed form, since we have a non-linear, ordinary differential equation (which is also coupled due to the $r$ term). There are 3 second-order ODEs, which can be expressed as 6 first-order ODEs. If you have access to a numerical integrator, all you need now are initial conditions on your position and velocity vectors (the 6 states), set your desired time, and away you go. You can propagate this problem backward in time as easily as forward in time.

With respect to orbital element and position velocity conversions, the unperturbed two-body problem essentially has one degree of freedom, that being your choice of anomaly in the plane. With one anomaly specified, you can calculate any of the other two. Given the mean anomaly $M$, you can iterate to find the eccentric anomaly $E$, and from there determine the true anomaly $\nu$ (with care to avoid quadrant ambiguity with the arctangent function). $$ M=E-e\sin{E} $$

$$ \tan{\frac{\nu}{2}} = \sqrt{\frac{1+e}{1-e}}\tan{\frac{E}{2}} $$

Propagating the anomalies allows one to more easily visualize the orbit, versus observing position/velocity vectors. There exist numerous open-source routines to easily convert (typically named rv2oe or rv2coe) between position/velocity and orbital elements. This is a good website for Vallado's collection of software.