Oliver Natt: Physik mit Python.

Lösung zu Aufgabe 6.6

Dynamik/Loesungen/flaechensatz.py
"""Animation des keplerschen Flächensatzes. """

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

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

# Anzahl der Zeitschritte, die für die Darstellung der Fläche
# des Fahrstrahls verwendet wird
dN = 40

# Massen der Sonne M [kg].
M = 1.9889e30

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

# Anfangspositionen des Planeten [m].
r0 = np.array([152.10e9, 0.0])
v0 = np.array([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 numerisch.
result = scipy.integrate.solve_ivp(dgl, [0, T], u0, rtol=1e-9,
                                   t_eval=np.arange(0, T, dt))
t = result.t
r, v = np.split(result.y, 2)

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

# Plotte die Bahnkurve des Himmelskörpers.
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 ein Polygon zur Darstellung der überstrichenen Fläche
# und füge dieses der Axes hinzu.
flaeche = mpl.patches.Polygon([[0, 0], [0, 0]], closed=True,
                              alpha=0.5, facecolor='red')
ax.add_artist(flaeche)

# Erzeuge zwei Textfelder für die Angabe der verstrichenen Zeit
# und die berechnete Fläche.
text_t = ax.text(0.01, 0.95, '', color='black',
                 transform=ax.transAxes)
text_A = ax.text(0.01, 0.90, '', color='black',
                 transform=ax.transAxes)


def polygon_flaeche(x, y):
    """Berechne die Fläche eines Polygons mit den angegebenen
    Koordinaten mithilfe der gaußschen Trapezformel. """
    return 0.5 * abs((y + np.roll(y, 1)) @ (x - np.roll(x, 1)))


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

    # Aktualisiere des Polygon und die Angabe der Fläche.
    if n >= dN:
        # Erzeuge ein (dN + 2) x 2 - Array. Als ersten Punkt
        # enthält dies die Position (0, 0) der Sonne und die
        # weiteren Punkte sind die dN + 1 Punkte der Bahnkurve
        # des Planeten.
        xy = np.zeros((dN + 2, 2))
        xy[1:, :] = r[:, (n - dN):(n + 1)].T / AE
        flaeche.set_xy(xy)

        # Berechne die Fläche des Polygons und gebe diese aus.
        A = polygon_flaeche(xy[:, 0], xy[:, 1])
        text_A.set_text(f'A = {A:.2e} AE²')
    else:
        # Zu Beginn der Animation kann noch keine Fläche
        # dargestellt werden.
        flaeche.set_xy([[0, 0], [0, 0]])
        text_A.set_text(f'')

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

    return planet, text_t, text_A, flaeche


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

plt.show()
Aufgaben

Lösung zu Aufgabe 6.6
