Omnimaga
General Discussion => Technology and Development => Computer Programming => Topic started by: ElementCoder on May 12, 2014, 09:25:35 am
-
I'm trying to make a two body simulation, but it doesn't seem to work :( Plotting x vs y should be circular or elliptical right? Am I choosing the wrong starting conditions or is my way of doing it just wrong?
#!/usr/bin/env python
'''
Two body simulation.
'''
from __future__ import division
import math
import matplotlib.pyplot as pyp
import numpy as np
# Gravitational constant (N*m^2*kg^-2)
G = 6.674e-11
# A giga year in seconds.
gyrts = math.pi * 1e7 * 1e9
def Fgrav(M, m, r):
''' Gravitational force of two bodies attracting each other.
'''
return -G * M * m / r**2
# Solar mass (kg)
M = 1.989e30
# Earth mass (kg)
m = 5.9721e24
# Initial conditions.
# Position (m).
x, y = 150e9, 0
X, Y = 0, 0
# Velocity (m/s).
vx, vy = 0, 30e3
Vx, Vy = 0, 0
time = 5*math.pi*1e7
t = 0
dt = 1000
xarr, yarr, rarr, tarr = [], [], [], []
while t < time:
# Update the positions.
x += vx * dt
y += vy * dt
# Calculate the separation vector.
dx = x - X
dy = y - Y
r2 = dx*dx + dy*dy + 1e-4
r = np.sqrt(r2)
#print 'dx, dy: ', dx, dy
# Calculate components of the gravitational force.
fgrav_x = Fgrav(M, m, r) * dx / r
fgrav_y = Fgrav(M, m, r) * dy / r
#print 'Fgrav: ', fgrav_x, fgrav_y
# Update the velocities.
vx += (fgrav_x / m) * dt
vy += (fgrav_y / m) * dt
Vx += (fgrav_x / M) * dt
Vy += (fgrav_y / M) * dt
#print 'vx, vy: ', vx, vy
X += Vx * dt
Y += Vy * dt
# Store all data.
xarr.append(dx); yarr.append(dy); rarr.append(r)
tarr.append(t)
t += dt
fig, (ax1, ax2) = pyp.subplots(2, 'nom')
ax1.plot(xarr, yarr)
ax1.set_xlabel('$x$ (m)'); ax1.set_ylabel('$y$ (m)')
ax2.plot(tarr, rarr)
ax2.set_xlabel('$t$ (s)'); ax2.set_ylabel('$r$ (m)')
fig.tight_layout()
pyp.show()
Edit: Updated with non-broken code.
-
The starting conditions are wrong (if you're simulating the earth and the sun):
- The initial positions: the real distance is ~150e9, but placing earth at (100e9, 50e9) gives us only a distance of 110e9. I'd choose (0, 0) for the sun and (150e9, 0) for earth.
- The initial speeds: why (3e3, 5e3)? If you choose (150e9, 0) as the starting position of earth, it's about (0, 30e3).
I don't have mathplotlib and numpy installed, so I can't test if it really works with my values. The results seem more reasonable, but I may be wrong :P
-
Those were just some test conditions.
I forgot a minus sign in the gravitational force. Everything works fine now :)
-
Ah lol, that would actually have been quite obvious from the results. Good to hear ;)