@Julio „s vynikající odpověď popisuje úhel dráhy letu a vysvětluje, že je úhel mezi tangenciálním směrem (kolmý na radiální vektor k centrálnímu tělu) a aktuálním vektorem rychlosti.
Nejprve jsem se pokusil získat úhel z tohoto výrazu, ale je zjevně špatný , protože $ \ arccos $ je sudá funkce a úhel může jít od $ – \ pi / 2 $ do $ \ pi / 2 $:
$$ \ arccos \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) – \ frac {\ pi} {2} \ \ \ \ text {(nesprávné!)} $$
Mám integrované oběžné dráhy pro jednotu GM ($ \ mu $) a SMA ($ a $) a počáteční vzdálenosti od 0,2 do 1,8. To činí období vždy $ 2 \ pi $. Když vykreslím výsledek své funkce, dostávám příliš mnoho kroutí.
Jaký výraz mohu použít k získání správného úhlu dráhy letu gama počínaje od stavových vektorů?
Revidovaný python pro chybnou část by byl oceněn, ale rozhodně není nutný pro odpověď.
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() Komentáře
- měla by tato otázka být TLDRed, což znamená, že se jednalo o chybu kódování, protože se stále zdá, že se ptáte, co ' je špatné se vzorcem
odpověď
Toto je problém, který sužuje skupiny lidí, kteří mají velmi dobré znalosti o orbitální dynamice, ale kteří se naučili používat různé učebnice: existují dvě různé definice „úhlu dráhy letu“ „!!
Kromě $ \ gamma $ je úhel mezi tangenciálním směrem a vektorem rychlosti také $ \ beta $, úhel mezi radiální směr a vektor rychlosti. Lidé často říkají „úhel dráhy letu“, aniž by říkali, jakou definici používají . Matoucí! (Jen jsem si všiml, že diagram v Juliově odpovědi také ukazuje $ \ beta $)
Pokud pracujete s $ \ beta $ místo $ \ gamma $, $ \ beta $ je dáno
$$ \ arccos \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {1} $$
od 0 („straight up“) do $ \ pi $ ( „rovnou dolů“). Při použití $ \ gamma $ je „straight up“ $ \ pi / 2 $ a „straight down“ je $ – \ pi / 2 $, takže při převodu $ \ beta $ na $ \ gamma $ stačí odečíst $ \ beta $ od $ \ pi / 2 $:
$$ \ gamma = \ pi / 2 – \ arccos \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {2} $$
Toto je ekvivalent
$$ \ gamma = \ arcsin \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {3} $$
Nejsem obeznámen s jazykem použili jste pro své výpočty a grafy, takže jsem se na váš algoritmus nepodíval, abych zjistil, proč je „příliš mnoho kroutí“.
Komentáře
- Díky! ' Přidal jsem značky (čísla) do rovnic. Řekli byste, že je příliš mnoho kroutí, nebo je to kroutící se chování ve skutečnosti rozumné? Vzhledem k tomu, že váš $ \ beta $ (ekv. 1) je stejný jako můj chybný $ \ gamma $ s výjimkou posunutí o polovinu pí, pak by měly být krouticí momenty v mém grafu stejné jako u správného grafu vašeho $ \ beta $ (ekv. 1).
- Vypadá to, že se mi příliš třese. ' to zkontroluji později.
- @uhoh, můj eq 1 je vlastně zápor vaší rovnice. Něco jiného není v pořádku. Samozřejmě víte, že něco je špatně, protože všechny vykreslené $ \ gamma $ s jsou záporné nebo nulové, což ' nemůže být výjimkou vnitřní spirály. Pro kepleriánskou excentrickou dráhu $ \ gamma $ by měla překročit nulu přesně dvakrát, při periapsi a apoapsi, a měla by být monotónní mezi extrémy, a to jak krátkou (extremum přes periapsis do druhého extrému), tak dlouhou (extremum přes apoapsis do druhého extrému) ) segmenty. ' Uvidím, jestli mohu nakreslit příklad toho, jak by měla vypadat křivka $ \ gamma $.
- Jejda, měl jsem říci výše, " Můj ekv. 2 je jen váš zápor. " Měla bych se odhlásit a jít spát!
- @uhoh " tangenciální " ke kouli soustředěné na primární straně, nikoli na oběžnou dráhu. Osobně ' d raději říkám " boční rychlost ", ale můj první profesor orbitální dynamiky ve Stanfordu používá " tangenciální ".
odpověď
Zjistil jsem chybu ve skriptu, byla to kvůli mému produktu „homebrew“ dot.Měl jsem extra kvadraturu:
dotted = ((xx*vv)**2).sum(axis=0) # WRONG dotted = (xx*vv).sum(axis=0) # Correct Takže pomocí tohoto plus @TomSpilker „s vynikající vysvětlení K výpočtu gama používám následující dvě metody:
Metoda 1:
$$ \ gamma_1 = \ arcsin \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {3} $$
Metoda 2:
Alternativní metoda brutální síly pro dvojitou kontrolu:
$$ \ 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 $$
Operace modulo je skutečně nutná pouze v počítačovém programu, protože každá theta pochází ze samostatné operace 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) Komentáře
- Opravdu nový $ \ gama $ plot je to, co jsem očekával. Hurá! Dobré odpalování.