Solar System Position and Velocity data

Question

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

Answer 1

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