地球轨道图 Python

Question

我正在用 Python 完成一个关于行星运动的项目，第一个任务是为地球绕太阳运行的轨道编写代码。这是我目前所拥有的：

def earth_orbit(rEar,v0):
    
    #r0 = xEar[0] = rEar 
    #v0 = vyEar[0] = np.sqrt(mu/rEar) 
    
    #Set parameters:
    N = 365                     # Earth days in a year
    dt = 1.00/N                 # Time Step: Fractions of a year - 1 Earth day (i.e. 1/365)
    mu = 4*np.pi**2             # Gravitational parameter 

    #Create an array, for all variables, of size N with all entries equal to zero:
    xEar = np.zeros((N,))
    yEar = np.zeros((N,))
    vxEar = np.zeros((N,))
    vyEar = np.zeros((N,))

    # Initial Conditions:
    xEar[0] = rEar                   # (x0 = r, y0 = 0) in AU
    vyEar[0] = v0                    #units: AU/yr

    #Implement Verlet Algorithm:
    for k in range(0,N-1):
        rEar = (xEar[k]**2+yEar[k]**2)**0.5
        vxEar[k+1] = vxEar[k] - (mu*xEar[k])/((rEar)**3)*dt
        xEar [k+1] = xEar[k] + vxEar[k+1]*dt
        vyEar[k+1] = vyEar[k] - (mu*yEar[k])/((rEar)**3)*dt
        yEar [k+1] = yEar[k] + vyEar[k+1]*dt

     #Plot:
    a = plt.plot(xEar, yEar, 'go', markersize = 1, label = 'Earth trajectory')
    plt.plot(0,0,'yo', label = 'Sun positon')        #yellow marker for the sun
    plt.plot(xEar[0],0,'bo', label = 'Earth initial positon')  #dark blue marker for earth's initial position
    plt.axis('equal')
    plt.xlabel ('x')
    plt.ylabel ('y')
    
    return a, xEar, yEar

这很好用，但是对于项目的后续部分，我被告知我不应该将 AU 用于此基本代码，因为它会造成困难。我尝试更改引力参数并将输入置于标准单位中，但图形变为仅 2 个点（xEar 和 vyEar 的值保持不变）而不是圆形轨道，我不知道为什么会这样?

Answer 1

由于引力常数的定义，当你改变长度单位时，你必须将mu的值改变为比例因子的立方。

为方便起见，我将mu 的定义拉到函数之外。

import numpy as np
import matplotlib.pyplot as plt

def earth_orbit(rEar, v0):

    #Set parameters:
    N = 365      # Earth days in a year
    dt = 1. / N  # Time Step: Fractions of a year - 1 Earth day (i.e. 1/365)

    #Create an array, for all variables, of size N with all entries equal to zero:
    xEar = np.zeros((N,))
    yEar = np.zeros((N,))
    vxEar = np.zeros((N,))
    vyEar = np.zeros((N,))

    # Initial Conditions:
    xEar[0] = rEar                   # (x0 = r, y0 = 0) in AU
    vyEar[0] = v0                    #units: AU/yr

    #Implement Verlet Algorithm:
    for k in range(0, N-1):
        rEar = (xEar[k]**2+yEar[k]**2)**0.5
        vxEar[k+1] = vxEar[k] - ((mu * xEar[k]) / (rEar**3)) * dt 
        xEar [k+1] = xEar[k] + vxEar[k+1]*dt
        vyEar[k+1] = vyEar[k] - ((mu * yEar[k]) / (rEar**3)) * dt 
        yEar [k+1] = yEar[k] + vyEar[k+1]*dt

     #Plot:
    a = plt.plot(xEar, yEar, 'go', markersize = 1, label = 'Earth trajectory')
    plt.plot(0,0,'yo', label = 'Sun positon')                  # yellow marker
    plt.plot(xEar[0],0,'bo', label = 'Earth initial positon')  # dark blue marker
    plt.axis('equal')
    plt.xlabel ('x')
    plt.ylabel ('y')

    return a, xEar, yEar


# average distance earth-sun (1 AU) in meter
au_to_m = 149_597_870_700.

mu = au_to_m**3 * 4 * np.pi**2  # Gravitational parameter 

earth_orbit(au_to_m, np.sqrt(mu / au_to_m));