abb2A.jpg [ 69.49 KiB | 6472 mal betrachtet ]

# 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()