Primerjava analitičnega in numeričnega rezultata - gravitacijska frača okoli Venere

Sedaj ko imamo simulator lahko preverimo natančnost analitičnega rezultata naloge iz vaj predmeta Naša in druga Osončja - gravitacijska frača okoli Venere.

Navodilo naloge¶

Satelit želimo poslati proti Veneri. Svoje potovanje začne v Zemljini orbiti ter po eliptičnem tiru potuje do Venerine orbite. Srečanje se zgodi pri pravi anomaliji $\Theta = -50^\circ$ . Satelit nato po hiperbolični orbiti potuje mimo Venere (najbližja točka je $h = 300\ \text{km}$ nad Venero) in se nato utiri v novo eliptično orbito okoli Sonca. Zanimata nas novi orbiti v primeru ko satelit potuje mimo Venere po osvetljeni ali temni strani. Poznamo naslednje vrednosti:

\begin{aligned} GM_\odot &= 1.3271\cdot 10^{11}\frac{\text{km}^3}{\text{s}^2},\\ a_Z &= 149.6\cdot 10^6\ \text{km},\\ a_V &= 108.2\cdot 10^6\ \text{km},\\ R_V &= 6050\ \text{km},\\ M_V &= 4.9\cdot 10^{24}\ \text{kg}. \end{aligned}

Matematično ozadje¶

V splošnem za poljubno stožnico lahko izračunamo naslednje parametre, ki popolnoma določijo obliko:

specifična orbitalna energija:
$\varepsilon = \frac{v^2}{2} - \frac{\mu}{r};\ \mu = G(m_1 + m_2),$
specifična vrtilna količina:
$\textbf{h} = \textbf{r}\times\textbf{v} = \frac{\textbf{L}}{m},$
ekscentrični vektor (ta kaže proti apsidni točki):
$\textbf{e} = \frac{\textbf{v}\times\textbf{h}}{\mu} - \frac{\textbf{r}}{r}\to e = |\textbf{e}|,$
ortonormirani sistem:
$\hat{\textbf{p}} = \frac{\textbf{e}}{e},\ \hat{\textbf{h}} = \frac{\textbf{h}}{h},\ \hat{\textbf{q}} = \hat{\textbf{h}}\times\hat{\textbf{p}}.$

Pozicijo poljubnega objekta tako lahko opišemo kot:

\textbf{r}(\nu) = \frac{p}{1 + e\cos\nu}(\cos\nu\hat{\textbf{p}} + \sin\nu\hat{\textbf{q}}).

Ideja rešitve¶

To nalogo smo na vajah rešili s kar nekaj poenostavitvami, česar si tu ne moremo prislužiti - tam smo namreč rekli da se satelit zabije v Venero in njegov odmik smo na roke nastavili tako, da je priletel na pravo višino nad površino. Tu tega ne moremo storiti tako, ampak moramo poiskati pravilen zamik Venere, da satelit zares leti po hiperboli na pravi višini. Satelit bo torej letel po enaki orbiti kot pri nalogi, Venero pa bomo za vsak primer ustrezno zamaknili, da bo satelit mimo nje šel po hiperboli, na pravilni višini od površine.

Da nalogo sploh lahko rešimo analitično moramo reševanje ločiti na tri dele: potovanje satelita do Venere po elipsi, oblet Venere po hiperboli in potovanje satelita okoli Sonca po novi elipsi. To moramo narediti, saj le tako lahko analitično rešimo problem (v vsakem delu gre tako za problem dveh teles). Tu takoj nastopi vprašanje, kje ločiti te tri primere - se pravi kaj je pogoj, da vemo ali smo v prvem, drugem ali tretjem delu naloge. Zaradi enostavnosti to naredimo tako, da med posameznima paroma ločimo glede na gravitacijsko silo Venere in Sonca na satelit (se pravi ločimo tam, kjer sta ta enaka).

Ko satelit prileti do Venere ima neko relativno hitrost, ki je odvisna od zamika Venere. To relativno hitrost moramo ujeti tako, da bo satelit šel po pravilni hiperboli okoli Venere - se pravi moramo rešiti sistem enačb, kjer iščemo zamik, da zadostimo pogojema.

Let do Venere¶

Začnimo preprosto. Satelit pošljemo po eliptični orbiti do Venere, kjer Veneri dopustimo poljubni začetni kot (to bo naša spremenljivka). Opis Venere je zelo preprost:

\textbf{r}_V = a_V\begin{bmatrix}\cos(\delta_V + \omega t) \\ \sin(\delta_V + \omega t)\end{bmatrix},\ \textbf{v}_V = \sqrt{\frac{GM}{a_V}}\begin{bmatrix}-\sin(\delta_V + \omega t) \\ \cos(\delta_V + \omega t)\end{bmatrix};\ \omega = \sqrt{\frac{GM}{a_V^3}}.

Opis satelita zahteva nekoliko več dela. Upoštevajmo enačbo elipse, kjer poznamo oddaljenosti pri dveh kotih:

\begin{aligned} r(\pi - \theta) &= \frac{a(1 - e^2)}{1 - e\cos(\pi - \theta)} = a_V,\\ r(0) &= \frac{a(1 - e^2)}{1 - e} = a(1 + e) = a_Z,\\ a_V &= \frac{a_Z}{1 + e}\frac{1 - e^2}{1 - e\cos(\pi - \theta)} = \frac{a_Z(1 - e)}{1 - e\cos(\pi - \theta)}. \end{aligned}

Takoj lahko izračunamo ekscentričnost naše orbite:

\begin{aligned} a_V(1 - e\cos(\pi - \theta)) &= a_Z(1 - e), \\ e(a_Z - a_V\cos(\pi - \theta)) &= a_Z - a_V,\\ e = \frac{a_Z - a_V}{a_Z - a_V\cos(\pi - \theta)}&\approx 0.1889,\\ a = \frac{a_Z}{1 + e} &\approx 125.8\cdot 10^6\ \text{km}. \end{aligned}

Sedaj lahko popolnoma opišemo pozicijo satelita, potrebujemo pa še vektor hitrosti. Posamezni komponenti hitrosti lahko izračunamo preko kota $\theta$ :

\begin{aligned} v_r &= \sqrt{\frac{G(M + m)}{a(1 - e^2)}}e\sin\theta,\\ v_\theta &= \sqrt{\frac{G(M + m)}{a(1 - e^2)}}(1 + e\cos\theta). \end{aligned}

Sedaj lahko popolnoma opišemo pozicijo satelita v orbiti s kotom $\theta$ , toda mi bi radi kot Venero, pozicijo opisali le s časom (tako z eno spremenljivko opišemo pozicijo obeh teles hkrati). Numerično moramo rešiti tako imenovano Keplerjevo enačbo za ekscentrično anomalijo in tako lahko izračunamo kot, posledično pa pozicijo in hitrost v odvisnosti od časa:

\begin{aligned} 2\pi\frac{t}{T} &= E - e\sin(E);\ T = 2\pi\sqrt{\frac{a^3}{G(M + m)}},\\ \tan\frac{\theta}{2} &= \sqrt{\frac{1 + e}{1 - e}}\tan\frac{E}{2}. \end{aligned}

Tako dobimo kot v odvisnosti od časa in lahko uporabimo vse enačbe, ki smo jih omenili do sedaj, da opišemo pozicijo in hitrost satelita. Kot vidimo že samo za reševanje enačbe kota v odvisnosti od časa potrebujemo numerično reševanje, torej od tu naprej ne gre le z enačbami. Ta korak končamo, ko je razdalja med satelitom in Venero taka, da je gravitacijski vpliv Venere in Sonca na satelit enak.

Oblet Venere¶

V trenutku, ko za satelit upoštevamo, da je v orbiti okoli Venere, moramo preiti iz sistema okoli Sonca, v sistem okoli Venere (hitrosti satelita odštejemo hitrost Venere). Satelitu nato lahko zelo preprosto izračunamo parametre, ki jih potrebujemo za hiperbolo. Tako lahko izračunam pozicijo satelita za poljubni kot na hiperboli. V začetni točki lahko preko skalarnega produkta zelo enostavno izračunam kot $\nu_0$ in posledično parameter $p$ , ter imam tako popolno opisani sistem. Postavim še zahtevo za rešitev:

r(0) = \frac{p}{1 + e} = R_V + h.

Rešitev bom tu seveda ponovno moral iskati numerično - spreminjati moram $\delta_V$ , dokler vrednost $r(0)$ ne zadosti zahtevi. Kot pri elipsi moram rešiti Keplerjevo enačbo (ki je sedaj v hiperbolični obliki), saj moram upoštevati premik Venere. Ko dosežem pogoj, da je vpliv Venere enak vplivu Sonca, izračunam hitrost satelita in preidem v Sončni sistem.

Orbita okoli Sonca¶

V zadnji točki obleta okoli Venere izračunam hitrost satelita in mu sedaj prištejem hitrost Venere (ponovno gremo v sistem Sonca). Podobno kot pri obletu Venere izračunam vse parametre in dobim tokrat dobim eliptično orbito okoli Sonca.

Omejitve rešitve¶

Kot omenjeno, smo na vajah nalogo zelo poenostavili, da smo jo lahko rešili popolnoma analitično. Kot vidimo, pri pravi “analitični” rešitvi potrebujemo več numeričnih rešitev, zato rešitev ni popolnoma pravilna (satelit ne gre točno $300\ \text{km}$ nad površino). Kljub temu je rešitev popolnoma določena, le vrednosti niso čisto take, kot smo jih zahtevali.

“Analitična” rešitev¶

Prvo moramo poiskati začetni $\delta_V$ , da bo satelit letel na pravi višini nad površino Venere. To najlažje naredimo s tako imenovano binarno metodo, kjer korake skaliramo glede na oddaljenost od željene vrednosti.

import numpy as np
import matplotlib.pyplot as plt
from simulatorv2 import Vec2

m_s = 500 # Masa satelita [kg]
m_S = 2e30 # Masa Sonca [kg]
m_V = 4.867e24 # Masa Venere [kg]
G = 6.6743e-11 # Gravitacijska konstanta [m^3kg^-1s^-2]
R_V = 6_050_000 # Radij Venere [m]
h_V = 300_000 # Višina na kateri naj satelit obleti Venero [m]
a_Z = 149.6e9 # Radij Zemljine orbite [m]
a_V = 108.2e9 # Radij Venerine orbite [m]


def solve_kepler_elliptic(M, e, tol=1e-12):
    E = np.array(M, dtype=float)
    
    while True:
        f = E - e * np.sin(E) - M
        df = 1 - e * np.cos(E)

        E_new = E - f / df
        if np.all(np.abs(E_new - E) < tol):
            return E_new

        E = E_new

def solve_kepler_hyperbolic(M, e, tol=1e-12):
    E = np.array(M, dtype=float)
    
    while True:
        f = E - e * np.sinh(E) - M
        df = 1 - e * np.cosh(E)

        E_new = E - f / df
        if np.all(np.abs(E_new - E) < tol):
            return E_new

        E = E_new

def true_anomaly_from_time_elliptic(t, T, e):
    M = 2 * np.pi * t / T
    E = solve_kepler_elliptic(M, e)
    theta = 2 * np.arctan2(
        np.sqrt(1 + e) * np.sin(E / 2),
        np.sqrt(1 - e) * np.cos(E / 2)
    )

    return theta

def conical(p, e, nu, p_vec=Vec2(1, 0), q_vec=Vec2(0, 1)):
    r = p / (1 - e * np.cos(nu))

    return (
        r * np.cos(nu) * p_vec +
        r * np.sin(nu) * q_vec
    )

def gravitational_force(m, M, R):
    return m * M * G / (R * R)

def transfer1(delta_V, tot=1e-12):
    a = 125.8e9
    T = 2 * np.pi * np.sqrt(np.power(a, 3) / (G * (m_s + m_S)))
    e = 0.1889
    r_p = a * (1 - e)

    t_steps = 100
    t = np.linspace(0, T/2, t_steps)
    
    nu_s = true_anomaly_from_time_elliptic(t, T, e)
    r_s = conical(a_Z * (1 - e), e, nu_s)

    nu_V = delta_V + true_anomaly_from_time_elliptic(t, T, 0)
    r_V = conical(a_V, 0, nu_V)

    for r_s_, r_V_ in zip(r_s, r_V):
        F_1 = gravitational_force(m_s, m_V, (r_s_ - r_V_).mag())
        F_2 = gravitational_force(m_s, m_S, r_s_.mag())

        if F_2 < F_1:
            return True
    return False
    
    v_s, v_V = [], []

    for r in r_V:
        n = r.normalize()
        v = np.sqrt(G * m_S / a_V)
        
        v_V.append(v * n.rotate(np.pi * 0.5))

    for nu in nu_s:
        f = np.sqrt(G * (m_S + m_s) / (a * (1 - e * e)))

        v_r = f * e * np.sin(nu)
        v_t = f * (1 + np.cos(nu))

        n = r.normalize()
        t = n.rotate(np.pi * 0.5)

        v_s.append(v_r * n + v_t * t)

delta_V = 0.32

0.3200000000000001