Я пытаюсь создать симулятор nbody и протестировать его. Мне нужны данные из реального мира. Хорошим тестом будет Солнечная система. Я столкнулся с программным обеспечением JPL Horizons , но если я не могу понять, как заставить его работать, оно дает список многих позиций для одного тела, тогда как я хотел бы одну позицию для многих тел.
Разве база данных Horizons не подходит для использования? Можете ли вы указать мне в правильном направлении. Если это так, можете ли вы показать мне несколько примеров того, как заставить его работать.
в идеале я хотел бы просто таблицу положений планет относительно солнца в декартовых координатах
Вы не указали, насколько точным должно быть ваше решение, но я однажды написал симуляцию n тел на питоне, солнечная система там, может быть, это поможет ?!
Если вы думаете, что могли бы использовать что-то из этого, но не понимаете этого полностью, просто дайте мне знать!
Для солнечной системы я использовал следующие значения:
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]))
Весь код (с немецкими аннотациями и интеграцией Эйлера (см. комментарии!)):
# -*- 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Кольцо
пользователь7971589
улица