Estoy tratando de crear un simulador de nbody y, para probarlo, me gustaría obtener algunos datos del mundo real. Una buena prueba sería el sistema solar. Me encontré con el software JPL Horizons , pero si parece que no puedo averiguar cómo hacerlo funcionar, da una lista de muchas posiciones para un cuerpo, mientras que me gustaría una posición para muchos cuerpos.
¿La base de datos de Horizons no es lo correcto para usar? Me puede apuntar en la dirección correcta. Si es así, ¿puede mostrarme algunos ejemplos de cómo hacerlo funcionar?
idealmente me gustaría una hoja de cálculo de las posiciones de los planetas en relación con el sol en coordenadas cartesianas
No dijiste qué tan precisa será tu solución, pero una vez escribí una simulación de n cuerpos en Python, el sistema solar está ahí, ¿tal vez ayude?
Si crees que podrías usar algo de esto pero no lo entiendes completamente, ¡avísame!
Para el sistema solar usé esos valores:
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]))
El código completo (con anotaciones alemanas e integración de Euler (¡ver comentarios!)):
# -*- 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 2 Anillo
usuario7971589
kale