The position of the Moon and of the planets

First, compute the eccentric anomaly, E, from M, the mean anomaly, and e, the eccentricity. As a first approximation, do (E and M in degrees):

E = M + e*(180/pi) * sin(M) * ( 1.0 + e * cos(M) )

or, if E and M are in radians:

E = M + e * sin(M) * ( 1.0 + e * cos(M) )

If e, the eccentricity, is less than about 0.05-0.06, this approximation is sufficiently accurate. If the eccentricity is larger, set E0=E and then use this iteration formula (E and M in degrees):

E1 = E0 - ( E0 - e*(180/pi) * sin(E0) - M ) / ( 1 - e * cos(E0) )

or (E and M in radians):

E1 = E0 - ( E0 - e * sin(E0) - M ) / ( 1 - e * cos(E0) )

For each new iteration, replace E0 with E1. Iterate until E0 and E1 are sufficiently close together (about 0.001 degrees). For comet orbits with eccentricites close to one, a difference of less than 1E-4 or 1E-5 degrees should be required.

If this iteration formula won't converge, the eccentricity is probably too close to one. Then you should instead use the formulae for near-parabolic or parabolic orbits.

Now compute the planet's distance and true anomaly:

xv = r * cos(v) = a * ( cos(E) - e )
yv = r * sin(v) = a * ( sqrt(1.0 - e*e) * sin(E) )

v = atan2( yv, xv )
r = sqrt( xv*xv + yv*yv )