问题描述
我几乎删除了最后的代码并开始了新的代码。我添加了一个名为 Object 的新类,它替代了名为 body_1 和 body_2 的列表。此外,所有计算现在都在 Object 类中完成。大多数以前存在的问题都通过这个过程得到了解决,但仍然存在一个问题。我相信它在 StartVeLocity() 函数中,它创建了启动 Leapfrog 算法所需的 v1/2。这应该给我一个地球静止轨道,但清晰可见,卫星在放大地球后很快就会逃逸。
代码是:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from object import Object
import numpy as np
class Simulation:
def __init__(self):
# Index: 0 Name,1 Position,2 VeLocity,3 Mass
body_1 = Object("Earth","g","r",np.array([[0.0],[0.0],[0.0]]),5.9722 * 10**24)
body_2 = Object("Satelite","b",np.array([[42164.0],[3075.4],5000.0)
self.bodies = [body_1,body_2]
def ComputePath(self,time_limit,time_step):
time_range = np.arange(0,time_step)
for body in self.bodies:
body.StartVeLocity(self.bodies,time_step)
for T in time_range:
for body in self.bodies:
body.Leapfrog(self.bodies,time_step)
def PlotObrit(self):
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
for body in self.bodies:
body.ReshapePath()
X,Y,Z = [],[],[]
for position in body.path:
X.append(position[0])
Y.append(position[1])
Z.append(position[2])
ax.plot(X,Z,f"{body.linecolor}--")
for body in self.bodies:
last_pos = body.path[-1]
ax.plot(last_pos[0],last_pos[1],last_pos[2],f"{body.bodycolor}o",label=body.name)
ax.set_xlabel("x-Axis")
ax.set_ylabel("y-Axis")
ax.set_zlabel("z-Axis")
ax.legend()
fig.savefig("Leapfrog.png")
if __name__ == "__main__":
sim = Simulation()
sim.ComputePath(0.5,0.01)
sim.PlotObrit()
import numpy as np
class Object:
def __init__(self,name,bodycolor,linecolor,pos_0,vel_0,mass):
self.name = name
self.bodycolor = bodycolor
self.linecolor = linecolor
self.position = pos_0
self.veLocity = vel_0
self.mass = mass
self.path = []
def StartVeLocity(self,other_bodies,time_step):
force = self.GetForce(other_bodies)
self.veLocity += (force / self.mass) * time_step * 0.5
def Leapfrog(self,time_step):
self.position += self.veLocity * time_step
self.veLocity += (self.GetForce(other_bodies) / self.mass) * time_step
self.path.append(self.position.copy())
def GetForce(self,other_bodies):
force = 0
for other_body in other_bodies:
if other_body != self:
force += self.Force(other_body)
return force
def Force(self,other_body):
G = 6.673 * 10**-11
dis_vec = other_body.position - self.position
dis_mag = np.linalg.norm(dis_vec)
dir_vec = dis_vec / dis_mag
for_mag = G * (self.mass * other_body.mass) / dis_mag**2
for_vec = for_mag * dir_vec
return for_vec
def ReshapePath(self):
for index,position in enumerate(self.path):
self.path[index] = position.reshape(3).tolist()
我知道身体 2 的位置必须乘以 1000 才能得到米,但如果我这样做,它只会直线飞行,而且不会有任何重力迹象。
解决方法
常数 G
的单位是 kg-m-sec。然而,卫星轨道的半径只能以公里为单位,否则轨道将在地球核心内部。然后以米/秒为单位的速度给出了一个近乎圆形的轨道,偏心率可以忽略不计。 (来自 Kepler law quantities 上的 math.SE 问题的代码)
import math as m
G = 6.673e-11*1e-9 # km^3 s^-2 kg^-1
M_E = 5.9722e24 # kg
R_E = 6378.137 # km
R_sat = 42164.0 # km from Earth center
V_sat = 3075.4/1000 # km/s
theta = 0
r0 = R_sat
dotr0 = V_sat*m.sin(theta)
dotphi0 = -V_sat/r0*m.cos(theta)
R = (r0*V_sat*m.cos(theta))**2/(G*M_E)
wx = R/r0-1; wy = -dotr0*(R/(G*M_E))**0.5
E = (wx*wx+wy*wy)**0.5; psi = m.atan2(wy,wx)
T = m.pi/(G*M_E)**0.5*(R/(1-E*E))**1.5
print(f"orbit constants R={R} km,E={E},psi={psi} rad")
print(f"above ground: min={R/(1+E)-R_E} km,max={R/(1-E)-R_E} km")
print(f"T={2*T} sec,{T/1800} h")
带输出
orbit constants R=42192.12133271948 km,E=0.0006669512550867562,psi=-0.0 rad
above ground: min=35785.863 km,max=35842.14320159004 km
T=86258.0162673565 sec,23.960560074265697 h
对于r(phi)=R/(1+E*cos(phi-psi))
在您的代码中实现这些更改并调用
sim.ComputePath(86e+3,600.0)
给出一个很好的圆形轨道