Oliver Natt: Physik mit Python.

Alle Programme, die in diesem Buch besprochen werden, und die Programme zu den Musterlösungen der Übungsaufgaben können Sie auch als ein zip-Archiv herunterladen.
Programm 6.10 Dynamik/planetenbahn3.py
"""Simulation einer Planetenbahn unter der Annahme, dass die
Sonne ruht. Die Anfangsbedingungen wurden so gewählt, dass eine
stark elliptische Umlaufbahn entsteht. """

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.animation
import scipy.integrate

# Konstanten 1 Tag [s] und 1 Jahr [s].
tag = 24 * 60 * 60
jahr = 365.25 * tag

# Eine Astronomische Einheit [m].
AE = 1.495978707e11

# Skalierungsfaktor für die Darstellung der Beschleunigung
# [AE / (m/s²)] und Geschwindigkeit [AE / (m/s)]
scal_a = 20
scal_v = 1e-5

# Simulationsdauer T und dargestellte Schrittweite dt [s].
T = 1 * jahr
dt = 1 * tag

# Masse der Sonne M [kg].
M = 1.9885e30

# Gravitationskonstante [m³ / (kg * s²)].
G = 6.674e-11

# Anfangsposition des Planeten [m].
r0 = np.array([152.10e9, 0.0])

# Anfangsgeschwindigkeit des Planeten [m/s].
v0 = np.array([0.0, 15e3])


def dgl(t, u):
    r, v = np.split(u, 2)
    a = - G * M * r / np.linalg.norm(r) ** 3
    return np.concatenate([v, a])


# Lege den Zustandsvektor zum Zeitpunkt t=0 fest.
u0 = np.concatenate((r0, v0))

# Löse die Bewegungsgleichung.
result = scipy.integrate.solve_ivp(dgl, [0, T], u0, rtol=1e-9,
                                   dense_output=True)
t_s = result.t
r_s, v_s = np.split(result.y, 2)

# Berechne die Interpolation auf einem feinen Raster.
t = np.arange(0, np.max(t_s), dt)
r, v = np.split(result.sol(t), 2)

# Erzeuge eine Figure und eine Axes.
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('x [AE]')
ax.set_ylabel('y [AE]')
ax.set_xlim(-0.2, 1.1)
ax.set_ylim(-0.6, 0.6)
ax.set_aspect('equal')
ax.grid()

# Plotte die Bahnkurve der Himmelskörper.
ax.plot(r_s[0] / AE, r_s[1] / AE, '.b')
ax.plot(r[0] / AE, r[1] / AE, '-b')

# Erzeuge Punktplots, für die Positionen der Himmelskörper.
planet, = ax.plot([0], [0], 'o', color='red')
sonne, = ax.plot([0], [0], 'o', color='gold')

# Erzeuge zwei Pfeile für die Beschleunigungsvektoren.
style = mpl.patches.ArrowStyle.Simple(head_length=6,
                                      head_width=3)
arrow_a = mpl.patches.FancyArrowPatch((0, 0), (0, 0),
                                      color='black',
                                      arrowstyle=style)
arrow_v = mpl.patches.FancyArrowPatch((0, 0), (0, 0),
                                      color='red',
                                      arrowstyle=style)

# Füge die Pfeil zur Axes hinzu.
ax.add_artist(arrow_a)
ax.add_artist(arrow_v)

# Füge ein Textfeld für die Angabe der verstrichenen Zeit hinzu.
text_t = ax.text(0.01, 0.95, '', color='blue',
                 transform=ax.transAxes)


def update(n):
    # Aktualisiere die Position des Himmelskörpers.
    planet.set_data(r[:, n].reshape(-1, 1) / AE)

    # Berechne die Momentanbeschleunigung und aktualisiere die
    # Vektorpfeile.
    u = np.concatenate((r[:, n], v[:, n]))
    u_punkt = dgl(t[n], u)
    a = np.split(u_punkt, 2)[1]

    arrow_a.set_positions(r[:, n] / AE,
                          r[:, n] / AE + scal_a * a)
    arrow_v.set_positions(r[:, n] / AE,
                          r[:, n] / AE + scal_v * v[:, n])

    # Aktualisiere das Textfeld für die Zeit.
    text_t.set_text(f't = {t[n] / tag:.0f} d')

    return planet, arrow_v, arrow_a, text_t


# Erzeuge das Animationsobjekt und starte die Animation.
ani = mpl.animation.FuncAnimation(fig, update, interval=30,
                                  frames=t.size)
plt.show()