地球轨道-跨越式Python

Question

我正在做这个家庭作业练习：地球轨道

我的阴谋没有显示出整个轨迹。我不知道这是我的方程式出了什么问题，或者是一个密谋的问题。

干杯!

地球轨道图

import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate as spi

G = 6.6738*10**-11
M = 1.9891*10**30
h = 3600
y = 1.4710*10**11
vx = 3.0287*20**4

def LeapFrog(f, t_start, t_stop, z0, h):

    t_vec = np.arange(t_start, t_stop, h) 
    n = len(t_vec)
    d = len(z0)  
    z_vec = np.zeros((n, d))
    z_vec[0,:] = z0
    z_half_step=z_vec[0 , :] + (1/2)*h*f(z0,t_vec[0]) 
    
    
    for i in range(0, n - 1):
        z_vec[i+1,:]=z_vec[i,:] + h*f(z_half_step, t_vec[i])
        z_half_step += h*f(z_vec[i+1,:], t_vec[i])   

    return t_vec, z_vec, 


def f(z,t):   
    
    x=z[0]
    y=z[1] 
    vx=z[2] 
    vy=z[3] 
    r=np.sqrt(x**2+y**2)

    dz=np.zeros(4)
    
    dz[0]=vx 
    dz[1]=vy
    dz[2]=-G*M*x/r**3
    dz[3]=-G*M*y/r**3

    return dz

t_start = 0
t_stop = 24*365*5
z0 = np.array([0,y,vx,0])
t_vec, z_vec = LeapFrog(f, t_start, t_stop, z0, h)

plt.plot(z_vec[:,0],z_vec[:,1], 'g', markersize=1, label='Earth trajectory')
plt.plot(0,0,'yo', label = 'Sun positon')  
plt.show()

Answer 1

好吧，有个很有趣的错误：

vx = 3.0287*20**4

这比你最想写的要大16倍，只比太阳逃逸速度少20%。

t_stop也是以小时为单位编写的，但其余代码假设是正常的SI秒(可以通过使用t_stop = 24 * 365 * 0.99 * h进行验证，这将导致几乎一个完整的轨道)。它需要乘以h。Sidenote:我喜欢使用astropy.units跟踪物理量！

不过，代码本身工作得很好！