Analytic solution to the Kepler problem from position + velocity initial conditions

Question

I am writing a javascript program (web page) that uses iterative simulation to show the motion and mechanics of a satellite in orbit around the earth.  So far, using only circular orbits, it has been working perfectly.  However, I am now adding arbitrary velocities to these satellites to see how that affects their motion and orbits and they do not seem to be moving in the ways that I expected.  My problem is that I do not know if this is a mistake in my programming or in my expectations.  So since I know that 2-body and simpler problems in gravity are solved/solvable, I figured that I would just look up the functions for this.
I have now spent the better part of the last two days trying to find the answer to what I believed must be a common question with a simple and straight-forward answer (I still cannot understand why this is so hard to find).  That is, what are the functions of position and velocity for an satellite (object with negligible mass) in orbit around a primary (unmoving) body at a specific point in time?  This can be in either two or three dimensions, or in polar coordinates.  For cartesian coordinates in two dimensions (preferred), the functions would have the following forms, either:
In vector/matrix form:

Xt = FX(m,X0,V0,t)
Vt = FV(m,X0,V0,t)

Where:
m       - the mass of the primary
t       - the current time
X0, V0  - satellite’s initial position and velocity vectors
Xt, Vt  - satellite’s the position and velocity at Time=t

I know of course that FV = FX’, so I can get by with just FV.  I can also work with these in parametric form and/or in polar coordinates as well.
As I mentioned above, I have spent a lot of time trying to find this, at more sites and Wikipedia pages than I can remember at this point (here are two: https://en.wikipedia.org/wiki/Orbital_mechanics, and https://en.wikipedia.org/wiki/Classical_central-force_problem#Specific_solutions) and also some questions on this site such as Using 2D position, velocity, and mass to determine the parametric position equations for an orbiting body, which specifically do not answer the question that I am asking.
Because of my frustration with all of these sites I want to be very clear what I am and am not looking for:  I am looking for an analytical, time-dependent function (or system of functions) that I can code into a program (or spreadsheet or calculator, etc.) that will tell me the position and velocity of a satellite around a primary at a given point in time.  Assuming that the primary is an ideal point-mass that doesn’t move (i.e., a Central Force), that the satellite’s mass and size are irrelevant and using idealized Newtonian gravity.
It is possible that the articles have such a solution somewhere in them, but I sure haven’t been able to find them.  For some reason, every time I search for “solution” in these long articles, I come up with things that look more like problems than actual solutions.  So, by analytical I mean, no integrals, derivatives, partial derivatives, differential equations, infinite series or other type of numerical solutions (technically, what I am already doing is a numerical solution, but I need to check it for accuracy).  I also am not looking for solutions in terms of other parameters such as eccentricity, major/minor axis, or angular momentum, as I do not have those values, what I have is the satellite's current position and velocity.
If for some reason it is not possible to provide exact time-dependent functions for these (and it really seems there should be because I have always heard that Two-Body and simpler problems are solvable), then I would accept 1) some reference/explanation for this when it's supposed to be a solved problem, and 2) a time-independent solution with equations of the form:
FX(m,X0,V0,Xt) = 0
FV(m,X0,V0,Vt) = 0

Again, parametric forms and/or polar coordinates would also work for me.  Although a time-dependent solution would be vastly preferred, I can at least use time-independent equations to check how far my simulation has deviated off the correct path/track.

Since the [computational-physics] and [software] tags have been added, I wanted to note a couple of things, in case they weren't already clear:

I am not a physicist, I am a programmer with a Math/CS background
Nonetheless, I do not think that this is a software/computational question, I included that information as background to explain what form of answer I was looking for.

FYI, for those who may be interested, I have put up this (still in development) website at http://www.rbarryyoung.com/EarthOrbitalSimulator.html.  The issue that prompted my concern is the motion of the missile after you press the Fire button.  It seems a bit too flattened to me, plus it's loop does not actually return to the orbiter, but slowly gets farther and farther away from it.  However, I am not sure if this is right or wrong.  (Also, there's no mobile layout for this page yet, so I don't know if it will work on phones or not)

Kyle Oman · Answer

The Kepler potential (i.e. for a point mass) is $Phi=-frac{GM}{r}$. It's convenient to make a simple change of variables $u=1/r$. The (specific) angular momentum is conserved along the orbit, it's value is easily obtained from your initial conditions as $L=|vec{r}times vec{v}|$. Note that since the orbit is in a plane, you can work in a 2D coordinate system, I'll use polar coordinates $(r,psi)$. You can always rotate the solution later.
The equation of motion in this potential is:
$$frac{{rm d}^2u}{{rm d}psi^2}+u=frac{GM}{L^2}.$$
The general solution of this equation is:
$$u(psi)=Ccos(psi-psi_0)+frac{GM}{L^2}.$$
$M$ is the mass of your point mass, and $L$ is the specific angular momentum defined above, which leaves two constants $C>0$ and $psi_0$ to be determined.
For the time dependence, solving $L=r^2dot{psi}$ gives:
$$t=frac{T_r}{2pi}(eta-esineta)$$
with $e=frac{CL^2}{GM}$, $eta$ is defined by the relationship $sqrt{1-e}tanleft(frac{1}{2}(psi-psi_0)right)=sqrt{1+e}tanleft(frac{1}{2}etaright)$, and $T_r=frac{a^2}{L}sqrt{1-e^2}$. Finally, $a=frac{L^2}{GM(1-e^2)}$.
Note that $t=0$ is defined to be at the pericentre (when the distance is $a(1-e)$). With all that, and your chosen "initial" distance and velocity vector, you should be able to determine all the constants and the time-evolution of the orbit. Easy enough, I guess...
If you're interested in unbound (parabolic & hyperbolic) orbits, it gets slightly more complicated.
Reference for all the above is Sec. 3.1 in Binney & Tremaine (2008), from which I've also borrowed the notation and all the equations.

Analytic solution to the Kepler problem from position + velocity initial conditions

One Answer

Add your own answers!

Ask a Question