This code is obscure. It has all to do with Kepler’s discovery that orbits are ellipses and sweep equal areas in equal times. In the routine Kepler we take an eccentricity e and local ‘mean anomaly’ L and return Cartesian position and velocity. The coordinates are in a Cartesian frame private to the planet where perihelion is at θ=0 and r=1/(1+e). The parameter L is also the local ‘mean anomaly’ to which ϖ must be added to get the conventional mean anomaly.

See this about the routine Kepler.

agm is angular momentum which should be and is preserved. td is total distance to foci which should be and is preserved.

This code (with logic) computes such orbits numerically starting from Newton’s differential equations. The code has been duplicated here and further developed in order to use the numerical solution to differential equations.

We use planetary masses from here.

Points of interest along the x axis for ellipse r = 1/(1+e∙cosθ): (0≤e<1)

−1/(1−e): Left end of ellipse. Aphelion. θ=π.
−2e/(1−e²): Alternate focus.
−e/(1−e²): Center of ellipse.
0: Origin. Primary. Sun. Focus.
1/(1+e): Right end of ellipse. Perihelion. θ=0.

Some references: How to compute planetary positions
Computing planetary positions
Celestial coordinate system