@Julio „s doskonała odpowiedź opisuje kąt toru lotu i wyjaśnia, że jest kątem między kierunkiem stycznym (prostopadłym do wektora promieniowego do ciała centralnego) a aktualnym wektorem prędkości.
Najpierw próbowałem uzyskać kąt z tego wyrażenia, ale jest to oczywiście błędne , ponieważ $ \ arccos $ jest funkcją parzystą i kąt może wynosić od $ – \ pi / 2 $ do $ \ pi / 2 $:
$$ \ arccos \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) – \ frac {\ pi} {2} \ \ \ \ text {(niepoprawne!)} $$
Zintegrowałem orbity dla GM ($ \ mu $) i SMA ($ a $) jedności i odległości początkowych od 0,2 do 1,8. To sprawia, że okres zawsze wynosi 2 $ \ pi $. Kiedy wykreślam wynik mojej funkcji, dostaję zbyt wiele ruchów.
Jakiego wyrażenia mogę użyć, aby uzyskać poprawną gamma kąta toru lotu, zaczynając od wektorów stanu?
Poprawiony Python z powodu błędnej części byłby mile widziany, ale z pewnością nie jest konieczny do udzielenia odpowiedzi.
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() Komentarze
- jeśli to pytanie być TLDRed, aby wskazać, że był to błąd w kodowaniu, ponieważ nadal wydaje się pytać, co ' jest nie tak ze wzorem
Odpowiedź
Jest to problem, który nękał grupy ludzi posiadających dużą wiedzę na temat dynamiki orbitalnej, ale którzy nauczyli się korzystać z różnych podręczników: istnieją dwie różne definicje „kąta toru lotu „!!
Oprócz $ \ gamma $, kąta pomiędzy kierunkiem stycznym a wektorem prędkości, jest $ \ beta $, kąt pomiędzy kierunek promieniowy i wektor prędkości. Ludzie często mówią „kąt toru lotu”, nie mówiąc jakiej definicji „używają . Mylące! (Właśnie zauważyłem, że diagram w odpowiedzi Julio również pokazuje $ \ beta $)
Jeśli pracujesz z $ \ beta $ zamiast $ \ gamma $, $ \ beta $ jest podane przez
$$ \ arccos \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {1} $$
, który przechodzi od 0 („prosto w górę”) do $ \ pi $ ( „prosto w dół”). Używając $ \ gamma $, „prosto w górę” to $ \ pi / 2 $, a „prosto w dół” to $ – \ pi / 2 $, więc przeliczając $ \ beta $ na $ \ gamma $ wystarczy odjąć $ \ beta $ od $ \ pi / 2 $:
$$ \ gamma = \ pi / 2 – \ arccos \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {2} $$
Jest to równoważne z
$$ \ gamma = \ arcsin \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {3} $$
Nie znam języka użyłeś do swoich obliczeń i wykresów, więc nie przyjrzałem się Twojemu algorytmowi, aby zobaczyć, dlaczego jest „zbyt wiele ruchów”.
Komentarze
- Dzięki! Dodałem ' tagi (liczby) do równań. Czy powiedziałbyś, że jest zbyt wiele ruchów, czy też takie zachowanie jest w rzeczywistości uzasadnione? Ponieważ twoje $ \ beta $ (eq 1) jest takie samo jak moje błędne $ \ gamma $ z wyjątkiem przesunięcia o połowę pi, to ruchy na moim wykresie powinny być takie same, jak we właściwym wykresie twojego $ \ beta $ (równ. 1).
- Wydaje mi się, że za dużo się kręci. ' sprawdzę to później.
- @uhoh, Właściwie to moje równanie 1 jest po prostu minusem twojego równania. Coś innego jest nie tak. Oczywiście wiesz, że coś jest nie tak, ponieważ wszystkie wykreślone $ \ gamma $ s są ujemne lub zerowe, co ' nie może być poza spiralą wewnętrzną. Dla keplerowskiej mimośrodowej orbity $ \ gamma $ powinno przecinać zero dokładnie dwa razy, w perycentrum i apocentrum, i być monotoniczne między ekstremami, zarówno dla krótkiej (ekstremum przez perycentrum do drugiego ekstremum), jak i długiej (ekstremum przez apocentrum do drugiego ekstremum ) segmenty. I ' zobaczę, czy mogę narysować przykład tego, jak powinna wyglądać krzywa $ \ gamma $.
- Ups, powinienem był powiedzieć powyżej, ” Mój eq. 2 to tylko twój minus. ” Powinienem się wylogować i iść do łóżka!
- @uhoh ” styczna ” do kuli wyśrodkowanej na stronie pierwotnej, a nie do orbity. Osobiście ' d wolę powiedzieć ” prędkość boczna „, ale mój pierwszy profesor dynamiki orbitalnej na Stanford użył ” stycznej „.
Odpowiedź
Znalazłem błąd w skrypcie, był on spowodowany moim iloczynem kropkowym „homebrew”.Miałem dodatkowe kwadraty:
dotted = ((xx*vv)**2).sum(axis=0) # WRONG dotted = (xx*vv).sum(axis=0) # Correct Więc używając tego plus @TomSpilker „s doskonałe wyjaśnienia Używam następujących dwóch metod do obliczenia współczynnika gamma:
Metoda 1:
$$ \ gamma_1 = \ arcsin \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {3} $$
Metoda 2:
Brute-force alternatywna metoda podwójnego sprawdzenia:
$$ \ 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 $$
Operacja modulo jest naprawdę potrzebna tylko w programie komputerowym, ponieważ każda theta pochodzi z oddzielnej operacji 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) Komentarze
- Rzeczywiście nowy $ \ działka gamma $ jest tym, czego się spodziewałem. Brawo! Dobre śledztwo.