Dateianhang:
abb2A.jpg [ 69.49 KiB | 6470 mal betrachtet ]
Für die Lagrange-Punkte ist eine Lösung des Dreikörperproblems möglich. Das in Heft 82 vorgestellte Programm analysiert, ob dort Körper stabil positioniert werden können. Die Grafik beruht auf einem mitrotierenden Koordinatensystem.
Das Modul plot.py wird benötigt.
Code:
# LagrangeStabil.py - Analyse der Lagrangepunkte auf Stabibilität
from turtle import *
from random import *
from math import *
from plot import *
# globale Variable
G = 6.67408e-11 #Gravitationskonstante
mProAE = 149597870700 #Astronomische Einheit
sekProTag = 86400
true = 1
false = 0
# true: kleine Abweichung, um die Stabilität zu prüfen
ungenau = true
# statische Variable für rkf5
rErrMin = 0
rErrMax = 0
dt = 0
N = 7 # Anzahl der Körper
r0 = [[0 for i in range(7)] for j in range(7)]
v0 = [[0 for i in range(7)] for j in range(7)]
a0 = [[0 for i in range(7)] for j in range(7)]
r1 = [[0 for i in range(7)] for j in range(7)]
v1 = [[0 for i in range(7)] for j in range(7)]
a1 = [[0 for i in range(7)] for j in range(7)]
r2 = [[0 for i in range(7)] for j in range(7)]
v2 = [[0 for i in range(7)] for j in range(7)]
a2 = [[0 for i in range(7)] for j in range(7)]
r3 = [[0 for i in range(7)] for j in range(7)]
v3 = [[0 for i in range(7)] for j in range(7)]
a3 = [[0 for i in range(7)] for j in range(7)]
r4 = [[0 for i in range(7)] for j in range(7)]
v4 = [[0 for i in range(7)] for j in range(7)]
a4 = [[0 for i in range(7)] for j in range(7)]
r5 = [[0 for i in range(7)] for j in range(7)]
v5 = [[0 for i in range(7)] for j in range(7)]
a5 = [[0 for i in range(7)] for j in range(7)]
# r0_4: Ergebnis für RKF 4. Ordnung zur Fehlerberechnung
r0_4 = [[0 for i in range(7)] for j in range(7)]
m = [0 for i in range(7)]
def startwerte():
Y = 0
f = 1.0009548 # Korrekturfaktor für L4L5 bei Jupitermasse
T = 1 # simulierte Zeit in Jupiterumläufen
P = 100 # Anzahl Zwischenwerte
if (ungenau):
Y = 0.001 # bei L1 bis L3 den Ort und die Bahngeschwindigkeit leicht verschieben
f = 1.0 # bei L4 und L5 den Korrekturfaktor weglassen
T = 120 # Hufeisen von L3 erfordert 120 Umläufe
P = 6000 # ... und ordentlich Zwischenwerte
# Massen
m[0] = 1.988475415966536e+30
m[1] = 1.898568695e+27
m[2] = 0
m[3] = 0
m[4] = 0
m[5] = 0
m[6] = 0
# Orte und Geschwindigkeiten
r0[0][0] = -742448884.05 # Sonne
r0[0][1] = 0
r0[0][2] = 0
v0[0][0] = 0
v0[0][1] = -12.461410754
v0[0][2] = 0
r0[1][0] = 777607551115.94 # Jupiter
r0[1][1] = 0
r0[1][2] = 0
v0[1][0] = 0
v0[1][1] = 13051.520863369
v0[1][2] = 0
rJ = r0[1][0] # der Jupiterradius ist konstant in dieser Rechnung
vJ = v0[1][1] # die Bahngeschwindigkeit auch
r0[2][0] = 0.9332557970715 * rJ # L1
r0[2][1] = Y * r0[1][1]
r0[2][2] = 0
v0[2][0] = -Y * v0[1][1]
v0[2][1] = 0.9332557970715 * vJ
v0[2][2] = 0
r0[3][0] = 1.0698510256405 * rJ # L2
r0[3][1] = Y * r0[1][1]
r0[3][2] = 0
v0[3][0] = -Y * v0[1][1]
v0[3][1] = 1.0698510256405 * vJ
v0[3][2] = 0
r0[4][0] = -1.0013526135998 * rJ # L3
r0[4][1] = Y * r0[1][1]
r0[4][2] = 0
v0[4][0] = -Y * v0[1][1]
v0[4][1] = -1.0013526135998 * vJ
v0[4][2] = 0
# L4 und L5
r0[5][0] = f * (0.5 * (r0[0][0] + r0[1][0]) - r0[0][0])
r0[5][1] = f * ((r0[0][0] + r0[1][0]) * sqrt(3) / 2)
r0[5][2] = 0
v0[5][0] = -r0[5][1] * vJ / rJ
v0[5][1] = r0[5][0] * vJ / rJ
v0[5][2] = 0
# L5
r0[6][0] = r0[5][0]
r0[6][1] = -r0[5][1]
r0[6][2] = 0
v0[6][0] = -r0[6][1] * vJ / rJ
v0[6][1] = r0[6][0] * vJ / rJ
v0[6][2] = 0
return rJ, T , P
def sq(x):
return x * x
def abstand(r, p, q):
return sqrt(sq(r[p][0] - r[q][0]) + sq(r[p][1] - r[q][1]) + sq(r[p][2] - r[q][2]))
def beschleunigung(N, m, r, a):
for p in range(N): # alle Beschleunigungen löschen
for k in range(3):
a[p][k] = 0
for p in range(N): # alle Planeten
for q in range(p + 1, N): # gegen alle anderen Planeten
R = abstand(r, p, q)
R3 = R * R * R
for k in range(3):
A = G * (r[q][k] - r[p][k]) / R3
a[p][k]+=m[q] * A
a[q][k]-=m[p] * A
def rkf5(tEnd, init):
global dt, rErrMin, rErrMax
t = 0
while(t < tEnd):
t = t + dt
beschleunigung(N, m, r0, a0)
for i in range(N):
for k in range(3):
r1[i][k] = r0[i][k] + 1 / 4 * dt * v0[i][k]
v1[i][k] = v0[i][k] + 1 / 4 * dt * a0[i][k]
beschleunigung(N, m, r1, a1)
for i in range(N):
for k in range(3):
r2[i][k] = r0[i][k] + 3 / 32 * dt * v0[i][k] + 9 / 32 * dt * v1[i][k]
v2[i][k] = v0[i][k] + 3 / 32 * dt * a0[i][k] + 9 / 32 * dt * a1[i][k]
beschleunigung(N, m, r2, a2)
for i in range(N):
for k in range(3):
r3[i][k] = r0[i][k] + 1932 / 2197 * dt * v0[i][k] + -7200 / 2197 * dt * v1[i][k] + \
7296 / 2197 * dt * v2[i][k]
v3[i][k] = v0[i][k] + 1932 / 2197 * dt * a0[i][k] + -7200 / 2197 * dt * a1[i][k] + \
7296 / 2197 * dt * a2[i][k]
beschleunigung(N, m, r3, a3)
for i in range(N):
for k in range(3):
r4[i][k] = r0[i][k] + 439 / 216 * dt * v0[i][k] - 8 * dt * v1[i][k] + 3680 / 513 * \
dt * v2[i][k] - 845 / 4104 * dt * v3[i][k]
v4[i][k] = v0[i][k] + 439 / 216 * dt * a0[i][k] - 8 * dt * a1[i][k] + 3680 / 513 * \
dt * a2[i][k] - 845 / 4104 * dt * a3[i][k]
beschleunigung(N, m, r4, a4)
for i in range(N):
for k in range(3):
r5[i][k] = r0[i][k] - 8 / 27 * dt * v0[i][k] + 2 * dt * v1[i][k] - 3544 / 2565 * \
dt * v2[i][k] + 1859 / 4104 * dt * v3[i][k] - 11 / 40 * dt * v4[i][k]
v5[i][k] = v0[i][k] - 8 / 27 * dt * a0[i][k] + 2 * dt * a1[i][k] - 3544 / 2565 * \
dt * a2[i][k] + 1859 / 4104 * dt * a3[i][k] - 11 / 40 * dt * a4[i][k]
beschleunigung(N, m, r5, a5)
for i in range(N):
for k in range(3):
r0_4[i][k] = r0[i][k]
# Ergebniss 4. Ordnung zur Fehlerberechnung
r0_4[i][k]+= dt * (25 / 216 * v0[i][k] + 0 * v1[i][k] + 1408 / 2565 * v2[i][k] + \
2197.0 / 4104 * v3[i][k] - 1 / 5 * v4[i][k] + 0 * v5[i][k])
# eigentliches Ergebnis 5. Ordnung
r0[i][k]+= dt * (16 / 135 * v0[i][k] + 0 * v1[i][k] + 6656 / 12825 * v2[i][k] + \
28561 / 56430 * v3[i][k] - 9 / 50 * v4[i][k] + 2 / 55 * v5[i][k])
v0[i][k]+= dt * (16 / 135 * a0[i][k] + 0 * a1[i][k] + 6656 / 12825 * a2[i][k] + \
28561 / 56430 * a3[i][k] - 9 / 50 * a4[i][k] + 2 / 55 * a5[i][k])
# Fehler ausrechnen. Benutzung erfordert eine while- statt einer
# for-Schleife
err = 0
for i in range(N): # welcher Körper gibt den größten Fehler?
localErr = 0
for k in range(3):
localErr+=sq(r0_4[i][k] - r0[i][k])
localErr = sqrt(localErr)
if (localErr > err):
err = localErr
if (init > 0): # Fehler wird in der ersten Runde so gesetzt,
rErrMin = err / 2 # wie es dem ersten dt entspricht
rErrMax = err
init = false
if (err < rErrMin): # Schrittweite verändern, falls der Fehler
dt*=1.2 # außerhalb der Grenzen
elif (err > rErrMax):
dt/=1.2
if (dt / sekProTag < 0.01):
Err = input("Schrittweitensteuerung versagt: dt erhöhen")
exit()
return t
# Hauptprogramm
rJ, T, P = startwerte()
initKoor(-1.3,1.3,0.2, -1.1,1.1,0.2)
dt = 2.5 # 2,5 Tage
# den anderen Teil der Startwerte in SI umwandeln
dt = dt * sekProTag
UJup = 2 * pi * rJ / v0[1][1]
T = T * UJup
# Punkte für Jupiter und Sonne
Plot(r0[0][0] / rJ,r0[0][1] / rJ,20, "gold")
Plot(r0[1][0] / rJ,r0[1][1] / rJ,8,"coral")
t = 0
init = true
while (t < T): # den gesamten Zeitbereich rechnen
t = t + rkf5(T / P, init) # einen Simulationsschritt
init = false # intern Fehlergrenze nur 1x setzen
xJ = r0[1][0]
yJ = r0[1][1]
for i in range(2,7): # Sonne, Jupiter, und L1-L5
x = r0[i][0]
y = r0[i][1]
# in mitrotierene Koordinaten umrechnen (Skalarprodukt)
rP = sqrt(x * x + y * y)
xP = xJ * x + yJ * y
xP = xP / rJ / rJ
yP = -yJ * x + xJ * y
yP = yP / rJ / rJ
if i==2: # L1
Plot(xP,yP,3,"blue")
if i==3: # L2
Plot(xP,yP,3,"red")
if i==4: # L3
Plot(xP,yP,3,"magenta")
if i>4:
Plot(xP,yP,3,"green")
update() # Zwischenergebnis anzeigen
update()
Registrieren
Anmelden
FAQ
Suche
