Vaja iz NDO - Sim-Astro

Sedaj ko imamo simulator simulirajmo rezultat naloge iz vaj predmeta Naša in druga Osončja - gravitacijska frača okoli Venere. 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 $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,\\ R_V &= 6050\ \text{km},\\ M_V &= 4.9\cdot 10^{24}\ \text{kg}. \end{aligned}

Nalogo moramo očitno rešiti v treh delih: prvo moramo izračunati vse parametre elipse do Venere, nato parametre hiperbole pri obletu Venere in na koncu parametre elipse po obletu.

Polet do Venere¶

Prvo satelit pošljemo po eliptični orbiti do Venere. Spomnimo se, da je prava anomalija definirana kot kot med smerjo proti periapsidi (najbližji točki gorišču v kateremu je Sonce) in smerjo proti trenutni legi telesa. Takoj lahko izračunamo ekscentričnost, saj namreč poznamo enačbo elipse iz izhodišča:

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

Iz prve enačbe izpostavimo $a$ in jo nesemo v drugo enačbo, ter tako dobimo:

\begin{aligned} a_V &= \frac{a_Z(1 + e)(1 - e)}{(1 + e)(1 + e\cos\Theta)} = \frac{a_Z(1 - e)}{1 + e\cos\Theta},\\ a_V + a_Ve\cos\Theta &= a_Z - a_Ze\to e = \frac{a_Z - a_V}{a_Z + a_V\cos\Theta} = 0.1889. \end{aligned}

Takoj lahko izračunamo še polosi elipse:

\begin{aligned} a &= \frac{a_Z}{1 + e} = 125.8\cdot 10^6\ \text{km},\\ c &= ea = 23.76\ \text{km},\\ b &= \sqrt{a^2 - c^2} = 123.5\cdot 10^6\ \text{km}. \end{aligned}

%load_ext jupyter_tikz

%%tikz -r
\usetikzlibrary{calc}
\begin{tikzpicture}
    \clip (-10,-10) rectangle (10, 10);

    \def\scale{7.5}
    \def\aZ{1 * \scale}
    \def\aV{0.723 * \scale}
    \def\RZ{5pt}
    \def\RV{5pt}
    \def\RS{10pt}
    \def\theta{130}

    \coordinate (A) at (\theta: \aV);
    \coordinate (B) at (0: \aZ);

    \def\cx{\aZ * 23.76 / 149.6}
    \def\rx{\aZ * 125.8 / 149.6}
    \def\ry{\aZ * 123.5 / 149.6}
    \draw[dashed] ({\cx + \rx}, 0) arc[start angle=0, end angle=180, x radius=\rx, y radius = \ry];
    \draw[dashed] (0, 0) -- (-\rx + \cx, 0);
    \draw[dashed] (0, 0) -- (A);
    \draw[dashed] (-30pt, 0) arc (180:\theta:30pt);
    \node[above left] at (90 - \theta/2 + \theta: 30pt) {$\Theta$};

    \draw[dashed, orange] (0, 0) circle (\aV);
    \draw[dashed, blue] (0, 0) circle (\aZ);
    \fill[orange] (A) circle (\RV);
    \fill[blue] (B) circle(\RZ);
    \node[above left=8pt] at (A) {V};
    \node[right=8pt] at (B) {Z};
    \fill[white] (0, 0) circle (\RS);
    \draw[thick] (0, 0) circle (\RS);
    \fill (0, 0) circle (1.5pt);
\end{tikzpicture}

Ko dosežemo Venero, lahko s parametri elipse izračunamo še hitrost satelita:

v_0 = \sqrt{GM_\odot\left(\frac{2}{a_V} - \frac{1}{a}\right)} = 37.39\ \frac{\text{km}}{\text{s}}.

Frača okoli Venere¶

Sedaj smo prileteli do Venere in moramo izračunati parametre hiperbole. Prvo izračunajmo vektorja hitrosti Venere in satelita, kar naredimo s pomočjo kota $\gamma$ (spomnimo se, da je to kot med smerjo proti gorišču in smerjo gibanja).

%%tikz -r
\usetikzlibrary{calc}
\begin{tikzpicture}
    \clip (-11,-1) rectangle (5, 7);

    \def\scale{7}
    \def\theta{140}
    \def\aV{0.723 * \scale}

    \coordinate (A) at (\theta:\aV);

    \draw[-stealth, dashed] (A) -- (0, 0);

    \def\rx{\aV*0.8}
    \def\ry{\aV*0.8}
    \draw[dashed] (A) + (\rx/4.5, -\ry/2) arc[start angle=220, end angle=\theta, x radius=\rx, y radius = \ry];
    \draw[-stealth, solid] (A) -- ++(.3, -2);
    \node[left=5pt] at ($(A)+(0.3, -1)$) {$\vec{v}_0$};
    \draw ($(A)+(1, -0.9)$) arc (-30:-110:20pt);
    \node at ($(A)+(0.75, -1.4)$) {$\gamma$};

    \draw[dashed, orange] (-\aV, 0) arc(180:\theta:\aV);
    \draw[-stealth, solid] (A) -- ++(-1, -1.3);
    \node at ($(A)+(-0.7, -0.4)$) {$\vec{v}_V$};
    \fill[orange] (A) circle (5pt);
    \node[above left=8pt] at (A) {V};

    \node[below right=8pt] at (0, 0) {S};
\end{tikzpicture}

Izračunamo lahko kot $\gamma$ in za prehod v sistem Venere izračunamo tangencialno in radialno hitrost satelita:

\begin{aligned} \sin\gamma &= \frac{\sqrt{\frac{2}{r_P} - \frac{1}{a}}}{\frac{2}{r} - \frac{1}{a}}\frac{1 + \varepsilon\cos\Theta}{1 + \varepsilon}\to\gamma = 82.72^\circ,\\ \vec{v}_{0,t} &= v_0\sin\gamma = 37.08\frac{\text{km}}{\text{s}},\\ \vec{v}_{0,r} &= v_0\cos\gamma = 4.74\frac{\text{km}}{\text{s}}. \end{aligned}

Izračunamo še hitrost Venere in ker ta kroži okoli Sonca, je njena hitrost kar tangencialna, zato lahko zelo preprosto preidemo v sistem Venere:

\begin{aligned} v_V &= \sqrt{\frac{GM\odot}{a_V}} = 35.02\frac{\text{km}}{\text{s}},\\ v_{1,t} &= (37.08 - 35.02)\frac{\text{km}}{\text{s}} = 2.06\frac{\text{km}}{\text{s}},\\ v_{1,r} &= 4.74\frac{\text{km}}{\text{s}}. \end{aligned}

Začnemo z računanjem parametrov hiperbole. Izračunamo lahko hitrost v neskončnosti in poznamo najbližjo razdaljo Veneri, zanima pa nas sprememba kota $\Theta_D$ :

\begin{aligned} |v_\infty| &= \sqrt{v_{1,t}^2 + v_{1,r}^2} = 5.16\frac{\text{km}}{\text{s}},\ r_h = a_h(e_h - 1) = R_V + 300\ \text{km},\\ v_\infty^2 &= \frac{GM_V}{a_h}\to e_h = \frac{r_hv_\infty^2}{GM_V} + 1 = 1.51\to a_h = 12450\ \text{km}. \end{aligned}

Spremembo kota izračunamo s pomočjo dveh enačb za parameter $b$ :

b^2 = a^2(e_h^2 - 1),\ b = a\cot\frac{\Theta_D}{2}\to\Theta_D = 2\tan\sqrt{\frac{1}{e_h^2 - 1}} = 82.9^\circ.

V obeh primerih začnemo z enakim kotom, hitrost pa se v neskončnosti ne spremeni. Kota izračunamo tako, da prištejemo ali odštejemo spremembo $\Theta_D$ :

\Theta_1 = \arctan\frac{v_{1,t}}{v_{1,r}} = 23.5^\circ\to\Theta_2 = \Theta_1 + \Theta_D = 106.4^\circ,\ \Theta_3 = \Theta_1 - \Theta_D = -59.4^\circ.

Izračunamo lahko še novi hitrosti v sončevem sistemu in primerjamo rezultata:

\begin{aligned} v_{2,r}&=|v_\infty|\cos\Theta_2=-1.45\frac{\text{km}}{\text{s}},\ v_{2,t}=|v_\infty|\sin\Theta_2=4.95\frac{\text{km}}{\text{s}}\to|\vec{v}_V + \vec{v}_2| = 39.99\frac{\text{km}}{\text{s}},\\ v_{3,r}&=|v_\infty|\cos\Theta_3=2.6\frac{\text{km}}{\text{s}},\ v_{3,t}=|v_\infty|\sin\Theta_3=-4.45\frac{\text{km}}{\text{s}}\to|\vec{v}_V + \vec{v}_3| = 30.68\frac{\text{km}}{\text{s}}. \end{aligned}

V primeru ko Venero satelit obleti po svetli strani gre torej za pospešek, ko po temni pa pojemek.

Simulacija¶

Da si malo bolj predstavljamo rezultat in da vidimo, kako se realnost razlikuje od idealističnih enačb, simulirajmo taki gravitacijski frači.

from simulatorv2 import *
import numpy as np

transfer_reference = ReferenceConic(eccentricity=0.1889, apoapsis=1, rotation=np.pi, t_0=0, t_1=np.pi, opacity=0.1)
sat = Body(
    name="Satelit",
    mass=500 / Units.SUN_MASS,
    radius=5 / Units.ASTRONOMICAL_UNIT,
    color=(200,200,200),
    position=Vec2(1,0),
    velocity=transfer_reference.velocityAtAngle(np.pi)
)
"""
# UPORABIL DA SEM ITERATIVNO NAŠEL REZULTATA
target_altitude = 300_000 / Units.ASTRONOMICAL_UNIT
guess = -52.0
step = 0.01
best = None

for _ in range(30):
    venus, sun = initialize_binary(Presets.venus(Vec2), Presets.sun(Vec2), Vec2(0.723, 0).rotate(np.deg2rad(guess)), Vec2(0, 0))
    sat = sat_og.copy()
    sim = Simulator(
        bodies=[sun, venus, sat],
        timestep=0.0001,
        steps=5000
    )
    result = sim.run()

    r = result.positions[:, 2] - result.positions[:, 1]
    d = np.linalg.norm(r, axis=1)
    altitude = d.min() - result.body_radii[1]
    error = altitude - target_altitude

    print(_, guess, error * Units.ASTRONOMICAL_UNIT)

    if best is None or abs(error) < abs(best[1]):
        best = (guess, error)

    step_size = np.arctan(abs(error * Units.ASTRONOMICAL_UNIT) / 5_000_000) / (np.pi * 0.5)
    if error < 0:
        guess += step * step_size
    else:
        guess -= step * step_size

print(best)
"""
r = -52.108266015541815
hr = 56.995+180

venus, sun = initialize_binary(Presets.venus(Vec2), Presets.sun(Vec2), Vec2(0.723, 0).rotate(np.deg2rad(r)), Vec2(0, 0))
sim = Simulator(
    bodies=[sun, venus, sat],
    timestep=0.0001,
    steps=5000,
    diagnostics=[
        TotalEnergy()
    ]
)
result = sim.run()

viewer = SimulationViewer(
    result,
    ViewerConfig(
        camera_mode="follow",
        follow_body="Venera",
        trail_length=1000,
        references=[
            ReferenceConic(opacity=0.1),
            ReferenceConic(periapsis=0.723, opacity=0.1, n=10000),
            ReferenceConic(periapsis=venus.radius + 300_000 / Units.ASTRONOMICAL_UNIT, opacity=0.1, follow_body="Venera"),
            ReferenceConic(
                periapsis=venus.radius + 300_000 / Units.ASTRONOMICAL_UNIT,
                eccentricity=1.51,
                max_distance=10 * venus.radius,
                opacity=0.1,
                follow_body="Venera",
                rotation=np.deg2rad(hr)
            ),
            transfer_reference
        ],
        diagnostics=[
            DiagnosticView(
                "Celotna Energija",
                ["Satelit"]
            ),
            DiagnosticView(
                "Celotna Energija",
                ["Venera"]
            )
        ],
    ),
)
viewer.show()

Loading...

from simulatorv2 import *
import numpy as np

transfer_reference = ReferenceConic(eccentricity=0.1889, apoapsis=1, rotation=np.pi, t_0=0, t_1=np.pi, opacity=0.1)
sat = Body(
    name="Satelit",
    mass=500 / Units.SUN_MASS,
    radius=5 / Units.ASTRONOMICAL_UNIT,
    color=(200,200,200),
    position=Vec2(1,0),
    velocity=transfer_reference.velocityAtAngle(np.pi)
)

r = -52.09225981719751
hr = -24.412

venus, sun = initialize_binary(Presets.venus(Vec2), Presets.sun(Vec2), Vec2(0.723, 0).rotate(np.deg2rad(r)), Vec2(0, 0))
sim = Simulator(
    bodies=[sun, venus, sat],
    timestep=0.0001,
    steps=5000,
    diagnostics=[
        TotalEnergy()
    ]
)
result = sim.run()

viewer = SimulationViewer(
    result,
    ViewerConfig(
        camera_mode="follow",
        follow_body="Venera",
        trail_length=1000,
        references=[
            ReferenceConic(opacity=0.1),
            ReferenceConic(periapsis=0.723, opacity=0.1, n=10000),
            ReferenceConic(periapsis=venus.radius + 300_000 / Units.ASTRONOMICAL_UNIT, opacity=0.1, follow_body="Venera"),
            ReferenceConic(
                periapsis=venus.radius + 300_000 / Units.ASTRONOMICAL_UNIT,
                eccentricity=1.51,
                max_distance=10 * venus.radius,
                opacity=0.1,
                follow_body="Venera",
                rotation=np.deg2rad(hr)
            ),
            transfer_reference
        ],
        diagnostics=[
            DiagnosticView(
                "Celotna Energija",
                ["Satelit"]
            ),
            DiagnosticView(
                "Celotna Energija",
                ["Venera"]
            )
        ],
    ),
)
viewer.show()

Loading...