Как исправить мой код для моделирования реалистичной задачи трех тел? Мой код отображает прямые линии, которые являются ⇐ Python

[Anonymous]

Это мой код: (я хочу построить в 3D хаотические и символические пути взаимодействия трех тел с учетом любых входных данных для массы, радиуса и начального положения трех тел) Я продолжаю получать прямые линии в качестве результатов на построенном 3D-графике. Может ли кто-нибудь сказать мне, что я сделал не так? Проблема здесь в физике, математике или программировании?
import math
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
import numpy as np

class Vector:
def __init__(self, x=0, y=0,z=0):
self.x = x
self.y = y
self.z = z

def showVector(self):
print(self.x,self.y,self.z)
def getVector(self):
vectors=[self.x,self.y,self.z]
return vectors
def getVecti(self):
return self.x
def getVectj(self):
return self.y
def getVectk(self):
return self.z
def magnitude(self):
m=math.sqrt((self.x)**2 + (self.y)**2 + (self.z)**2)
return m
def changeVector(self,vector):
self.x=vector.x
self.y=vector.y
self.z=vector.z
def addVector(self,vector1,vector2):
self.x=vector1.x+vector2.x
self.y=vector1.y+vector2.y
self.z=vector1.z+vector2.z
self.changeVector(Vector(self.x,self.y,self.z))
def timesVector_withNum(self,vector,Num):
self.x=Num*vector.x
self.y=Num*vector.y
self.z=Num*vector.z
self.changeVector(Vector(self.x,self.y,self.z))
def subtractVector(self,v1,v2):
self.x=v1.x-v2.x
self.y=v1.y-v2.y
self.z=v1.z-v2.z
self.changeVector(Vector(self.x, self.y, self.z))
#def crossProduct(self,i,j,k,x,y,):
# self.x=j*z-k*y
# self.y=k*x-i*z
# self.z=i*y-j*x
#self.changeVector(Vector(self.x, self.y, self.z))
#return Vector(self.x, self.y, self.z)
def Vector_divideNum(self,vector,n):
self.x=vector.x/n
self.y=vector.y/n
self.z=vector.z/n
self.changeVector(Vector(self.x,self.y,self.z))

def check_float(variable):
check=True
print(variable, end='')
variable=input()
while check:
try:
variable=float(variable)
check=False
except:
variable=input('Invalid input, enter another input: ')
return variable

class objects:
def __init__(self, mass, position, radius, g, velocity=Vector(), acceleration=Vector()):
self.mass=mass
self.position=position
self.r=radius
self.g=g
self.v=velocity
self.a=acceleration

def getMass(self):
return self.mass
def getPosition(self):
return self.position
def get_r(self):
return self.r
def get_g(self):
return self.g
def getVelocity(self):
return self.v
def getAcceleration(self):
return self.a

G=6.67430*10**(-11)

timeInterval=1000

#object 1: (PlanetA)
mass1="Enter mass of object 1: "
mass1=check_float(mass1)

r1="Enter radius of object 1:"
r1=check_float(r1)

p1='Enter centre position i value: '#p for position
p1=check_float(p1)

p2='Enter centre position j value: '
p2=check_float(p2)

p3='Enter centre position k value: '
p3=check_float(p3)

position1=Vector(p1,p2,p3)

#object 2: (PlanetB)

mass2="Enter mass of object 2: "
mass2=check_float(mass2)

r2="Enter radius of object 2: "
r2=check_float(r2)

p1='Enter centre position i value: '
p1=check_float(p1)

p2='Enter centre position j value: '
p2=check_float(p2)

p3='Enter centre position k value: '
p3=check_float(p3)

position2=Vector(p1,p2,p3)

#Object 3 : PlanetC

mass3="Enter mass of object 3: "
mass3=check_float(mass3)

r3="Enter radius of object 3: "
r3=check_float(r3)

p1='Enter centre position i value: '
p1=check_float(p1)

p2='Enter centre position j value: '
p2=check_float(p2)

p3='Enter centre position k value: '
p3=check_float(p3)

position3=Vector(p1,p2,p3)

#F_21=-G (m_1 m_2)/|r_21 |^2 r ̂_21
#r ̂_21=(r_2-r_1)/mod(r_2-r_1)

def Force_on_2(mass1,mass2,mass3,position1,position2,position3): #Here I used equation: Force on 2 = -G*mass1*mass2*vector distance/(distance cubed)-G*mass3*mass2*vector distance/(distance cubed)

total=Vector()
F_21=Vector()
F_23=Vector()
r=Vector()
r.subtractVector(position2,position1)
r_magnitude=r.magnitude()
if r_magnitude!=0:
F_21.timesVector_withNum(r,-G*(mass1*mass2)/(r_magnitude**3))

r=Vector()
r.subtractVector(position2,position3)
r_magnitude=r.magnitude()
if r_magnitude!=0:
F_23.timesVector_withNum(r,-G*(mass3*mass2)/(r_magnitude**3))
total.addVector(F_21,F_23)
return total

planetA=objects(mass1, position1,r1,G*mass1/(r1**2),Vector()) #You can try any initial velocity vector you want by entering numbers in the form Vector(a,b,c), where (a,b,c) is velocity vector

planetB=objects(mass2, position2,r2,G*mass2/(r2**2),Vector())

planetC=objects(mass3, position3,r3,G*mass3/(r3**2),Vector())

FonA=Vector() #Force on A

Acceleration_A=Vector()
Acceleration_B=Vector()
Acceleration_C=Vector()
x1=[planetA.position.x]
y1=[planetA.position.y]
z1=[planetA.position.z]

x2=[planetB.position.x]
y2=[planetB.position.y]
z2=[planetB.position.z]

x3=[planetC.position.x]
y3=[planetC.position.y]
z3=[planetC.position.z]

for i in range(1000000):

FonA=Force_on_2(planetB.mass,planetA.mass,planetC.mass,planetB.position,planetA.position,planetC.position)
Initial_Acceleration_A=Vector(FonA.x/planetA.mass, FonA.y/planetA.mass, FonA.z/planetA.mass)

planetA.position.changeVector(Vector(planetA.position.x+planetA.v.x*timeInterval+0.5*Initial_Acceleration_A.x*(timeInterval**2),
planetA.position.y+planetA.v.y*timeInterval+0.5*Initial_Acceleration_A.y*(timeInterval**2),
planetA.position.z+planetA.v.z*timeInterval+0.5*Initial_Acceleration_A.z*(timeInterval**2)))

FonA=Force_on_2(planetB.mass,planetA.mass,planetC.mass,planetB.position,planetA.position,planetC.position)
Final_Acceleration_A=Vector(FonA.x/planetA.mass, FonA.y/planetA.mass, FonA.z/planetA.mass)
planetA.v.changeVector(Vector(planetA.v.x+0.5*(Initial_Acceleration_A.x+Final_Acceleration_A.x)*timeInterval,
planetA.v.y+0.5*(Initial_Acceleration_A.y+Final_Acceleration_A.y)*timeInterval,
planetA.v.z+0.5*(Initial_Acceleration_A.z+Final_Acceleration_A.z)*timeInterval))

if i%1000==0:
x1.append(planetA.position.x)
y1.append(planetA.position.y)
z1.append(planetA.position.z)

#B
FonB=Force_on_2(planetA.mass,planetB.mass,planetC.mass,planetA.position,planetB.position,planetC.position)
Initial_Acceleration_B=Vector(FonB.x/planetB.mass, FonB.y/planetB.mass, FonB.z/planetB.mass)

planetB.position.changeVector(Vector(planetB.position.x+planetB.v.x*timeInterval+0.5*Initial_Acceleration_B.x*(timeInterval**2),
planetB.position.y+planetB.v.y*timeInterval+0.5*Initial_Acceleration_B.y*(timeInterval**2),
planetB.position.z+planetB.v.z*timeInterval+0.5*Initial_Acceleration_B.z*(timeInterval**2)))

FonB=Force_on_2(planetA.mass,planetB.mass,planetC.mass,planetA.position,planetB.position,planetC.position)
Final_Acceleration_B=Vector(FonB.x/planetB.mass, FonB.y/planetB.mass, FonB.z/planetB.mass)
planetB.v.changeVector(Vector(planetB.v.x+0.5*(Initial_Acceleration_B.x+Final_Acceleration_B.x)*timeInterval,
planetB.v.y+0.5*(Initial_Acceleration_B.y+Final_Acceleration_B.y)*timeInterval,
planetB.v.z+0.5*(Initial_Acceleration_B.z+Final_Acceleration_B.z)*timeInterval))

if i%1000==0:
x2.append(planetB.position.x)
y2.append(planetB.position.y)
z2.append(planetB.position.z)

FonC=Force_on_2(planetA.mass,planetC.mass,planetB.mass,planetA.position,planetC.position,planetB.position)
Initial_Acceleration_C=Vector(FonC.x/planetC.mass, FonC.y/planetC.mass, FonC.z/planetC.mass)

planetC.position.changeVector(Vector(planetC.position.x+planetC.v.x*timeInterval+0.5*Initial_Acceleration_C.x*(timeInterval**2),
planetC.position.y+planetC.v.y*timeInterval+0.5*Initial_Acceleration_C.y*(timeInterval**2),
planetC.position.z+planetC.v.z*timeInterval+0.5*Initial_Acceleration_C.z*(timeInterval**2)))

FonC=Force_on_2(planetA.mass,planetC.mass,planetB.mass,planetA.position,planetC.position,planetB.position)
Final_Acceleration_C=Vector(FonC.x/planetC.mass, FonC.y/planetC.mass, FonC.z/planetC.mass)
planetC.v.changeVector(Vector(planetC.v.x+(Initial_Acceleration_C.x+Final_Acceleration_C.x)*timeInterval/2,
planetC.v.y+(Initial_Acceleration_C.y+Final_Acceleration_C.y)*timeInterval/2,
planetC.v.z+(Initial_Acceleration_C.z+Final_Acceleration_C.z)*timeInterval/2))

if i%1000==0:
x3.append(planetC.position.x)
y3.append(planetC.position.y)
z3.append(planetC.position.z)

fig=plt.figure()
ax=plt.axes(projection='3d')
ax.plot(np.array(x1),np.array(y1),np.array(z1),color='red',lw=1,markersize=3,alpha=0.9)
ax.plot(np.array(x2),np.array(y2),np.array(z2),color='blue',lw=1,markersize=3,alpha=0.9)
ax.plot(np.array(x3),np.array(y3),np.array(z3),color='black',lw=1,markersize=3,alpha=0.9)
plt.show()


Подробнее здесь: https://stackoverflow.com/questions/798 ... -plots-str