La excelente respuesta de @Julio describe un ángulo de trayectoria de vuelo y explica que es el ángulo entre la dirección tangencial (perpendicular al vector radial al cuerpo central) y el vector de velocidad actual.
Primero intenté obtener el ángulo de esta expresión, pero obviamente es incorrecto, ya que es una función par y el ángulo puede ir de a :
He integrado órbitas para GM ( ) y SMA ( ) de la unidad y distancias iniciales de 0,2 a 1,8. Eso hace que el período siempre . Cuando trazo el resultado de mi función, obtengo demasiados cambios.
¿Qué expresión puedo usar para obtener el ángulo gamma correcto de la ruta de vuelo a partir de los vectores de estado?
Se agradecería Python revisado por la parte errónea, pero ciertamente no es necesario para una respuesta.
def deriv(X, t):
x, v = X.reshape(2, -1)
acc = -x * ((x**2).sum())**-1.5
return np.hstack((v, acc))
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint
halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
T = twopi
time = np.linspace(0, twopi, 201)
a = 1.0
rstarts = 0.2 * np.arange(1, 10)
vstarts = np.sqrt(2./rstarts - 1./a) # from vis-viva equation
answers = []
for r, v in zip(rstarts, vstarts):
X0 = np.array([r, 0, 0, v])
answer, info = ODEint(deriv, X0, time, full_output= True)
answers.append(answer.T)
gammas = []
for a in answers:
xx, vv = a.reshape(2, 2, -1)
dotted = ((xx*vv)**2).sum(axis=0)
rabs, vabs = [np.sqrt((thing**2).sum(axis=0)) for thing in (xx, vv)]
gamma = np.arccos(dotted/(rabs*vabs)) - halfpi
gammas.append(gamma)
if True:
plt.figure()
plt.subplot(4, 1, 1)
for x, y, vx, vy in answers:
plt.plot(x, y)
plt.plot(x[:1], y[:1], '.k')
plt.plot([0], [0], 'ok')
plt.title('y vs x')
plt.subplot(4, 1, 2)
for x, y, vx, vy in answers:
plt.plot(time, x, '-b')
plt.plot(time, y, '--r')
plt.title('x (blue) y (red, dashed)')
plt.xlim(0, twopi)
plt.subplot(4, 1, 3)
for x, y, vx, vy in answers:
plt.plot(time, vx, '-b')
plt.plot(time, vy, '--r')
plt.title('vx (blue) vy (red), dashed')
plt.xlim(0, twopi)
plt.subplot(4, 1, 4)
for gamma in gammas:
plt.plot(time, gamma)
plt.title('gamma?')
plt.xlim(0, twopi)
plt.show()
Este es un problema que ha afectado a grupos de personas muy conocedoras de la dinámica orbital pero que aprendieron usando diferentes libros de texto: ¡hay dos definiciones diferentes de "ángulo de trayectoria de vuelo"!
Además de , el ángulo entre la dirección tangencial y el vector velocidad, hay , el ángulo entre la dirección radial y el vector velocidad. La gente a menudo dice "ángulo de trayectoria de vuelo" sin decir qué definición están usando . ¡Confuso! (Acabo de notar que el diagrama en la respuesta de Julio también muestra )
si trabajas con en lugar de , es dado por
que va de 0 ("hacia arriba") a ("hacia abajo"). Usando , "hacia arriba" es y "hacia abajo" es , por lo que convertir a solo restas de :
Esto es equivalente a
No estoy familiarizado con el lenguaje que usó para sus cálculos y gráficos, por lo que no he mirado su algoritmo para ver por qué hay "demasiadas oscilaciones".
Encontré el error en el script, se debió a mi producto de punto "homebrew". Tuve un cuadrado extra:
dotted = ((xx*vv)**2).sum(axis=0) # WRONG
dotted = (xx*vv).sum(axis=0) # Correct
Entonces, usando esto más las excelentes aclaraciones de @TomSpilker, he usado los siguientes dos métodos para calcular gamma:
Método 1:
Método 2:
Un método alternativo de fuerza bruta para verificar dos veces:
La operación de módulo solo es realmente necesaria en el programa de computadora, ya que cada theta proviene de una operación arctan2 separada:
gammas_1, gammas_2 = [], []
for a in answers:
xx, vv = a.reshape(2, 2, -1)
dotted = (xx*vv).sum(axis=0)
rabs, vabs = [np.sqrt((thing**2).sum(axis=0)) for thing in (xx, vv)]
gamma_1 = np.arcsin(dotted/(rabs*vabs)) # Per Tom Spilker's answer Eq. 3
theta_r = np.arctan2(xx[1], xx[0])
theta_v = np.arctan2(vv[1], vv[0])
theta_tanj = theta_r + halfpi
gamma_2 = theta_tanj - theta_v
gamma_2 = np.mod(gamma_2 + pi, twopi) - pi
gammas_1.append(gamma_1)
gammas_2.append(gamma_2)
plt.figure()
plt.subplot(2, 1, 1)
for gamma_1 in gammas_1:
plt.plot(time, gamma_1)
plt.title('gammas_1', fontsize=16)
plt.subplot(2, 1, 2)
for gamma_2 in gammas_2:
plt.plot(time, gamma_2)
plt.title('gammas_2', fontsize=16)
usuario20636