利用Numpy,Matplotlib，axartist可视化地球轨道

（本章不讨论物理问题 只有程序的思路）

上面为程序设计思路 可直接跳转至下面完整代码

为了简单的模拟地球的受力情况 我们假设太阳系中只有地球和太阳两个物体

于是乎 我们需要计算的变量以及记录的常量只有 星球半径 星球质量 星球间距离 引力常量等

为了更好的计算这些数据 我们使用一个字典 其中存放着所有参量的列表

以及为了实现以太阳为中心的二维俯视图 我们使用求解极坐标的方式画出图像

并且将极坐标划分为四个象限 根据不同象限改变角度的大小

在后面调试的过程中 我们使用openpyxl这个包实现了和excel的交互

每次运行完的数据都会被自动存储到一个名为“resultsText的xlsx中 绝对路径需要自己修改

代码如下

import random, numpy as np, pandas as pd
from IPython.core.display import display
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import openpyxl
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation
import time

resultsFilePath = r"C:\\Users\\17006\\Desktop\\CSDN\\CSDN PlanetaryOrbital\\resultsText.txt"
resultsTextName = "resultsText.xlsx"
#     Basic variables:
orbitalColor = "black"
pointColor = "black"
periodPointColor = "red"
ax = plt.subplot(111, projection="polar")  # Draw the figure:
#    initial values:
G = 6.67e-11
t_perTime = 60000
t_realtime = 10
sun_mass = 1.989e+30
interval_limit = 1.1
specialPointSize = 15
textSize = 10
AU_M_constant = 1.49597871e+11
pauseTime = 0.001
pi_Value = math.pi
omiga = 7.292e-5
#    ListDatabase:
dataBase = {"SxList": [0.9832687 * AU_M_constant],
            "SyList": [0.0], "rList": [0.9832687 * AU_M_constant],
            "r0_3List": [1 / ((0.9832687 * AU_M_constant) ** 3)],
            "VxList": [0], "VyList": [30750],
            "axList": [-0.0061314726564], "ayList": [0], "tList": [0]}


# defined functions:
def translation_AU_to_Meter(AU):
    return AU * 1.49597871e+11


def roundFunction(x):
    return [round(a, 3) for a in x]


def getAngle(Ylist, Xlist):
    return np.arctan(Ylist / Xlist)


def getLength(Ylist, Xlist):
    return (Ylist ** 2 + Xlist ** 2) ** 0.5


def drawDynamicOrbital(i):
    ax.scatter(changeAngle(angleList)[i], lengthList[i], c=pointColor, s=pointSize)
    plt.pause(0.00001 / dataTimes)


def changeAngle(angleList):
    transAngleList = [0]
    index = 0
    for angle in angleList:
        angle = abs(angle)
        if dataBase["SxList"][index] < 0 and dataBase["SyList"][index] > 0:
            transAngleList.append(pi_Value - angle)
        if dataBase["SxList"][index] < 0 and dataBase["SyList"][index] < 0:
            transAngleList.append(pi_Value + angle)
        if dataBase["SxList"][index] > 0 and dataBase["SyList"][index] < 0:
            transAngleList.append(2 * (pi_Value) - angle)
        if dataBase["SxList"][index] > 0 and dataBase["SyList"][index] > 0:
            transAngleList.append(angle)
        index += 1
    return transAngleList


while True:
    dataCounts = input("Enter dataCounts:")
    try:
        dataTimes = int(dataCounts)
        pointSize = 1500 / (dataTimes)
        break
    except:
        print("Try again,enter dataCounts:")  # get the times of data should be calculated

for i in range(0, dataTimes):
    dataBase["SxList"].append((dataBase["SxList"][i] + t_perTime * (dataBase["VxList"][i])))
    dataBase["SyList"].append((dataBase["SyList"][i] + t_perTime * (dataBase["VyList"][i])))
    dataBase["rList"].append((math.sqrt(dataBase["SxList"][i + 1] ** 2 + dataBase["SyList"][i + 1] ** 2)))
    dataBase["r0_3List"].append((dataBase["rList"][i] ** 3) ** (-1))

    dataBase["axList"].append(-G * sun_mass * dataBase["SxList"][i] / (dataBase["rList"][i] ** 3))
    dataBase["ayList"].append(-G * sun_mass * dataBase["SyList"][i] / (dataBase["rList"][i] ** 3))

    dataBase["VxList"].append((dataBase["VxList"][i]) + t_perTime * dataBase["axList"][i + 1])
    dataBase["VyList"].append((dataBase["VyList"][i]) + t_perTime * dataBase["ayList"][i + 1])

    dataBase["tList"].append(i * t_perTime)  # Calculate and store all the data

SxListTranslated = np.array(dataBase["SxList"])
SyListTranslated = np.array(dataBase["SyList"])
angleList = getAngle(SyListTranslated, SxListTranslated)
lengthList = (SxListTranslated ** 2 + SyListTranslated ** 2) ** 0.5

df = pd.DataFrame(dataBase)
df.to_excel(resultsTextName)  # Save the results

##ani = FuncAnimation(fig,drawDynamicOrbital,frames = len(angleList),interval = 500,blit=False, repeat=False)#frame is a Iterative integer that would be passed to func which is drawDynamicOrbital
ax.scatter(angleList[0], lengthList[0], c=periodPointColor, s=specialPointSize)  # draw the original point
ax.text(angleList[0] - 0.1, lengthList[0] + 0.01, 'Earth', size=textSize, color="red")
ax.scatter(0, 0, c=periodPointColor, s=specialPointSize + 2)  # draw the sun
ax.text(0.2, 0.2, 'Sun', size=textSize, color="red")
ax.scatter(changeAngle(angleList), lengthList, c=pointColor, s=pointSize)
# draw the dynamic graph
##for orbitalPeriod in range(len(angleList)):
##    drawDynamicOrbital(orbitalPeriod)
plt.show()
# Stage 2 在使用了刚开始的模拟数据之后 我们发现这颗行星的运动轨迹会偏离原始运动轨道 我们开始怀疑 是不是因为初始值没有符合真实情况的原因 模拟数据可以拟合之后 我们带入了水星单独一个行星的数据 发现轨道无法回归 轨道会越来越偏移原始位置 因此 我们引入了时间作为第三个维度 画出了三维的图像

图为1000个60,000s之后的结果