@Julio “s resposta excelente descreve o ângulo da trajetória de voo e explica que é o ângulo entre a direção tangencial (perpendicular ao vetor radial ao corpo central) e o vetor de velocidade atual.
Primeiro tentei obter o ângulo a partir desta expressão, mas está obviamente errado , uma vez que $ \ arccos $ é uma função par e o ângulo pode ir de $ – \ pi / 2 $ a $ \ pi / 2 $:
$$ \ arccos \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) – \ frac {\ pi} {2} \ \ \ \ text {(incorreto!)} $$
Eu integrei órbitas para GM ($ \ mu $) e SMA ($ a $) de unidade e distâncias iniciais de 0,2 a 1,8. Isso torna o período sempre $ 2 \ pi $. Quando ploto o resultado de minha função, recebo muitos movimentos.
Que expressão posso usar para obter a gama correta do ângulo da trajetória de vôo a partir de vetores de estado?
Python revisado para a parte errada seria apreciado, mas certamente não é necessário para uma resposta.
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()
Comentários
- deve esta pergunta ser TLDR para indicar que foi um erro de codificação, pois ainda parece perguntar o que ‘ está errado com a fórmula
Resposta
Este é um problema que atormentou grupos de pessoas com muito conhecimento sobre dinâmica orbital, mas que aprenderam usando diferentes livros: existem duas definições diferentes de “ângulo da trajetória de vôo “!!
Além de $ \ gamma $, o ângulo entre a direção tangencial e o vetor de velocidade, existe $ \ beta $, o ângulo entre a direção radial e o vetor de velocidade. As pessoas costumam dizer “ângulo da trajetória de voo” sem dizer qual definição estão usando . Confuso! (Acabei de notar que o diagrama na resposta de Julio também mostra $ \ beta $)
Se você trabalhar com $ \ beta $ em vez de $ \ gamma $, $ \ beta $ é dado por
$$ \ arccos \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {1} $$
que vai de 0 (“direto para cima”) a $ \ pi $ ( “diretamente para baixo”). Usando $ \ gamma $, “direto para cima” é $ \ pi / 2 $ e “direto para baixo” é $ – \ pi / 2 $, então convertendo $ \ beta $ em $ \ gamma $ você apenas subtrai $ \ beta $ de $ \ pi / 2 $:
$$ \ gamma = \ pi / 2 – \ arccos \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {2} $$
Isso é equivalente a
$$ \ gamma = \ arcsin \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {3} $$
Não estou familiarizado com o idioma você usou para seus cálculos e gráficos, então não olhei seu algoritmo para ver por que há “muitos wiggles”.
Comentários
- Obrigado! Eu ‘ adicionei tags (números) às equações. Você diria que há muitos wiggles, ou esse comportamento de wiggling é de fato razoável? Uma vez que seu $ \ beta $ (eq 1) é o mesmo que meu $ \ gamma $ errôneo, exceto por um deslocamento de meio pi, então os meneios em meu gráfico devem ser os mesmos que aqueles em um gráfico adequado de seu $ \ beta $ (eq 1).
- Parece que muitos meneios para mim. Eu ‘ verificarei isso mais tarde.
- @uhoh, Na verdade, minha eq 1 é apenas o negativo de sua equação. Outra coisa está errada. Claro que você sabe que algo está errado porque todos os $ \ gamma $ s plotados são negativos ou zero, o que pode ‘ não ser, exceto por uma espiral interna. Para uma órbita excêntrica Kepleriana $ \ gamma $ deve cruzar zero exatamente duas vezes, no periapsis e apoapsis, e ser monotônico entre os extremos, tanto para o curto (extremo através do periapsis para o outro extremo) e longo (extremo através do apoapsis para o outro extremo) ) segmentos. Eu ‘ verei se consigo desenhar um exemplo de como a curva $ \ gamma $ deve se parecer.
- Ops, eu deveria ter dito acima, ” Minha eq. 2 é apenas o seu negativo. ” Devo fazer logoff e ir para a cama!
- @uhoh ” tangencial ” a uma esfera centrada no primário, não na órbita. Pessoalmente, eu ‘ d prefiro dizer ” velocidade lateral “, mas meu primeiro Orbital Dynamics em Stanford usou ” tangencial “.
Resposta
Encontrei o erro no script, era devido ao meu produto escalar “homebrew”.Tive um quadrado extra:
dotted = ((xx*vv)**2).sum(axis=0) # WRONG dotted = (xx*vv).sum(axis=0) # Correct
Portanto, usando este mais @TomSpilker “s excelentes esclarecimentos Eu uso os dois métodos a seguir para calcular gama:
Método 1:
$$ \ gamma_1 = \ arcsin \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {3} $$
Método 2:
Um método alternativo de força bruta para verificar novamente:
$$ \ theta_r = \ arctan2 (y, x) $$
$$ \ theta_v = \ arctan2 (vy, x) $$
$$ \ theta_ {tanj} = \ theta_r + \ frac {\ pi} {2} $$
$$ \ gamma_2 = \ theta_ {tanj} – \ theta_v $$
$$ \ gamma_ {2mod} = \ mod (\ gamma_2 + \ pi, 2 \ pi) – \ pi $$
A operação do módulo só é realmente necessária no programa de computador, pois cada teta vem de uma operação separada do arctan2:
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)
Comentários
- Certamente o novo $ \ gamma $ plot é o que eu esperava. Hooray! Boa investigação.