在围绕中心物体的引力场中计算物体轨迹的跳跃算法 (Python 3.8.2)

Leapfrog algorithm to compute a objects trajectory in a Gravitational field around a central body (Python 3.8.2)

提问人:Prometheus 提问时间:3/2/2021 最后编辑:Lutz LehmannPrometheus 更新时间:3/12/2021 访问量:790

问:

我几乎删除了最后一个代码并开始新的代码。我添加了一个名为 Object 的新类,它取代了名为 body_1 和 body_2 的列表。此外,所有计算现在都是在 Object 类中完成的。以前的大多数现有问题都通过这个过程得到了解决,但仍有一个问题是抵制的。我相信它在 StartVelocity() 函数中,该函数创建启动 Leapfrog 算法所需的 v1/2。这应该给我一个地球静止轨道,但可以清楚地看到,卫星在穿越地球后很快就逃脱了。enter image description here

代码是:

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]]),
                        np.array([[0.0], [0.0], [0.0]]), 
                        5.9722 * 10**24)

        body_2 = Object("Satelite", "b", "r",
                        np.array([[42164.0], [0.0], [0.0]]),
                        np.array([[0.0], [3075.4], [0.0]]),
                        5000.0)
        
        self.bodies = [body_1, body_2]

    def ComputePath(self, time_limit, time_step):

        time_range = np.arange(0, time_limit, 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, Y, 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, other_bodies, 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 才能得到米,但如果我这样做,它只会直线飞行,并且不会有引力的迹象。

python 数值方法 轨道力学 verlet 积分

评论

0赞 Sam Mason 3/2/2021
这篇文章可能更适合 codereview.stackexchange.com
0赞 Prometheus 3/2/2021
codereview.stackexchange.com 要求我有一个正确工作的代码来回答一个问题@SamMason

答:

1赞 Lutz Lehmann 3/11/2021 #1

常数以千米秒为单位。然而,卫星轨道的半径只能以公里为单位有意义,否则轨道将在地核内。然后以米/秒为单位的速度给出一个近圆形的轨道,偏心率可以忽略不计。(来自开普勒定律量的 math.SE 问题代码G)

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)

给出一个很好的圆形轨道

enter image description here