Dado r/a, ¿cuáles son los límites en la dirección en que un cuerpo en órbita podría moverse (por ejemplo, ángulo sólido frente a r/a)?

No hay límites en la dirección.

La ecuación Vis-viva te dará una velocidad. Suponiendo masas puntuales y siguiendo la mecánica clásica, a la ecuación de Vis-Viva no le importa en absoluto en qué dirección apunta su velocidad; Es simplemente una ecuación basada en cómo la energía orbital total (que es la misma para todas las órbitas con el mismo semieje mayor alrededor del mismo cuerpo) debe distribuirse entre la energía potencial gravitatoria y la energía cinética.

Para las órbitas keplerianas, las únicas restricciones sobre$r$y$a$son:

• $r$será un valor positivo.
• $a$debe ser distinto de cero.
• Si$a$es positivo (lo que significa una órbita elíptica),$r$nunca excederá$2a$(Si$r$=$2a$, estás viendo la apoapsis de la elipse lineal degenerada)
• Si$a$es negativo (lo que significa una trayectoria hiperbólica),$r$puede ser cualquier valor positivo que elija.

Dicho de otro modo, por la ecuación vis-viva, dada una distancia radial$r$y un semieje mayor$a$alrededor de un cuerpo gravitatorio define un valor de velocidad orbital$v$. En las condiciones newtonianas ideales de dos cuerpos, independientemente de la dirección en la que apunte, esa velocidad$v$, siempre estarás en una Órbita/Trayectoria Kepleriana.

¡La respuesta complementaria que confirma que @notovny es correcta!

Si bien vis-viva te da la velocidad, ¡aparentemente todas las direcciones parecen ser posibles!

Parece que me he desconcertado esta vez.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

def deriv(X, t):
    x, v = X.reshape(2, -1)
    acc = -x * ((x**2).sum())**-1.5
    return np.hstack((v, acc))

halfpi, pi, twopi   = [f*np.pi for f in (0.5, 1, 2)]
r       = 1.0
answerz = []
titles  = []

for r_over_a in (0.7, 1.4):
    titles.append('r/a = ' + str(round(r_over_a, 2)))
    answers = []
    a     = r / r_over_a
    T     = twopi * np.sqrt(a**3)
    times = np.linspace(0, T, 1001)
    v0    = np.sqrt(2./r - 1./a)

    thetas = np.linspace(0, pi, 8)[:-1] # make the result odd to avoid singularity

    for theta in thetas:
        s, c = [f(theta) for f in (np.sin, np.cos)]
        X0   = np.array([r, 0, s*v0, c*v0])
        answer, info = ODEint(deriv, X0, times, full_output=True)
        answers.append(answer)
    answerz.append(answers)

if True:
    fig = plt.figure()
    for i, (title, answers) in enumerate(zip(titles, answerz)):
        ax  = fig.add_subplot(2, 1, i+1)
        for a in answers:
            x, y = a.T[:2]
            ax.plot(x, y)
        ax.plot([0], [0], 'oy', markersize=12)
        ax.set_aspect('equal')
        ax.set_title(title, fontsize=16)
    plt.show()