Problem dveh teles

Kot smo navajeni, opišimo sistem dveh teles. Gibanje nebesnih teles je prvi seveda opisal Newton z gravitacijsko silo, ki pravi, da je sila drugega na prvo telo enaka:

\vec{F}_{12} = G\frac{m_1m_2}{r_{12}^2}\frac{\vec{r}_{12}}{r_{12}};\ \vec{r}_{12}\equiv\vec{r}_2 - \vec{r}_1,

kjer je $G = 6,674\cdot 10^{-11}\text{m}^3\ \text{kg}^{-1}\text{s}^{-2}$ gravitacijska konstanta. Da popolnoma opišemo stanje takega sistema očitno potrebujemo 6 spremenljivk (vsako telo ima namreč tri proste koordinate). Ekvivalentno lahko vpeljemo dva nova vektorja - težišče in vektor razlike - ki opišeta enak sistem:

\vec{R} = \frac{m_1\vec{r}_1 + m_2\vec{r}_2}{m_1 + m_2},\ \vec{r} = \vec{r}_2 - \vec{r}_1.

Z obravnavamo v Lagrangeovem formalizmu, ki zahteva kar nekaj izpeljave, dobimo končni rezultat za obliko orbite v polarnih koordinatah kot:

\frac{1}{r} = C\left[1 + e\cos(\theta - \theta')\right];\ e = \sqrt{1 + \frac{2EL^2}{mk^2}}.

To je enačba stožnice, kjer opazimo, da je najbolj pomemben parameter $e$ , ki mu pravimo ekscentričnost. Opazimo, da je ta predvsem odvisen od energije in vrtilne količine sistema. To pomeni, da če imamo sistem z dvema telesoma in jima ne spreminjamo mase, je zelo pomembno kako se v določenih točkah ti dve telesi premikata (se pravi tudi če imamo dva primera, kjer sta hitrosti enaki, je zelo pomembna tudi orientacija vektorjev teh hitrosti). Da iz izhodišča dobimo oblike, ki gredo skozi isto točko, moramo radij pomnožiti za faktor $(1 + e)$ (efektivno zmanjšamo obliko, toda to naredimo uniformno).

import numpy as np
import matplotlib.pyplot as plt

theta = np.linspace(-np.pi + 0.01, np.pi - 0.01, 2000)
plt.figure(figsize=(6,6))

for e in [0, 0.5, 1.0, 1.5]:
    r = (1 + e) / (1 + e * np.cos(theta))
    
    mask = np.isfinite(r) & (r > 0)
    theta_valid = theta[mask]
    r_valid = r[mask]

    x = r_valid * np.cos(theta_valid)
    y = r_valid * np.sin(theta_valid)
    plt.plot(x, y, label=f"e = {e}")

plt.grid()
plt.legend()
plt.title("Stožnice z različnimi $e$")

plt.xlim(-4, 4)
plt.ylim(-4, 4)

plt.show()

Te oblike se imenujejo stožnice, ker jih dobimo s presekom ravnine skozi osnovni stožec, kjer nam ekscentričnost podaja tangens kota normale ravnine.

%matplotlib widget

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
from matplotlib.gridspec import GridSpec

u = np.linspace(0, 2*np.pi, 200)
v = np.linspace(0, 5, 200)
U, V = np.meshgrid(u, v)

Xcone = V * np.cos(U)
Ycone = V * np.sin(U)
Zcone = V

fig = plt.figure(figsize=(8, 10))
gs = GridSpec(2, 1, height_ratios=[3, 1], figure=fig)

ax = fig.add_subplot(gs[0], projection="3d")
ax2 = fig.add_subplot(gs[1])

plt.subplots_adjust(bottom=0.12)

ax.plot_surface(Xcone, Ycone, Zcone, alpha=0.25)

ax.set_xlim(-6, 6)
ax.set_ylim(-6, 6)
ax.set_zlim(0, 6)

ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title("Stožec in ravnina")

e_vals = np.linspace(0, 3, 300)
theta_vals = np.degrees(np.arctan(e_vals))

ax2.plot(e_vals, theta_vals, label=r"$\theta(e)=\arctan(e)$")
ax2.set_xlabel("e")
ax2.set_ylabel("θ [°]")
ax2.set_title("Kot normale ravnine")
ax2.grid(True)

dot, = ax2.plot([0], [0], "ro")

def intersection(e):

    theta = np.arctan(e)

    t = np.linspace(0, 2*np.pi, 2000)
    c = 2

    denom = 1 - np.sin(t) * np.tan(theta)

    denom[np.abs(denom) < 1e-3] = np.nan

    r = c / denom

    x = r * np.cos(t)
    y = r * np.sin(t)
    z = r

    mask = np.isfinite(r)

    x, y, z = x[mask], y[mask], z[mask]

    mask2 = z >= 0

    return x[mask2], y[mask2], z[mask2]

e0 = 0
xi, yi, zi = intersection(e0)

line, = ax.plot(xi, yi, zi, linewidth=2, color="black")

slider_ax = fig.add_axes([0.2, 0.01, 0.6, 0.04])

slider = Slider(
    slider_ax,
    "e",
    0.0,
    3.0,
    valinit=e0,
    valstep=0.05
)

def update(val):

    global line

    e = slider.val

    line.remove()
    xi, yi, zi = intersection(e)
    line, = ax.plot(xi, yi, zi, linewidth=2, color="black")

    dot.set_data([e], [np.degrees(np.arctan(e))])

    fig.canvas.draw_idle()

slider.on_changed(update)

plt.show()

Kot vidimo, se orbite z večanjem ekscentričnosti oddaljujejo od krožnice. Zaradi enostavnosti najlažje vrtilno količino orbite izračunamo v najbližji točki - z drugim imenom ji rečemo tudi teme ali vrh - saj sta tam radialni vektor in hitrostni vektor vedno pravokotna. Če pogledamo enačbo za stožnico opazimo, da je ta zelo enostavne oblike v dveh primerih: ko je $e=0$ (takrat ji rečemo krožnica) in ko je $e=1$ (takrat ji rečemo parabola). V splošnem vezanim/zaključenim oblikam pravimo elipse, prostim/nezaključenim pa hiperbole.

Testna simulacija¶

Za problem dveh teles lahko ponovno napovemo kako naj bi se gibali telesi zato to uporabimo za še en test simulatorja - pošljimo Zemljo okoli Sonca po trajektorijah na zgornji sliki. Za katerokoli orbito velja formula:

v^2 = GM\left(\frac{2}{r} - \frac{1}{a}\right)\to v^2 = 4\pi^2\left(2 - \frac{1}{a}\right);\ r_p = a(1 - e) = 1\to a = \frac{1}{1 - e}.

Sedaj lahko zelo enostavno dobimo hitrost za vsak primer in to lahko prikažemo v simulatorju.

from simulatorv2 import *
from numpy import sqrt

results = []

G = NewtonianGravity().G

b1 = Presets.sun(Vec2)
b2 = Presets.earth(Vec2)
b2.position = Vec2(1.0, 0)

for e in [0, 0.5, 1.0, 1.5]:
    v = sqrt(G * (1 + e))
    b2.velocity = Vec2(0, v)

    sim = Simulator(
        bodies=[b1, b2],
        timestep=0.001,
        steps=1000
    )

    results.append(sim.run())

viewer = SimulationViewer(
    results,
    ViewerConfig(
        fixed_center=(0, 0),
        trail_length=10000,
        view_size=4,
        references=[
            ReferenceConic(0.0, opacity=0.5),
            ReferenceConic(0.5, opacity=0.5),
            ReferenceConic(1.0, opacity=0.5),
            ReferenceConic(1.5, opacity=0.5),
        ]
    )
)

viewer.show()

sun = Presets.sun(Vec3)

earth = Body(
    name="Earth",
    mass=3.003e-6,
    radius=6_371_000 / Units.ASTRONOMICAL_UNIT,
    color=(100,149,237),
    position=Vec3(1,0,0),
    velocity=Vec3(0,2*np.pi,0)
)

moon = Body(
    name="Moon",
    mass=3.7e-8,
    radius=1_737_000 / Units.ASTRONOMICAL_UNIT,
    color=(180,180,180),
    position=Vec3(1.2,0,0.1),
    velocity=Vec3(0,2*np.pi*0.9,0)
)


sim = Simulator(
    bodies=[
        sun,
        earth,
        moon
    ],
    timestep=1/365,
    steps=1000,
    integrator=RungeKutta,
    diagnostics=[
        TotalEnergy(),
        AngularMomentum()
    ]
)


result = sim.run(
    name="3D orbit"
)


viewer = SimulationViewer(
    result,
    ViewerConfig(
        view_size=1.5,
        camera_mode="fixed",
        fixed_center=(0,0,0),
        diagnostics=[
            DiagnosticView(
                diagnostic="Celotna Energija",
                entities=["vsi"]
            ),
            DiagnosticView(
                diagnostic="Vrtilna Količina",
                entities=["vsi"]
            )
        ]
    )
)

viewer.show()

Napoved se torej ujema z numeričnim rezultatom, zato sklepamo, da je simulator zaupljiv in ga lahko uporabljamo še za kakšne bolj komplicirane simulacije.