Given an orbit's semi-major axis a, eccentricity e, and a radial distance r, you can recover the true anomaly: where the body sits along its orbit. Here's the compact Julia routine I use.

The math

From the orbit equation, the eccentric anomaly E satisfies cos E = (a - r) / (a·e), and the true anomaly f follows via the half-angle relation. Wrapping it in atan keeps it stable across quadrants.

Two things in that line are doing more work than they look. The first is that acos returns a value in [0, π] only, so it cannot tell you whether the body is climbing from perihelion to aphelion or falling back. Radius alone is ambiguous: every distance between perihelion and aphelion is visited twice per orbit. That's what the direction argument is for, and it is the caller's job to know the sign of the radial velocity.

The second is the clamp. Analytically (a - r)/(a·e) lives in [-1, 1], but feed it an r computed in floating point at exactly perihelion and it lands on 1.0000000000000002, and acos returns a domain error instead of zero. Clamping is not sloppiness; it's the correct handling for a quantity whose bounds are guaranteed by the maths and not by IEEE 754.

The code

function radial_kepler(r, a, e; direction=1)
    e >= 1 && error("Only elliptical orbits (e < 1) supported")
    cosE = clamp((a - r) / (a * e), -1.0, 1.0)
    E = acos(cosE)
    f = 2 * atan(sqrt((1 + e) / (1 - e)) * tan(E / 2))
    return direction * f
end

a, e, r = 1.0, 0.1, 0.9        # r = a(1 - e): perihelion
println("True anomaly (deg): ", rad2deg(radial_kepler(r, a, e)))

At perihelion the true anomaly is zero, exactly what comes back. A nice sanity check.

Note that a circular orbit divides by zero here, since e = 0 puts a·e in the denominator. That's not a bug so much as a statement that the question is meaningless: on a circle every point is at the same radius and there is no perihelion to measure from.

Why Julia

For math-heavy numeric work Julia hits a sweet spot: the syntax reads like the equations, and it runs at near-C speed without dropping into another language. The reason it can is that the compiler specializes each function on the concrete argument types at the first call and emits native code through LLVM. Call radial_kepler with Float64 arguments and you get a tight float routine; call it with a dual number from an autodiff package and you get a compiled derivative of the same source, with no changes to it.