# Solar System Position and Velocity data

I'm trying to create an nbody simulator and to test it I'd like some real world data. A good test would be the solar system. I've come across the JPL Horizons software but if cant seem to figure out how to get it to work, it gives a list of many positions for one body, whereas I would like one position for many bodies.

Is the Horizons database not the correct thing to be using? can you point me in the right direction. If it is, can you show me some examples of how to get it to work.

ideally I would just like a spreadsheet of the positions of the planets relative to the sun in Cartesian coordinates

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

• FWIW, your code uses simple Euler integration, which quickly accumulates errors unless the time step is tiny. You get much better results using a more accurate integrator, preferably a symplectic integrator (which conserves energy), like Leapfrog or Verlet. – PM 2Ring Nov 18 '18 at 4:39
• Appreciate you sharing this! Eventually though I ended up painstakingly copy pasting the information from the horizons website into a text file I could read in. Gonna probably have to write a script that does this for me if i wish to do it with more bodies! – user7971589 Nov 18 '18 at 5:46
• You are absolutely right, @2Ring ! Even using a simple Runge-Kutta Method would lower the error immensely! I'd like to highlight one of the most striking STABLE constellation: The infinity loop: Infinity Loop. 3 Objects are circling around each other in a stable 8 or $\infty$ form. It is analytically stable, but not found anywhere in the universe as far as I know. – kalle Nov 18 '18 at 12:11