In fortran, my nbody code worked flawlessly. When I translated it over to python I messed up somewhere. The code still successfully reads in a .txt file of bodies and still prints out a .txt file that records their position and velocity over each timestep. However, all my bodies seem to slingshot into the great unknown. I know this because when comparing the plots (using mathematica) I notice that my fortran Jupiter generated plot (First image below) is perfectly elliptical and stays within the given distances expected of Jupiter and my python generated plot (second image below) goes all the way out to a magnitude of 10^15 which is absurd. I'm assuming it's a simple fix in the math somewhere but I just can't seem to figure it out.
import numpy as np
import pandas as pd
import os
import math as mt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
nbodies = int(input("Input the number of bodies: "))
dt = int(input("Input the number of timesteps: "))
ntimesteps = int(input("Max timestep: "))
G = 6.6743E-11
def getFile():
path = os.getcwd()
global fileName
fileName = input("What file? ")
inputFile = open(path + "\\" + str(fileName),"r")
return(inputFile)
def hasNumbers(inputString):
for char in inputString:
if(char == '-'):
return True
try:
float(char)
return True
except:
return False
def readFile(inputFile):
newdata = []
data = []
while True:
line = inputFile.readline()
#if the end of the file is reached
if not line:
break;
rawdata = line.split(" ")
for i in rawdata:
if hasNumbers(i) == True:
i = i.rstrip()
newdata.append(i)
if newdata:
data.append(newdata)
newdata = []
inputFile.close()
return(data)
def emptyFile(fileName, datatype):
path = os.getcwd()
outfile = open(path + "\\" + fileName[:-4] + "_" + datatype + ".txt", "w")
return(outfile)
inputFile = getFile()
Orbital = emptyFile(fileName, "Orbital")
nData = readFile(inputFile)
m = []
x = []
y = []
z = []
vx = []
vy = []
vz = []
for i in nData:
if hasNumbers(i) == True:
m.append(float(i[0]))
x.append(float(i[1]))
y.append(float(i[2]))
z.append(float(i[3]))
vx.append(float(i[4]))
vy.append(float(i[5]))
vz.append(float(i[6]))
mass = np.array(m)
xpos = np.array(x)
ypos = np.array(y)
zpos = np.array(z)
xvel = np.array(vx)
yvel = np.array(vy)
zvel = np.array(vz)
Fx = np.zeros(nbodies)
Fy = np.zeros(nbodies)
Fz = np.zeros(nbodies)
def CalcForce(axnew, aynew, aznew):
axnew = np.zeros(nbodies)
aynew = np.zeros(nbodies)
aznew = np.zeros(nbodies)
for i in range(nbodies):
for j in range(nbodies):
if i!=j:
xdiff = x[j] - x[i]
ydiff = y[j] - y[i]
zdiff = z[j] - z[i]
r2 = (xdiff**2)+(ydiff**2)+(zdiff**2)
C = (G*m[j]*m[i])/((r2)**(1.5))
Fx[i] = Fx[i]+C*(xdiff)
Fx[j] = Fx[j]-C*(xdiff)
Fy[i] = Fy[i]+C*(ydiff)
Fy[j] = Fy[j]-C*(ydiff)
Fz[i] = Fz[i]+C*(zdiff)
Fz[j] = Fz[j]-C*(zdiff)
axnew[i] = np.divide(Fx[i],m[i])
axnew[j] = np.divide(Fx[j],m[j])
aynew[i] = np.divide(Fy[i],m[i])
aynew[j] = np.divide(Fy[j],m[j])
aznew[i] = np.divide(Fz[i],m[i])
aznew[j] = np.divide(Fz[j],m[j])
else:
continue
return axnew, aynew, aznew
axnew, aynew, aznew = CalcForce(axnew, aynew, aznew)
for k in range(ntimesteps):
xpos = xpos + xvel * dt + 0.5 * axnew * dt**2
ypos = ypos + yvel * dt + 0.5 * aynew * dt**2
zpos = zpos + zvel * dt + 0.5 * aznew * dt**2
axold = axnew
ayold = aynew
azold = aznew
xvel = xvel + 0.5 * (axold + axnew) * dt
yvel = yvel + 0.5 * (ayold + aynew) * dt
zvel = zvel + 0.5 * (azold + aznew) * dt
np.savetxt(Orbital, np.transpose([np.arange(1, nbodies+1), np.full(nbodies, k*dt), xpos, ypos, zpos]),fmt="%2i %20i %15.5E %15.5E %15.5E")
Above I have my newest version of the python code. The "CalcForce.axnew" called accelerations worked when there was no If statement for i!=j. It doesn't anymore and I'm naïve as to why. Removing it and going back to trying to return the accelerations just returns 0s when called outside the definition.
Here's useful input to copy to a .txt file. It includes 11 bodies from the Sun to a comet (including Pluto which I use for a good time comparison). It's done in all 3 cartesian coordinates based on each body's orbital inclination with total linear momentum conserved.
m x y z vx vy vz
1.988e30 -4.490e5 0.000 0.000 0.000 -1.600e1 -0.4247
3.285e23 5.790e10 0.000 0.000 0.000 4.701e4 5.772e3
4.867e24 1.080e11 0.000 0.000 0.000 3.496e4 2.077e3
5.972e24 1.494e11 0.000 0.000 0.000 2.980e4 0.000
6.390e23 2.280e11 0.000 0.000 0.000 2.399e4 0.775e3
1.898e27 7.783e11 0.000 0.000 0.000 1.306e4 0.299e3
5.683e26 1.434e12 0.000 0.000 0.000 9.671e3 0.421e3
8.681e25 2.871e12 0.000 0.000 0.000 6.799e3 0.091e3
1.024e26 4.495e12 0.000 0.000 0.000 5.427e3 0.168e3
1.309e22 7.376e12 0.000 0.000 0.000 4.532e3 1.399e3
2.200e14 5.300e12 0.000 0.000 0.000 0.951e3 0.309e3
My fortran code for comparison is
Module nbody
integer, parameter:: n = 100
integer:: nbodies, nTimesteps, i, k
double precision:: m(n), x(n), y(n), z(n), vx(n), vy(n), vz(n)
double precision:: axold(n), ayold(n), azold(n), axnew(n), aynew(n), aznew(n)
double precision:: Fx(n), Fy(n), Fz(n)
double precision:: dt, G = 6.6743D-11, r2, C, S = 1.989D30
character:: infile*40, outfile*40
character(len=20), parameter :: fmt0 = "(1X,I15)",fmt1 = "(1X,E30.10)",fmt2 = "(1X,E30.10,/)"
end module
Program Gravitas
use nbody
print*,'input the number of bodies'
read*, nbodies
print*, 'input filename'
read*, infile
print*, 'output filename'
read*, outfile
open(unit=42, file=infile, status='unknown')
open(unit=43, file=outfile, status='unknown')
read(42,*)
Do i=1, nbodies
read(42,*) m(i), x(i), y(i), z(i), vx(i), vy(i), vz(i)
end do
print*, 'input timestep(dt)'
read*, dt
print*, 'Tmax'
read*, nTimesteps
Do i=1, nbodies
WRITE(43,fmt1,advance="no")real(k)*dt
WRITE(43,fmt0,advance="no")i
WRITE(43,fmt1,advance="no")x(i)
WRITE(43,fmt1,advance="no")y(i)
WRITE(43,fmt2,advance="no")z(i)
end do
call calcforces
Do k = 1, nTimesteps
x = x + vx * dt + 0.5 * axnew * dt**2
y = y + vy * dt + 0.5 * aynew * dt**2
z = z + vz * dt + 0.5 * aznew * dt**2
axold = axnew
ayold = aynew
azold = aznew
call calcforces
vx = vx + 0.5 * (axold + axnew) * dt
vy = vy + 0.5 * (ayold + aynew) * dt
vz = vz + 0.5 * (azold + aznew) * dt
Do i = 1, nbodies
WRITE(43,fmt1,advance="no")real(k)*dt
WRITE(43,fmt0,advance="no")i
WRITE(43,fmt1,advance="no")x(i)
WRITE(43,fmt1,advance="no")y(i)
WRITE(43,fmt2,advance="no")z(i)
end do
Subroutine calcforces
use nbody
Fx=0.0
Fy=0.0
Fz=0.0
do i=1, nbodies
do j=(i+1), (nbodies)
xdiff = x(j)-x(i)
ydiff = y(j)-y(i)
zdiff = z(j)-z(i)
r2 = (xdiff**2)+(ydiff**2)+(zdiff)**2
C = (G*m(j)*m(i))/((r2)**1.5)
Fx(i)=Fx(i)+C*(xdiff)
Fx(j)=Fx(j)-C*(xdiff)
Fy(i)=Fy(i)+C*(ydiff)
Fy(j)=Fy(j)-C*(ydiff)
Fz(i)=Fz(i)+C*(zdiff)
Fz(j)=Fz(j)-C*(zdiff)
end do
end do
axnew = Fx/m
aynew = Fy/m
aznew = Fz/m
end subroutine