轨道轨道模拟器

Question

我已经编写了一个简单的程序来在地球-月球系统中进行轨迹模拟，它还有很长的路要走--我正在努力使它更加面向类，并且正在考虑实现一个比直接时间步长积分法更好的系统。

我想把重点放在让它更面向班级。我不知道应该在什么地方划一个好的界限，究竟有多少东西应该由类/函数驱动，有多少应该在主块中。我也在寻找关于如何使它更快的反馈，就像现在我的时间戳很低的时候，它需要很长时间才能运行。此外，我正在努力寻找一种实现行星的好方法，因为我现在拥有它们的方式不是这样的，我可以创建一个planet类，并初始化我将要使用的所有类，但是我宁愿只导入/使用它们。

下面是一些输出示例：

控制台：(旧输出)

Days remaining: 28.84, (96.14% left)
Tics per second: 49286 est time remaining: 00:00:50
Days remaining: 27.69, (92.28% left)
Tics per second: 49555 est time remaining: 00:00:48
Days remaining: 26.53, (88.43% left)
Tics per second: 49653 est time remaining: 00:00:46
Days remaining: 25.37, (84.57% left)
Tics per second: 49770 est time remaining: 00:00:44
Days remaining: 24.21, (80.71% left)
Tics per second: 49777 est time remaining: 00:00:42
Days remaining: 23.06, (76.85% left)
Tics per second: 49789 est time remaining: 00:00:40
Days remaining: 21.90, (72.99% left)
Tics per second: 49844 est time remaining: 00:00:37
Days remaining: 20.74, (69.14% left)
Tics per second: 49901 est time remaining: 00:00:35
Days remaining: 19.58, (65.28% left)
Tics per second: 49796 est time remaining: 00:00:33
Days remaining: 18.43, (61.42% left)
Tics per second: 49725 est time remaining: 00:00:32
Days remaining: 17.27, (57.56% left)
Tics per second: 49682 est time remaining: 00:00:30
Days remaining: 16.11, (53.70% left)
Tics per second: 49632 est time remaining: 00:00:28
Days remaining: 14.95, (49.85% left)
Tics per second: 49645 est time remaining: 00:00:26
Days remaining: 13.80, (45.99% left)
Tics per second: 49553 est time remaining: 00:00:24
Days remaining: 12.64, (42.13% left)
Tics per second: 49541 est time remaining: 00:00:22
Days remaining: 11.48, (38.27% left)
Tics per second: 49519 est time remaining: 00:00:20
Days remaining: 10.32, (34.41% left)
Tics per second: 49523 est time remaining: 00:00:18
Days remaining: 9.17, (30.56% left)
Tics per second: 49538 est time remaining: 00:00:15
Days remaining: 8.01, (26.70% left)
Tics per second: 49512 est time remaining: 00:00:13
Days remaining: 6.85, (22.84% left)
Tics per second: 49532 est time remaining: 00:00:11
Days remaining: 5.69, (18.98% left)
Tics per second: 49501 est time remaining: 00:00:09
Days remaining: 4.54, (15.12% left)
Tics per second: 49518 est time remaining: 00:00:07
Days remaining: 3.38, (11.27% left)
Tics per second: 49516 est time remaining: 00:00:05
Days remaining: 2.22, (7.41% left)
Tics per second: 49533 est time remaining: 00:00:03
Days remaining: 1.06, (3.55% left)
Tics per second: 49526 est time remaining: 00:00:01
Total Tics per second: 49509.11
Total time: 00:00:52

情节：

蓝线是飞行器(轨道被月球改变了)，绿线是卫星随时间变化的位置。

( 这里在一个git，请原谅我正在尝试开始读博士和狮身人面像，但不是很顺利)。

Orbital.py (主位)

import time as t
import math
import sys

from astropy.time import Time
import numpy as np

#Custum objects live in sub directory
sys.path.append('./Objects/')
# Import Objects
from Objects import Craft
import moon
import earth

# Generate output for plotting a sphere
def drawSphere(xCenter, yCenter, zCenter, r):
    # Draw sphere
    u, v = np.mgrid[0:2 * np.pi:20j, 0:np.pi:10j]
    x = np.cos(u) * np.sin(v)
    y = np.sin(u) * np.sin(v)
    z = np.cos(v)
    # Shift and scale sphere
    x = r * x + xCenter
    y = r * y + yCenter
    z = r * z + zCenter
    return (x, y, z)

def plot(ship,planets):
    """3d plots earth/moon/ship interaction"""
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    # Initilize plot
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.set_xlabel('X Km')
    ax.set_ylabel('Y Km')
    ax.set_zlabel('Z Km')
    ax.set_xlim3d(-500000, 500000)
    ax.set_ylim3d(-500000, 500000)
    ax.set_zlim3d(-500000, 500000)
    
    ax.plot(xs=ship.hist[0][0::10],
            ys=ship.hist[1][0::10],
            zs=ship.hist[2][0::10],
            zdir='z', label='ys=0, zdir=z')
    
    # Plot planet trajectory
    for planet in planets:
        ax.plot(xs=moon.hist[0],
                ys=moon.hist[1],
                zs=moon.hist[2],
                zdir='z', label='ys=0, zdir=z')
    
    # Plot Earth (plot is moon position relative to earth)
    # also plotting to scale
    (xs, ys, zs) = drawSphere(0, 0, 0, 6367.4447)
    ax.plot_wireframe(xs, ys, zs, color="r")
    
    plt.show()

def sim(startTime, endTime, step, ship, planets):
    """Runs orbital simulation given ship and planet objects as well as start/stop times"""

    # Caculate moon planet update rate (1/10th as often as the craft)
    plan_step = int(math.ceil(((endTime - startTime) / step)/100))

    # Initilize positions of planets
    for planet in planets:
                planet.getRelPos(startTime)
                planet.log()

    start = t.time()
    for i, time in enumerate(np.arange(startTime, endTime, step)):
        # Every plan_step update the position estmation
        if (i % plan_step == 0):
            for planet in planets:
                planet.getRelPos(time)
                planet.log()

        # Update craft_vol
        for planet in planets:
            ship.force_g(planet.mass, planet.pos[0], planet.pos[1], planet.pos[2])
        ship.update()

        # Log the position of the ship every 1000 steps
        if (i % 1000) == 0:
            # Append the position to the lists
            ship.log()

        # Print status update every 100,000 steps
        if (i % 100000) == 0:
            if t.time()-start > 1:
                print "Days remaining: {0:.2f}, ({1:.2f}% left)".format((endTime - time), ((1-((time - startTime) / (endTime - startTime)))*100))
                # Caculate estmated time remaining
                # Total tics
                tot_tics = int((endTime - startTime) / step)
                # Tics per second till now
                tics_s = int(math.ceil(i/(t.time()-start)))
                # Estmated remaining seconds
                sec_rem = (tot_tics-i)/tics_s
                m, s = divmod(sec_rem, 60)
                h, m = divmod(m, 60)
                print "Tics per second: {0} est time remaining: {1:0>2}:{2:0>2}:{3:0>2}".format(tics_s, h, m, s)
    print "Total Tics per second: {0:.2f}".format(((endTime - startTime) / step)/(t.time()-start))
    tot_s = t.time()-start
    m, s = divmod(tot_s, 60)
    h, m = divmod(m, 60)
    print "Total time: {0:0>2.0f}:{1:0>2.0f}:{2:0>2.0f}".format(h, m, s)

if __name__ == "__main__":

    # Delta time for simulations (in days)
    del_t = np.longdouble(1.0) / (24.0 * 60.0 * 60.0)

    ship = Craft(del_t, x=35786, y=1, z=1, v_x=0, v_y=4.5, v_z=0, mass=12)
    planets = [earth,moon]

    # Initilize simulation time
    start_time = Time('2015-09-10T00:00:00')
    end_time = Time('2015-10-10T00:00:00')

    # Initilize start date/time (Julian)
    time = start_time.jd
    
    sim(start_time.jd, end_time.jd, del_t, ship, planets)
    plot(ship,[moon])

./Orbital/Objects.py

import math
import numpy as np
from jplephem.spk import SPK
import os

class Craft(object):

    def __init__(self, delt_t, x=0.0, y=0.0, z=0.0,
                v_x=0.0, v_y=0.0, v_z=0.0, mass=0):
        # Pos is in km
        self.pos = np.array([np.longdouble(x),
                            np.longdouble(y),
                            np.longdouble(z)])
        # Vol is in km/s
        self.vol = np.array([np.longdouble(v_x),
                            np.longdouble(v_y),
                            np.longdouble(v_z)])
        # Mass is in kg
        self.mass = mass
        # Delta_t is in days
        self.del_t = np.longdouble(delt_t) * 86400

        # Initilize non input vars
        # Initilize history of positions
        self.hist = [[np.longdouble(x)],
                    [np.longdouble(y)],
                    [np.longdouble(z)]]
        # Initilize force array
        self.f = np.array([np.longdouble(0),
                        np.longdouble(0),
                        np.longdouble(0)])
        # Initilize constants
        # Gravitational constant (6.674*10-11 N*(m/kg)2)
        self.G_const = np.longdouble(0.00000000006674)

    def force_g(self, body_mass, body_x=0.0, body_y=0.0, body_z=0.0):
        # Caculate x, y, z distances
        xdist = (self.pos[0] - body_x)
        ydist = (self.pos[1] - body_y)
        zdist = (self.pos[2] - body_z)
        # Caculate vector distance
        d = math.sqrt((xdist**2) + (ydist**2) + (zdist**2)) * 1000
        # Caculate comman force of gravity
        g_com = ((self.G_const * body_mass * self.mass) / (d**3))

        # Update forces array
        self.f -= np.array([g_com * (xdist * 1000),
                            g_com * (ydist * 1000),
                            g_com * (zdist * 1000)])

    # Function to step simulation based on force caculations
    def update(self):
        # Step velocity
        self.vol += np.array([(((self.f[0] / self.mass) * self.del_t)) / 1000,
                            (((self.f[1] / self.mass) * self.del_t)) / 1000,
                            (((self.f[2] / self.mass) * self.del_t)) / 1000])
        # Step position
        self.pos += np.array([(self.del_t * self.vol[0]),
                            (self.del_t * self.vol[1]),
                            (self.del_t * self.vol[2])])
        # Reset force profile
        self.f = np.array([np.longdouble(0),
                        np.longdouble(0),
                        np.longdouble(0)])
    
    def VV_update(self):
        """Parameters:
        r is a numpy array giving the current position vector
        v is a numpy array giving the current velocity vector
        dt is a float value giving the length of the integration time step
        a is a function which takes x as a parameter and returns the acceleration vector as an array"""
        r_new = r + v*dt + a(r)*dt**2/2
        v_new = v + (a(r) + a(r_new))/2 * dt

    def dist(self, body_x=0.0, body_y=0.0, body_z=0.0):
        return math.sqrt(((self.pos[0] - body_x)**2) +
                        ((self.pos[1] - body_y)**2) +
                        ((self.pos[2] - body_z)**2))

    def log(self):
        self.hist[0].append(self.pos[0])
        self.hist[1].append(self.pos[1])
        self.hist[2].append(self.pos[2])

./Objects/earth.py

import numpy as np
from jplephem.spk import SPK
import os

if not os.path.isfile('de430.bsp'):
    raise ValueError('de430.bsp Was not found!')
kernel = SPK.open('de430.bsp')
# Initilize mass
mass = np.longdouble(5972198600000000000000000)
# Initilize history of positions
hist = [[],[],[]]
pos = np.array([np.longdouble(0),
                np.longdouble(0),
                np.longdouble(0)])

def getPos(time):
    """Returns the earths position relative to the solar system barycentere
    """
    return kernel[3, 399].compute(time)
def getRelPos(time):
    global pos
    """Returns relitive position of the earth (which is the origin
    in an earth moon system)"""
    pos = np.array([np.longdouble(0),
                    np.longdouble(0),
                    np.longdouble(0)])
    return pos
def log():
    global pos
    """log current position"""
    hist[0].append(pos[0])
    hist[1].append(pos[1])
    hist[2].append(pos[2])

./Objects/moon.py

import numpy as np
from jplephem.spk import SPK
import os

if not os.path.isfile('de430.bsp'):
    raise ValueError('de430.bsp Was not found!')
kernel = SPK.open('de430.bsp')
# Initilize mass
mass = np.longdouble(7.34767309 * 10**22)
# Initilize history of positions
hist = [[],[],[]]
pos = []

def getPos(time):
    """Returns the moons position relative to the solar system barycentere"""
    global pos
    pos = kernel[3, 301].compute(time)
    return pos
def getRelPos(time):
    """Returns relitive position of the moon (relative to the earth)"""
    global pos
    pos = kernel[3, 301].compute(time) - kernel[3, 399].compute(time)
    return np.array([np.longdouble(pos[0]),
                    np.longdouble(pos[1]),
                    np.longdouble(pos[2])])
def log():
    """log current position"""
    global pos
    global hist
    hist[0].append(pos[0])
    hist[1].append(pos[1])
    hist[2].append(pos[2])

(注意:如果您想运行这段代码，您将需要de430.bsp文件，您可以得到该这里)

我的更新代码可以找到这里。

Answer 1

下面是为运行给定次数而修改的代码，结果表明，只有大约3000次就足以获得类似的数字。你需要跑多少次它确实取决于不同物体的速度和转动速度(就像船吊索在月球或地球上飞行时一样)。

当然，当从250万步减少到3000步时，时间变得更好了.

planets.py和Objects.py文件与我的另一个答案一样。

文件：Orbital.py

import time as t
import math
import sys
from jplephem.spk import SPK
import os
from astropy.time import Time
import numpy as np

#Custum objects live in sub directory
sys.path.append('./Objects/')

# Import Objects
from Objects import Craft
from planets import Planet

import matplotlib.pyplot as plt
# The following is needed for the projection='3d' to work in plot()
from mpl_toolkits.mplot3d import Axes3D


def profile(func): return func

# Generate output for plotting a sphere
def drawSphere(xCenter, yCenter, zCenter, r):
    # Draw sphere
    u, v = np.mgrid[0:2 * np.pi:20j, 0:np.pi:10j]
    x = np.cos(u) * np.sin(v)
    y = np.sin(u) * np.sin(v)
    z = np.cos(v)
    # Shift and scale sphere
    x = r * x + xCenter
    y = r * y + yCenter
    z = r * z + zCenter
    return (x, y, z)

def plot(ship,planets):
    """3d plots earth/moon/ship interaction"""

    # Initialize plot
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.set_xlabel('X Km')
    ax.set_ylabel('Y Km')
    ax.set_zlabel('Z Km')
    ax.set_xlim3d(-500000, 500000)
    ax.set_ylim3d(-500000, 500000)
    ax.set_zlim3d(-500000, 500000)

    ax.plot(xs=ship.trajectory[0][0::10],
            ys=ship.trajectory[1][0::10],
            zs=ship.trajectory[2][0::10],
            zdir='z', label='ys=0, zdir=z')

    # Plot planet trajectory
    for planet in planets:
        ax.plot(xs=planet.trajectory[0],
                ys=planet.trajectory[1],
                zs=planet.trajectory[2],
                zdir='z', label='ys=0, zdir=z')

    # Plot Earth (plot is moon position relative to earth)
    # also plotting to scale
    (xs, ys, zs) = drawSphere(0, 0, 0, 6367.4447)
    ax.plot_wireframe(xs, ys, zs, color="r")

    plt.show()



@profile
def run_simulation(start_time, end_time, step, ship, planets):
    """Runs orbital simulation given ship and planet objects as well as start/stop times"""

    # Calculate moon & planet update rate (1/10th as often as the craft)
    planet_step_rate = int(math.ceil(((end_time - start_time) / step)/100))

    # Initialize positions of planets
    for planet in planets:
        planet.update_position(start_time)
        planet.log_position()

    start = t.time()
    total_time = end_time - start_time
    # Total tics
    total_steps = int((end_time - start_time) / step)

    print("Days in simulation: {:6.2f}, number of steps: {}".format(
            total_time, total_steps))    

    for i, time in enumerate(np.arange(start_time, end_time, step)):

        # Every planet_step_rate update the position estmation
        if (i % planet_step_rate == 0):
            for planet in planets[1:]:
                planet.update_position(time)
                planet.log_position()

        # Update craft velocity and force vectors related to planets
        for planet in planets:
            ship.force_g(planet)

        # Update ship position according to sum of affecting forces
        ship.update()

        # Log the position of the ship every 1000 steps
        if (i % 1000) == 0:
            # Append the position to the lists
            ship.log()

        # Print status update every 100,000 steps
        if (i % 100000) == 0:
            if t.time()-start > 1:
                time_since_start = time - start_time

                # Calculate estmated time remaining
                remaining_days_of_simulation = end_time - time
                remaining_percentage_of_simulation = (1 - time_since_start / total_time) * 100

                # Tics per second till now
                steps_pr_second = int(math.ceil( i / (t.time()-start)))
                # Estmated remaining seconds
                remaining_simulation_time_in_secs = (total_steps - i) / steps_pr_second

                # Split remaining_simulation_time_in_secs into hours, minutes and second
                minutes, seconds = divmod(remaining_simulation_time_in_secs, 60)
                hours, minutes = divmod(minutes, 60)
                print("  Simulation remaining: {:6.2f} days, {:5.2f}%.  Script statistics - steps/s: {:6}, est. time left: {:0>2}:{:0>2}:{:0>2}".format(
                      remaining_days_of_simulation, remaining_percentage_of_simulation,
                      steps_pr_second, hours, minutes, seconds))

    end = t.time()              
    elapsed = end - start
    print("Final steps per second: {0:.2f}".format((total_time / step)/elapsed))

    minutes, seconds = divmod(elapsed, 60)
    hours, minutes = divmod(minutes, 60)
    print("Total running time: {:0>2.0f}:{:0>2.0f}:{:0>2.0f}".format(
          hours, minutes, seconds))

def main():

    # Delta time for simulations (in days)
    delta_t = np.longdouble(1.0) / (24.0 * 60.0 * 60.0)

    ship = Craft(delta_t, x=35786, y=1, z=1, v_x=0, v_y=4.5, v_z=0, mass=12)

    # Initialize simulation time
    simulation_start = Time('2015-09-10T00:00:00')
    simulation_end = Time('2015-10-10T00:00:00')

    if not os.path.isfile('de430.bsp'):
        raise ValueError('de430.bsp Was not found!')
    kernel = SPK.open('de430.bsp')

    planets = [Planet(kernel, 399, 399, np.longdouble(5972198600000000000000000)),
               Planet(kernel, 301, 399, np.longdouble(7.34767309 * 10**22))]

    run_simulation(simulation_start.jd, simulation_end.jd, delta_t, ship, planets)
    plot(ship, [planets[1]])

if __name__ == "__main__":
    main()

这将在每一步更新行星和飞船，并打印历史上的所有位置。可以在船舶更新周围添加一个for循环，使其更新的频率稍微高一些，这几乎是没有好处的，因为船的计算是耗时的。

定时输出：

$ time python Orbital.py
Days in simulation:  30.00, number of steps: 5000
  Simulation remaining:  30.00 days, 100.00%.  Script statistics - remaining rounds: 5000
  Simulation remaining:  27.00 days,  90.00%.  Script statistics - steps/s:   2225, est. time left: 00:00:02
  Simulation remaining:  24.00 days,  80.00%.  Script statistics - steps/s:   2237, est. time left: 00:00:01
  Simulation remaining:  21.00 days,  70.00%.  Script statistics - steps/s:   2231, est. time left: 00:00:01
  Simulation remaining:  18.00 days,  60.00%.  Script statistics - steps/s:   2236, est. time left: 00:00:01
  Simulation remaining:  15.00 days,  50.00%.  Script statistics - steps/s:   2230, est. time left: 00:00:01
  Simulation remaining:  12.00 days,  40.00%.  Script statistics - steps/s:   2236, est. time left: 00:00:00
  Simulation remaining:   9.00 days,  30.00%.  Script statistics - steps/s:   2232, est. time left: 00:00:00
  Simulation remaining:   6.00 days,  20.00%.  Script statistics - steps/s:   2230, est. time left: 00:00:00
  Simulation remaining:   3.00 days,  10.00%.  Script statistics - steps/s:   2226, est. time left: 00:00:00
Final steps per second: 2225.10
Total running time: 00:00:02

real    0m6.939s
user    0m3.143s
sys 0m0.352s

这是7秒，包括编译Python和访问那个大型de430.bsp文件。而旋转到与OP旋转几乎相同的图像给出了这样的图像：

关于精度

的增编

在注释中，您说明您想要精确性，因此希望计算更多。在执行250万步和3000步之间有一种平衡，这与各种物体的速度和距离有关。距离越远，他们对彼此的影响就越小。

关于准确性的另一点是，您声称使用地球作为您的参考点，而您的代码实际上并没有这样做，以下是原因：

• 你的飞船开始于一个给定的点(35786，1，1)，这大概是地球在模拟启动时的位置。但你不能纠正船的位置相对于移动的地球位置
• 在计算月球位置时，你使用的是[3, 399]和[3, 301]，但你并不认为这是根据地球S重心对太阳系重心(又名0，3)的定位。参见https://pypi.python.org/pypi/jplephem/，关于“了解火星相对于地球的位置”的段落。
• 如果我没有弄错的话，地球的位置也是如此。

所以，如果你为了精确，你可能想重新评估数学，以及你的基准应该是什么，也就是如果你应该使用地球重心或者地球的实际位置，或者可能是太阳系的重心。在我提供的代码中，我假设你的数学是正确的，所以我没有抵消这些不准确之处。