Проблема с графиком земной орбиты с использованием python

Question

Я пытаюсь написать код для орбиты Земли в СИ с использованием интегратора симплекти c, моя попытка следующая:

import numpy as np
import matplotlib.pyplot as plt

#Set parameters
G = 6.67348e-11
mEar = 5.972e24
mSun = 1.989e30

def earth_orbit(x0, y0, vx0, vy0, N):
    dt = 1/N                    #timestep
    pos_arr = np.zeros((N,2))   #empty array to store position
    vel_arr = np.zeros((N,2))   #empty array to store velocities

    #Initial conditions
#     x0 = x
#     y0 = y
#     vx0 = vx
#     vy0 = vy
    pos_arr[0] = (x0,y0)        #set the intial positions in the array
    vel_arr[0] = (vx0,vy0)      #set the initial velocities in the array

    #Implement Verlet Algorithm
    for k in range (N-1):
        pos_arr[k+1] = pos_arr[k] + vel_arr[k]*dt                       #update positions
        force = -G * mSun * mEar * pos_arr[k+1] / (np.linalg.norm(pos_arr[k+1])**3)   #force calculation 
        vel_arr[k+1] = vel_arr[k] + (force/mEar) * dt                     #update velocities

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

    return pos_arr, vel_arr

earth_orbit(149.59787e9, 0, 0, 29800, 1000)

Результат - 2 точки, и я могу ' Не выяснить, проблема ли это в агрегате или в расчетах?

Answer 1

Отобразить траекторию

pos_arr содержит координаты x и y в своих столбцах. Таким образом, для отображения всей траектории можно использовать plt.plot(pos_arr[:,0], pos_arr[:,1]). Я бы предпочел использовать plt.plot(*pos_arr.T) как более короткую альтернативу. Строку, отображающую траекторию, необходимо заменить на:

plt.plot(*pos_arr.T, 'g', label = 'Earth trajectory')

Изменить временной шаг

Здесь временной шаг ( в секундах ) выбран как 1/N, где N - количество итераций. Итак, общая продолжительность симуляции равна timestep * N = 1 second! Для N=1000 вы можете вместо этого попробовать timestep = 3600*12 (полдня), чтобы общая продолжительность была чуть меньше 1,5 лет. Я предлагаю передать duration в качестве параметра функции earth_orbit, а затем установить timestep как duration / N.

def earth_orbit(x0, y0, vx0, vy0, N=1000, duration=3.15e7):
    dt = duration / N
    ...