Applying Kepler’s rules is tricky. Many years ago a teacher whose name I do not recall told us that people who worked in this area were highly conventional in their notation. They stuck to a rigid notation for names of variables and even what quantities to compute. The conventional way is confusing but variations make it worse. I am beginning to understand that now. There are many ways to do this calculation but few advantages among them. I think the stuff below is slightly unconventional and that adds confusion to an already confusing situation. I think it is right but I shall try to make it conventional without introducing errors.

We consider a central attractor at the origin whose mass is 1/(gravitational constant). Actually we will assume units where that constant is 1 and so the central mass is 1 and that mass causes another small mass at (x, y) to be accelerated toward the origin with magnitude 1/(xx+yy). Said otherwise the gravitational potential is −1/r where r is the distance to the origin. The period of a circular orbit with radius 1 is 2π.

We consider the set of orbits given by the formula:
r = 1/(1 + e∙cos θ)
where e is the eccentricity. This curve is an ellipse with a focus at the origin.

With eccentricity e, the semi-major axis length a = 1/(1−e2) and the semi-minor axis length (ys) = 1/√(1−e2). To achieve equal area in equal time (equal ΔL) imagine a unit circle in a new plane which is an affine transform of the orbit. The center of the ellipse goes to the origin of the new plane. The perihelion and aphelion go respectively to (1, 0) and (−1, 0). The foci of the ellipse go to (e, 0) and (−e, 0). Horizontal scale factor is 1/a = 1−e2 and vertical factor is 1/√a = 1/(ys). A source of confusion is that the total area to be swept out depends on e among our set of orbits. An affine transform leaves areas proportional and it is easier to compute areas in the circle. The Sun goes to (e, 0). We need to consider the area swept out by the radius vector to (cosφ, sinφ) but deducting the area of the triangle spanned by (0, 0), (e, 0), (cosφ, sinφ). The respective areas are φ/2 and e∙(sinφ)/2. The triangle’s area is negative when π<φ<2π. The circle in the new plane has area π but the mean anomaly L which is an argument to Kepler runs from 0 to 2π for one traversal of the entire orbit. (That factor of 2 cost me a day!) We need to find φ so that L = φ − e∙sinφ. This is uniquely possible when −1 < e < 1 and 0 ≤ L < 2π. Then 0 ≤ φ < 2π. L and φ are monotonic increasing functions of each other. Likewise the relationship between angles θ and φ. The point (cosφ, sinφ) transforms back to the point on the orbit matching the given mean anomaly L.

In the code the variable q converges to a triple (cosφ, sinφ, φ) so that L = φ − e∙sinφ. This φ is not the θ from 1/(1+e∙cosθ)! The code is in a convenient position to return the velocities as well.

There are more efficient ways to solve for φ but we don’t spend much time here.

The 2nd argument to routine Kepler is named L which conventionally refers to ‘mean longitude’. To compute cartesian positions and velocities from ‘the elements’ the 2nd argument is actually L−ϖ where ϖ is the ‘longitude of perihelion’. The velocities that Kepler returns are scaled as if they are in the plane where we found the original orbit: r = 1/(1 + e∙cos θ). Kepler(e, L) = Kepler(e, L + 2π). When e>0 L is not t+constant because the period of the orbit is 2π(1−e2)−3/2 ≠ 2π. The velocity that Kepler returns (Kepler(e, L).xd) reflects this longer period. The caller may need to scale the returned velocity up accordingly.