Statement of the Problem
The problem you want to solve is called the Kepler problem. In your formulation of the problem, you're starting out with the Cartesian orbital state vectors (also called Cartesian elements): that is, the initial position and velocity.
As you have discovered, the only way to propagate the Cartesian elements forward in time is by numerical integration. This works OK, but it can be slow if you want high accuracy, and there are some numerical problems (errors caused by rounding[slowly accumulate, and many integrators cause energy drift). You can get around some of these problems by using a higher-order integrator (Runge-Kutta is a popular one) which allows you to take larger steps for the same level of accuracy, or get better accuracy for the same step size. However, this is kind of overkill for simple simulation.
If your simulation can be treated as a two-body problem, then things simplify dramatically. The two-body problem is a good simplification if the simulation objects are mainly affected by a single, large object. For example, the Earth traveling around the Sun or a spacecraft traveling in a low Earth orbit are well-modeled as a two body problem; however, a spacecraft traveling from the Earth to the Moon is not (more on that later).
Since you're trying to model the positions of the planets with medium accuracy, reduction to the two-body problem should work for you.
Definition of Terms
The traditional solution to the two-body problem involves a different way of representing the position of the orbiting body, called the Keplerian orbital elements (also just called orbital elements). Instead of specifying position and speed, they specify six different parameters of the orbit (if you just want to get to the code, you can skip this part):
Semi-major axis, $a$: Half the maximum diameter of the ellipitical orbit, ( = cicle radius if the orbit is circular). The energy and period of the orbit depend only on $a$. The semi-latus rectum $\ell$, the "width" of the orbit, can be a better choice for orbits that are close to parabolic (like for asteroids) or that change from elliptical to hyperbolic (like interplanetary spacecraft). The two are related by $\ell=a(1-e^2)$.
Eccentricity, $e$: The "pointiness* of the orbit. Ranges from $e=0$ for a perfectly circular orbit, to $e=1$ for a parabolic orbit, to $e>1$ for hyperbolic orbits. Mercury is the most eccentric planet with $e \approx 0.2$. Earth-orbiting spacecraft usually have $e<0.01$.
Aside: From $e$ and $a$ we can determine the most far and the closest points in the orbit, the apoapsis and periapsis (together apsides):
$$ r_a = a(1+e) \\ r_p = a(1-e) $$
The naming of these points is a little funny: apoapsis and periapsis are the generic terms, but orbits around particular bodies have specific terms: a spacecraft around the Earth has an apogee and perigee, while the Earth (in orbit around the Sun) has an aphelion and perihelion.
The two parameters $a$ and $e$ are enough to determine the shape of the orbit. The next three parameters define the orientation of the orbit relative to a coordinate system consisting of a reference plane, and a reference direction (parallel to the plane).
For almost all orbits in the solar system, the coordinate system used is the ecliptic coordinate system. The reference plane is the ecliptic plane, the plane of the Earth's orbit around the Sun. The reference direction is the vernal equinox point, the direction from the Earth to the Sun at the moment of the spring equinox. Since both these references drift slowly over time, we must specify a particular time at which these references are defined, called the epoch. The most common is J2000, noon on January 1, 2000 (UTC).
Earth-centered orbits often use the equatorial coordinate system, whose reference plane is the equator of the Earth. The situation with the epoch is a little complicated, so I won't get into it here.
Following parameters locate the orbit w.r.t. Earth orbit:
Inclination, $i$: the angle between the plane of the orbit and the reference plane. Inclination between 90 and 180 degrees refers to a retrograde orbit, one that orbits "backwards" from the usual direction.
Longitude of the ascending node, $\Omega$: the ascending node is where the orbit crosses from below the reference plane to above it. (It's at the intersection between the orbital plane and the reference plane) $\Omega$ is the angle between this point and the reference direction, measured counterclockwise.
Argument of periapsis, $\omega$: the angle between the ascending node and periapsis (the lowest point in the orbit). For orbits with very low inclination where the location of the ascending node is difficult to determine (since it's the intersection between two almost parallel planes), we instead use the longitude of periapsis $\varpi = \Omega + \omega$.
The sixth parameter defines the position of the object in its orbit. There are a couple different choices, but the most common is:
- Mean anomaly, $M$: an "imaginary" angle that is zero at periapsis and increases at a constant rate of 360 degrees per orbit.
The rate at which $M$ changes is called the mean motion, $n$, equal to $2\pi/T$. Usually you have a measurement of $M$ at a particular epoch $t_0$, called (unsurprisingly) the mean anomaly at epoch, $M_0$.
Just like the argument of periapsis, for low-inclination orbits we use a related value, the mean longitude, $L=\varpi + M$.
The actual angle between the orbiting body and periapsis is called the true anomaly, $\nu$. This is the angle we need to compute the position of the body. Unfortunately there is no way to directly compute $\nu$ from $M$. Instead we first solve for the eccentric anomaly $E$:
$$
M = E - e \sin E
$$
This is called Kepler's equation, and it cannot be solved analytically. Once we have $E$ though, there is a relatively simple expression for $\nu$.
Computing Position from Orbital Elements
We'll perform this computation in three steps: first, we'll solve Kepler's equation. Second, we'll compute the 2d position of the body in the orbital plane. Lastly, we'll rotate our 2d position into 3d coordinates. I'll give some "pseudocode" in Javascript for most of these tasks.
I'll assume that you're using a set of elements like these from JPL's web site. These use $L$ and $\varpi$ instead of $M$ and $\omega$. The table gives two values for each of the elements; the second is the time derivative. If you use the values in this table you should use the derivatives as well.
Calculate the time $t$ in centuries from J2000:
// month is zero-indexed, so 0 is January
var tMillisFromJ2000 = Date.now() - Date.UTC(2000, 0, 1, 12, 0, 0);
var tCenturiesFromJ2000 = tMillisFromJ2000 / (1000*60*60*24*365.25*100);
Now we calculate the current values of each of the orbital parameters. For example, the semimajor axis of Earth, using the values from Table 1 (valid from 1800–2500):
// a0 = 1.00000261; adot = 0.00000562
var a = a0 + adot * tCenturiesFromJ2000;
(Note that the values are actually given for "EM Barycenter," the center-of-mass of the Earth-Moon system. The Earth is around 4600 kilometers from the barycenter in the opposite direction from the Moon. If you want to correct this inaccuracy you'll need to simulate the motion of the Moon as well, but that's probably overkill.)
Table 2a gives elements that are accurate from 3000 BC to 3000 AD; however, if you use the elements from table 2a, you must supplement them with corrections to $L$ from Table 2b! For example, here is computing the longitude of Saturn:
// L0 = 34.33479152; Ldot = 3034.90371757
// b = -0.00012452
// c = 0.06064060
// s = -0.35635438
// f = 38.35125000
var L = L0 + Ldot * tCenturiesFromJ2000
+ b * Math.pow(tCenturiesFromJ2000, 2)
+ c * Math.cos(f * tCenturiesFromJ2000)
+ s * Math.sin(f * tCenturiesFromJ2000);
We don't need to explicitly compute the mean motion and add it to $L$, since both tables include it in $\dot L$.
Now we're ready to compute $M$ and $\omega$ (w):
var M = L - p \\ p is the longitude of periapsis
var w = p - W \\ W is the longitude of the ascending node
On to step 2: we need to solve the Kepler equation:
$$
M = E - e \sin E
$$
We can solve this numerically using Newton's method. Solving the Kepler equation is equivalent to finding the roots of $f(E) = E - e \sin E - M$. Given $E_i$, an estimate of $E$, we can use Newton's method to find a better estimate:
$$
E_{i+1} = E_i - f(E_i) / f'(E_i) \\
f'(E) = 1 - e \cos E
$$
Since the nonlinear part $e \sin E$ is very small, we can start with the estimate $E=M$. Our code looks something like this:
E = M;
while(true) {
var dE = (E - e * Math.sin(E) - M)/(1 - e * Math.cos(E));
E -= dE;
if( Math.abs(dE) < 1e-6 ) break;
}
Now there are two ways to compute the position from the eccentric anomaly. We can first compute the true anomaly and radius (the position of the object in polar coordinates), and then convert to rectangular coordinates; however, if we apply a bit of geometry we can instead compute the coordinates directly from $E$:
var P = a * (Math.cos(E) - e);
var Q = a * Math.sin(E) * Math.sqrt(1 - Math.pow(e, 2));
(P and Q form a 2d coordinate system in the plane of the orbit, with +P pointing towards periapsis.)
Finally, we can rotate these coordinates into the full 3d coordinate system:
// rotate by argument of periapsis
var x = Math.cos(w) * P - Math.sin(w) * Q;
var y = Math.sin(w) * P + Math.cos(w) * Q;
// rotate by inclination
var z = Math.sin(i) * y;
y = Math.cos(i) * y;
// rotate by longitude of ascending node
var xtemp = x;
x = Math.cos(W) * xtemp - Math.sin(W) * y;
y = Math.sin(W) * xtemp + Math.cos(W) * y;
(x, y, and z will be in units of AU.)
And you're done!
A few tips:
If you want to calculate the velocity as well, you can do it at the same time as you calculate $P$ and $Q$, then rotate it in the same way.
$$
\dot M = n = \dot L \\
\dot M = \dot E - e (\cos E) \dot E \\
\dot E = \dot M / (1 - e \cos E) \\
\dot P = -a (\sin E) \dot E \qquad \dot Q = a (\cos E) \dot E \sqrt{1 - e^2}
$$
Note I don't include any of the derivatives (except $\dot L$) in this calculation, since they don't affect the outcome much. You could code this as:
var vP = - a * Math.sin(E) * Ldot / (1 - e * Math.cos(E));
var vQ = a * Math.cos(E) * Math.sqrt(1 - e*e) * Ldot / (1 - e * Math.cos(E));
Note that the velocities will be in AU per century.
If you are updating the positions very frequently, you could use the previous value of $E$ to seed Newton's method, and do a fixed number of iterations (probably just one would suffice). Note however that you need to keep that value of $E$ local to each object!
You can also just use a fixed number of iterations for the initial solution. Even for $e=0.2$, after three iterations the error in $E$ is only about $10^{-13}$, and after four iterations the error is smaller than the rounding error of an IEEE double up to $e=0.42$.
If you want more information you can search online, but if you're really interested you should read an introductory text on orbital mechanics. I personally recommend Fundamentals of Astrodynamics by Bate, Mueller, and White (pdf). My dad used this book back when he was in college, and I found it to be more readable than my college textbook. You'd be interested in Chapter 4, Position and Velocity as a Function of Time.
