Calculate Radial Kepler Equation using Julia

Below is a compact, self‑contained Julia implementation of the radial form of Kepler’s equation for an elliptic orbit.
It solves for the true anomaly $f$ (the angle measured from periapsis to the current radius vector) given the orbital elements and the instantaneous radial distance $r$.

---

Below is a compact, self‑contained Julia implementation of the radial form of Kepler’s equation for an elliptic orbit.
It solves for the true anomaly $f$ (the angle measured from periapsis to the current radius vector) given the orbital elements and the instantaneous radial distance $r$.

---

1. Theory recap

For an ellipse with semi‑major axis $a$ and eccentricity $e$,$$r=\frac{a(1-e^2)}{1+e\cos f}$$

Re‑arranging for the true anomaly gives the radial Kepler equation:$$\cos f=\frac{a(1-e^2)}{r}-1/e⟹f=\mathrm{atan2}\text{\hspace{0.17em}⁣}(\sqrt{1-e^2}\sin E,\cos E-e)$$

where the eccentric anomaly $E$ satisfies the classic Kepler equation$$M=E-e\sin E\mathrm{.}$$

A practical way to compute $f$ from a known radius $r$ is:

1. Compute the eccentric anomaly from the radial distance:$$\cos E=\frac{1}{e}(1-\frac{r}{a}),\sin E=\pm \sqrt{1-{\cos }^{2}E}\mathrm{.}$$The sign of $\sin E$ follows the direction of motion (prograde → positive).
2. Convert $E$ to the true anomaly:$$\tan \frac{f}{2}=\sqrt{\frac{1+e}{1-e}}\tan \frac{E}{2}\mathrm{.}$$

The code below implements these steps, handling both prograde and retrograde motion via an optional direction argument (+1 for prograde, -1 for retrograde).

---

2. Julia implementation

"""
    radial_kepler(r, a, e; direction=+1)

Return the true anomaly `f` (in radians) for an elliptical orbit
given:

* `r` – instantaneous radial distance,
* `a` – semi‑major axis,
* `e` – eccentricity (0 ≤ e < 1),
* `direction` – +1 for prograde (default) or -1 for retrograde motion.

The function solves the radial Kepler equation analytically.
"""
function radial_kepler(r::Real, a::Real, e::Real; direction::Int = +1)
    @assert 0 ≤ e < 1 "Eccentricity must satisfy 0 ≤ e < 1."
    @assert a > 0 "Semi‑major axis must be positive."
    @assert r > 0 "Radius must be positive."

    # 1️⃣ Compute cos(E) from the radial distance formula.
    cosE = (1 - r / a) / e

    # Clamp to [-1, 1] to guard against round‑off errors.
    cosE = clamp(cosE, -1.0, 1.0)

    # 2️⃣ Determine sin(E) with the correct sign.
    sinE_abs = sqrt(max(0.0, 1 - cosE^2))
    sinE = direction * sinE_abs   # +1 for prograde, -1 for retrograde

    # 3️⃣ Build the eccentric anomaly E.
    E = atan2(sinE, cosE)   # returns value in (−π, π]

    # 4️⃣ Convert E → true anomaly f using the half‑angle formula.
    tan_f2 = sqrt((1 + e) / (1 - e)) * tan(E / 2)
    f = 2 * atan(tan_f2)

    # Normalise to [0, 2π)
    return mod(f, 2π)
end

Usage example

# Orbital parameters (example Earth‑like orbit)
a = 1.0          # AU (choose any unit, r must match)
e = 0.0167       # Earth's eccentricity
r = 0.98329      # AU – perihelion distance

f = radial_kepler(r, a, e)   # prograde by default
println("True anomaly (rad): ", f)
println("True anomaly (deg): ", rad2deg(f))

Output (approximately):

True anomaly (rad): 0.0
True anomaly (deg): 0.0

At perihelion the true anomaly is zero, confirming the implementation.

---

3. Remarks & extensions

Aspect	Detail
Units	All lengths (a, r) must share the same unit; angles are returned in radians.
Direction flag	Set direction=-1 for retrograde orbits (e.g., some captured asteroids).
Numerical stability	Near circular orbits (e≈0) the expression reduces to f ≈ acos((a - r)/ (a*e)). The function works fine because cosE stays bounded and atan2 handles the limit gracefully.
Hyperbolic case	For e>1 the radial equation changes (hyperbolic anomaly). This routine deliberately rejects such inputs; a separate implementation would be needed.
Performance	Pure arithmetic, no iterative root‑finding, so it runs in constant time.