Ich versuche, einen nbody-Simulator zu erstellen, und um ihn zu testen, möchte ich einige Daten aus der realen Welt. Ein guter Test wäre das Sonnensystem. Ich bin auf die JPL Horizons -Software gestoßen , aber wenn ich anscheinend nicht weiß, wie ich sie zum Laufen bringen kann, gibt sie eine Liste mit vielen Positionen für einen Körper, während ich eine Position für viele Körper haben möchte.
Ist die Horizons-Datenbank nicht das Richtige? können Sie mich in die richtige Richtung weisen. Wenn ja, können Sie mir einige Beispiele zeigen, wie man es zum Laufen bringt.
Idealerweise möchte ich nur eine Tabelle mit den Positionen der Planeten relativ zur Sonne in kartesischen Koordinaten
Sie haben nicht angegeben, wie genau Ihre Lösung sein soll, aber ich habe einmal eine n-Körper-Simulation in Python geschrieben, das Sonnensystem ist da, vielleicht hilft es?!
Wenn Sie denken, Sie könnten etwas davon gebrauchen, es aber nicht vollständig verstehen, lassen Sie es mich einfach wissen!
Für das Sonnensystem habe ich diese Werte verwendet:
m = np.array(([1],[1/6023600],[1/408524],[1/332946.038],[1/3098710],[1/1047.55],[1/3499],[1/22962],[1/19352])) #Masse der Objekte
r = np.array(([0,0],[0.4,0],[0,0.7],[1,0],[0,1.5],[5.2,0],[0,9.5],[19.2,0],[0,30.1]))
v = np.array(([0,0],[0,-np.sqrt(1/0.4)],[np.sqrt(1/0.7),0],[0,-1],[np.sqrt(1/1.5),0],[0,-np.sqrt(1/5.2)],[np.sqrt(1/9.5),0],[0,-np.sqrt(1/19.2)],[np.sqrt(1/30.1),0]))
Der gesamte Code (mit deutschen Anmerkungen und Euler-Integration (siehe Kommentare!)):
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 25 17:08:35 2018
r(t_k) known.
t_k = t0 + k * dt
t_(k+1/2) = t0 + (k+1/2) * dt
1. next Position
r(t_(k+1)) = r(t_k) + dt * dr(t_k+1/2)/dt
2. acceleration at Position at t_(k+1)
a(t_(k+1)) = -G*SUM[m_j * (r_i-r_j)/||r_i-r_j||³]
3. velocity at t_(k+3/2)
v(t_(k+3/2)) = v(t_(k+1/2)) + dt * a(t_(k+1))
@author: kalle
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
##Animation
fig, ax = plt.subplots()
xdata, ydata = [], []
ln, = plt.plot([], [], 'o', animated=True, color="red")
def init():
ax.set_xlim(-5,5)#-31,31 für Sonnensystem, sonst -5,5
ax.set_ylim(-5,5)#-31,31 für Sonnensystem
return ln,
#update plot
def update(i):
Pos=r_list[i]
for j in range(len(Pos)):
x_vals=Pos[j][0]
y_vals=Pos[j][1]
xdata.append(x_vals)
ydata.append(y_vals)
ln.set_data(xdata, ydata)
return ln,
#Simulationsbedingungen
ni = 150 #Anzahl Iterationsschritte
dt = 0.1 #Zeitinterval
#Vorgegebene initiale Objekteigenschaften
"""
Gravitationskraft propto 1/r^2
Massen im gesuchten Bsp. identisch. m1=m2=1
Für Kreisbahn muss Beschleunigung konstant sein, d.h. der Abstand der beiden
Massen relativ zueinander darf sich nicht ändern.
Auf die Masse selbst kommt es ebenfalls an, da zu leichte Körper zu weit
voneinander entfernt sonst einfach aneinander vorbeifliegen.
"""
""" verschiedene Systeme
Doppelsterne:
m = np.array(([1],[1])) #Masse der Objekte
r = np.array(([0,0.5],[0,-0.50]))#[0,1],[0,0]
v = np.array(([0,0.5],[0,-0.50]))#[-0.7071,0],[0.7071,0]
Doppelstern mit kleinem weit entfernten Begleitstern (leicht instabil, drift, Begleitstern zu nah / schwer)
m = np.array(([1],[1],[0.01])) #Masse der Objekte
r = np.array(([0,1],[0,0],[3,0]))
v = np.array(([-0.7071,0],[0.7071,0],[00,0.8]))
Sonnensystem:
Massen in Sonnenmassen,
Distanzen in AE,
v in Keplergeschwindigkeit v=sqrt(GM/r) mit M=Sonnenmasse, bzw. Masse des schweren Körpers.
m = np.array(([1],[1/6023600],[1/408524],[1/332946.038],[1/3098710],[1/1047.55],[1/3499],[1/22962],[1/19352])) #Masse der Objekte
r = np.array(([0,0],[0.4,0],[0,0.7],[1,0],[0,1.5],[5.2,0],[0,9.5],[19.2,0],[0,30.1]))
v = np.array(([0,0],[0,-np.sqrt(1/0.4)],[np.sqrt(1/0.7),0],[0,-1],[np.sqrt(1/1.5),0],[0,-np.sqrt(1/5.2)],[np.sqrt(1/9.5),0],[0,-np.sqrt(1/19.2)],[np.sqrt(1/30.1),0]))
Infinity Loop
m = np.array(([1],[1],[1])) #Masse der Objekte
r = np.array(([0.97000436,-0.24308753],[-0.97000436,0.24308753],[0,0]))
v = np.array(([0.93240737/2,0.86473146/2],[0.93240737/2,0.86473146/2],[-0.93240737,-0.86473146]))
"""
#Infinity Loop
m = np.array(([1],[1],[1])) #Masse der Objekte
r = np.array(([0.97000436,-0.24308753],[-0.97000436,0.24308753],[0,0]))
v = np.array(([0.93240737/2,0.86473146/2],[0.93240737/2,0.86473146/2],[-0.93240737,-0.86473146]))
r_list = []
r_list.append(r)
nk = len(m) #Anzahl der Objekte
#Initiale Werte
a = np.zeros((nk,2))
for i in range(0,nk):
for j in range(0,nk):
if i==j:
continue
a[i,0] = a[i,0] - m[j]*(r[i,0]-r[j,0])/np.power(np.linalg.norm(r[i,:]-r[j,:]),3)
a[i,1] = a[i,1] - m[j]*(r[i,1]-r[j,1])/np.power(np.linalg.norm(r[i,:]-r[j,:]),3)
v = v + 0.5*dt*a
#Iteration
for q in range(0,ni):
r = r + dt*v
r_list.append(r)
for i in range(0,nk):
a[i,:]=0
for j in range(0,nk):
if i==j:
continue
a[i,0] =a[i,0] - m[j]*(r[i,0]-r[j,0])/np.power(np.linalg.norm(r[i,:]-r[j,:]),3)
a[i,1] =a[i,1] - m[j]*(r[i,1]-r[j,1])/np.power(np.linalg.norm(r[i,:]-r[j,:]),3)
#a[:,:] = -a[:,:]
v = v + dt*a
ani = FuncAnimation(fig, update, frames=range(len(r_list)),
init_func=init, blit=True)
ani.save('basic_animation.mp4', fps=60, extra_args=['-vcodec', 'libx264'])
plt.show()
PM 2Ring
Benutzer7971589
Kalle