Edit: I made a test model with 2 bodies. The sun and earth or w/e. But I now have the problem of zero divisions every time that one of the coordinates is zero.
import numpy as np from numpy import exp, sqrt, pi, linspace from time import clock,time import pylab as pylab S1 = 1E3 #Mass Star M1 = 1 #Mass Planet Rs = np.array([0, 0]) #position R1 = np.array([0,200]) Vs = np.array([0,0]) #velocity V1 = np.array([10,0]) As = np.array([0,0]) #acceleration A1 = np.array([0,0]) dt = 0.1 #timestep N = 1000 #steps G = 1 #G const for i in range(1,N): deltaR = Rs - R1 Rsq = deltaR**2 F12 = -G*M1*S1/(Rsq) F21 = -F12 As = F21/S1 A1 = F12/M1 Vs += As*dt V1 += A1*dt Rs += Vs*dt + 0.5*As*dt**2 R1 += V1*dt + 0.5*A1*dt**2 x = [ Rs, R1] y = [ Rs, R1] pylab.ion() pylab.scatter( x, y ) pylab.axis([-200, 200, -200, 200]) pylab.plot pylab.clf()
So I want to simulate the solar system but want to start simple with one orbiting body. However, I never did anything like this before and was wondering if anyone here could give me some hints.
Clearly the first thing we need to use is the law of gravitation and Kepler. Then I'd like to throw everything into an array and just calculate the timesteps. However, I am not sure on how to get the position of the planet/moon as a function of time. (Can I just integrate the acceleration via the law of gravitation?)