You didn't state how accurate your solution shall be, but I wrote once a n-body simulation in python, the solar system is there, maybe it helps?!
If you think you could use something of this but don't understand it fully, just let me know!
For the solar system I used those values:
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]))
The whole code (with german annotations and Euler integration (see comments!)):
# -*- 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()