@Julio „s răspuns excelent descrie unghiul traiectului de zbor și explică faptul că este unghiul dintre direcția tangențială (perpendiculară pe vectorul radial față de corpul central) și vectorul vitezei curente.
Am încercat mai întâi să obțin unghiul din această expresie, dar este în mod evident greșit , deoarece $ \ arccos $ este o funcție uniformă și unghiul poate merge de la $ – \ pi / 2 $ la $ \ pi / 2 $:
$$ \ arccos \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) – \ frac {\ pi} {2} \ \ \ \ text {(incorect!)} $$
Am integrat orbite pentru GM ($ \ mu $) și SMA ($ a $) de unitate și distanțe de pornire de la 0,2 la 1,8. Aceasta face ca perioada să fie întotdeauna de 2 $ \ pi $. Când trasez rezultatul funcției mele, primesc prea multe clătinări.
Ce expresie pot folosi pentru a obține unghiul corect al traiectoriei de zbor pornind de la vectorii de stare?
Python revizuit pentru partea eronată ar fi apreciat, dar cu siguranță nu este necesar pentru un răspuns.
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() Comentarii
- în cazul în care această întrebare fiți TLDRed pentru a indica faptul că a fost o eroare de codare, deoarece încă pare să întrebați ce ‘ este greșit cu formula
Răspuns
Aceasta este o problemă care a afectat grupuri de oameni foarte cunoscuți despre dinamica orbitală, dar care au învățat folosind manuale diferite: există două definiții diferite ale „unghiului traiectului de zbor” „!!
În plus față de $ \ gamma $, unghiul dintre direcția tangențială și vectorul viteză, există $ \ beta $, unghiul dintre direcția radială și vectorul viteză. Oamenii spun adesea „unghiul căii de zbor” fără să spună ce definiție „folosesc . Confuz! (Am observat doar că diagrama din răspunsul lui Julio arată și $ \ beta $)
Dacă lucrați cu $ \ beta $ în loc de $ \ gamma $, $ \ beta $ este dat de
$$ \ arccos \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {1} $$
care merge de la 0 („direct în sus”) la $ \ pi $ ( „drept în jos”). Folosind $ \ gamma $, „drept în sus” este $ \ pi / 2 $ și „drept în jos” este $ – \ pi / 2 $, deci convertind $ \ beta $ în $ \ gamma $ trebuie doar să scădeți $ \ beta $ din $ \ pi / 2 $:
$$ \ gamma = \ pi / 2 – \ arccos \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {2} $$
Aceasta este echivalentă cu
$$ \ gamma = \ arcsin \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {3} $$
Nu cunosc limba pe care l-ați folosit pentru calcule și grafice, așa că nu m-am uitat la algoritmul dvs. pentru a vedea de ce există „prea multe mișcări”.
Comentarii
- Mulțumesc! Am ‘ am adăugat etichete (numere) la ecuații. Ați spune că există prea multe clătinări sau este de fapt rezonabil acest comportament de mișcare? Deoarece $ \ beta $ (eq 1) este același cu $ \ gamma $ eronat al meu, cu excepția unui offset de jumătate pi, atunci mișcările din complotul meu ar trebui să fie aceleași ca cele dintr-un complot corespunzător al $ \ beta $ dvs. (ec. 1).
- Mi se pare prea multe mișcări. ‘ voi verifica asta mai târziu.
- @uhoh, De fapt, echivalentul meu 1 este doar negativ al ecuației tale. Altceva nu este în regulă. Bineînțeles că știi că ceva este greșit, deoarece toate $ \ gamma $ s reprezentate sunt negative sau zero, ceea ce nu poate fi ‘ cu excepția unei spirale interioare. Pentru o orbită excentrică Kepleriană $ \ gamma $ ar trebui să traverseze zero exact de două ori, la periapsis și apoapsis, și să fie monoton între extrema, atât pentru scurt (extrem prin periapsis către celălalt extrem), cât și pentru lung (extrem prin apoapsis către celălalt extrem ) segmente. ‘ voi vedea dacă pot desena un exemplu de cum ar trebui să arate curba $ \ gamma $.
- Hopa, ar fi trebuit să spun mai sus, ” Eq. 2 este doar negativul dvs. ” Ar trebui să mă deconectez și să mă culc!
- @uhoh ” tangențială ” la o sferă centrată la primar, nu la orbită. Personal, ‘ prefer să spun ” viteza laterală „, dar primul meu profesor de dinamică orbitală la Stanford a folosit ” tangențială „.
Răspuns
Am găsit eroarea în script, se datora produsului meu dot „homebrew”.Am avut un cadru suplimentar:
dotted = ((xx*vv)**2).sum(axis=0) # WRONG dotted = (xx*vv).sum(axis=0) # Correct Deci, folosind acest plus @TomSpilker „s clarificări excelente Am folosit următoarele două metode pentru a calcula gama:
Metoda 1:
$$ \ gamma_1 = \ arcsin \ left (\ frac {\ mathbf {r \ centerdot v}} {| \ mathbf {r} | \ | \ mathbf {v} |} \ right) \ tag {3} $$
Metoda 2:
O metodă alternativă de forță brută pentru verificarea dublă:
$$ \ 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 $$
Operația modulo este necesară doar în programul de computer, deoarece fiecare theta provine dintr-o operațiune arctan2 separată:
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) Comentarii
- Într-adevăr, noul $ \ complot gamma $ este ceea ce mă așteptam. Ura! Sănătate bună.