如何在已经绘制的椭圆上绘制行星轨道作为时间的函数？

Question

我正在尝试为我的课程创建一个程序，用户在 16 年 1 月 1 日之后输入一个时间，该程序将根据需要多次运行行星的周期，然后绘制轨道图像行星在正确的位置。 我的程序不需要太复杂，所以我忽略了长时间的变化等。我只需要知道一种方法（可能使用积分）来映射椭圆上的一个点并考虑行星的变化速度帐户。我知道只要在同一时间间隔对两个点进行采样，由两个点和一个焦点绘制的区域保持不变，因此我对整合线方程的想法可能很有用，但我有点不合时宜深度。

这是我目前的代码，感谢您提供的任何帮助：

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
#Set axes aspect to equal as orbits are almost circular; hence square is needed
ax = plt.figure(0).add_subplot(111, aspect='equal')

#Setting the title, axis labels, axis values and introducing a grid underlay
#Variable used so title can indicate user inputed date
plt.title('Inner Planetary Orbits at[user input date]')
plt.ylabel('x10^6 km')
plt.xlabel('x10^6 km')
ax.set_xlim(-300, 300)
ax.set_ylim(-300, 300)
plt.grid()

#Creating the point to represent the sun at the origin (not to scale), 
ax.scatter(0,0,s=200,color='y')
plt.annotate('Sun', xy=(0,-30))

#Implementing ellipse equations to generate the values needed to plot an ellipse
#Using only the planet's min (m) and max (M) distances from the sun
#Equations return '2a' (the ellipses width) and '2b' (the ellipses height)
def OrbitLength(M, m):
    a=(M+m)/2
    c=a-m
    e=c/a
    b=a*(1-e**2)**0.5
    print(a)
    print(b)
    return 2*a, 2*b

#This function uses the returned 2a and 2b for the ellipse function's variables
#Also generating the orbit offset (putting the sun at a focal point) using M and m
def PlanetOrbit(Name, M, m):
    w, h = OrbitLength(M, m)
    Xoffset= ((M+m)/2)-m
    Name = Ellipse(xy=((Xoffset),0), width=w, height=h, angle=0, linewidth=1, fill=False)
    ax.add_artist(Name)

#These are the arguments taken from hyperphysics.phy-astr.gsu.edu/hbase/solar/soldata2.html
#They are the planet names, max and min distances, and their longitudinal angle
#Also included is Halley's Comet, used to show different scale  and eccentricity
PlanetOrbit('Mercury', 69.8, 46.0)
PlanetOrbit('Venus', 108.9, 107.5)
PlanetOrbit('Earth', 152.1, 147.1)
PlanetOrbit('Mars', 249.1, 206.7)
PlanetOrbit("Halley's Comet",45900,88)

Answer 1

在this gamedev question 找到了类似的东西。那里有人链接this wikipedia article，它描述了一种逐步计算位置（在日心极坐标中）作为时间函数的方法。这些方程当然包含牛顿的引力常数和太阳的质量。其中一些只能通过数值方法解决。

上述方法的实现：

from math import *

EPSILON = 1e-12
def solve_bisection(fn, xmin,xmax,epsilon = EPSILON):
  while True:
      xmid = (xmin + xmax) * 0.5
      if (xmax-xmin < epsilon):
        return xmid
      fn_mid = fn(xmid)
      fn_min = fn(xmin)
      if fn_min*fn_mid < 0:
          xmax = xmid
      else:
          xmin = xmid
    
'''
Found something similar at this gamedev question:
https://gamedev.stackexchange.com/questions/11116/kepler-orbit-get-position-on-the-orbit-over-time?newreg=e895c2a71651407d8e18915c38024d50

Equations taken from:
https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion#Position_as_a_function_of_time
'''
def SolveOrbit(rmax, rmin, t):
  # calculation precision
  epsilon = EPSILON
  # mass of the sun [kg]
  Msun = 1.9891e30
  # Newton's gravitational constant [N*m**2/kg**2]
  G = 6.6740831e-11
  # standard gravitational parameter
  mu = G*Msun
  # eccentricity
  eps = (rmax - rmin) / (rmax + rmin)
  # semi-latus rectum
  p = rmin * (1 + eps)
  # semi/half major axis
  a = p / (1 - eps**2)
  # period
  P = sqrt(a**3 / mu)
  # mean anomaly
  M = (t / P) % (2*pi)
  # eccentric anomaly
  def fn_E(E):
    return M - (E-eps*sin(E))
  E = solve_bisection(fn_E, 0, 2*pi)
  # true anomaly
  # TODO: what if E == pi?
  theta = 2*atan(sqrt((((1+eps)*tan(E/2)**2)/(1-eps))))
  # if we are at the second half of the orbit
  if (E > pi):
    theta = 2*pi - theta
  # heliocentric distance
  r = a * (1 - eps * cos(E))
  return theta, r

def DrawPlanet(name, rmax, rmin, t):
  SCALE = 1e9
  theta, r = SolveOrbit(rmax * SCALE, rmin * SCALE, t)
  x = -r * cos(theta) / SCALE
  y = r * sin(theta) / SCALE
  planet = Circle((x, y), 8)
  ax.add_artist(planet)

所以整个程序看起来像这样：

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse, Circle

#Set axes aspect to equal as orbits are almost circular; hence square is needed
ax = plt.figure(0).add_subplot(111, aspect='equal')

#Setting the title, axis labels, axis values and introducing a grid underlay
#Variable used so title can indicate user inputed date
plt.title('Inner Planetary Orbits at[user input date]')
plt.ylabel('x10^6 km')
plt.xlabel('x10^6 km')
ax.set_xlim(-300, 300)
ax.set_ylim(-300, 300)
plt.grid()

#Creating the point to represent the sun at the origin (not to scale), 
ax.scatter(0,0,s=200,color='y')
plt.annotate('Sun', xy=(0,-30))

#Implementing ellipse equations to generate the values needed to plot an ellipse
#Using only the planet's min (m) and max (M) distances from the sun
#Equations return '2a' (the ellipses width) and '2b' (the ellipses height)
def OrbitLength(M, m):
    a=(M+m)/2
    c=a-m
    e=c/a
    b=a*(1-e**2)**0.5
    print(a)
    print(b)
    return 2*a, 2*b

#This function uses the returned 2a and 2b for the ellipse function's variables
#Also generating the orbit offset (putting the sun at a focal point) using M and m
def PlanetOrbit(Name, M, m):
    w, h = OrbitLength(M, m)
    Xoffset= ((M+m)/2)-m
    Name = Ellipse(xy=((Xoffset),0), width=w, height=h, angle=0, linewidth=1, fill=False)
    ax.add_artist(Name)

    
    
from math import *

EPSILON = 1e-12
def solve_bisection(fn, xmin,xmax,epsilon = EPSILON):
  while True:
      xmid = (xmin + xmax) * 0.5
      if (xmax-xmin < epsilon):
        return xmid
      fn_mid = fn(xmid)
      fn_min = fn(xmin)
      if fn_min*fn_mid < 0:
          xmax = xmid
      else:
          xmin = xmid
    
'''
Found something similar at this gamedev question:
https://gamedev.stackexchange.com/questions/11116/kepler-orbit-get-position-on-the-orbit-over-time?newreg=e895c2a71651407d8e18915c38024d50

Equations taken from:
https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion#Position_as_a_function_of_time
'''
def SolveOrbit(rmax, rmin, t):
  # calculation precision
  epsilon = EPSILON
  # mass of the sun [kg]
  Msun = 1.9891e30
  # Newton's gravitational constant [N*m**2/kg**2]
  G = 6.6740831e-11
  # standard gravitational parameter
  mu = G*Msun
  # eccentricity
  eps = (rmax - rmin) / (rmax + rmin)
  # semi-latus rectum
  p = rmin * (1 + eps)
  # semi/half major axis
  a = p / (1 - eps**2)
  # period
  P = sqrt(a**3 / mu)
  # mean anomaly
  M = (t / P) % (2*pi)
  # eccentric anomaly
  def fn_E(E):
    return M - (E-eps*sin(E))
  E = solve_bisection(fn_E, 0, 2*pi)
  # true anomaly
  # TODO: what if E == pi?
  theta = 2*atan(sqrt((((1+eps)*tan(E/2)**2)/(1-eps))))
  # if we are at the second half of the orbit
  if (E > pi):
    theta = 2*pi - theta
  # heliocentric distance
  r = a * (1 - eps * cos(E))
  return theta, r

def DrawPlanet(name, rmax, rmin, t):
  SCALE = 1e9
  theta, r = SolveOrbit(rmax * SCALE, rmin * SCALE, t)
  x = -r * cos(theta) / SCALE
  y = r * sin(theta) / SCALE
  planet = Circle((x, y), 8)
  ax.add_artist(planet)
    
#These are the arguments taken from hyperphysics.phy-astr.gsu.edu/hbase/solar/soldata2.html
#They are the planet names, max and min distances, and their longitudinal angle
#Also included is Halley's Comet, used to show different scale  and eccentricity
PlanetOrbit('Mercury', 69.8, 46.0)
PlanetOrbit('Venus', 108.9, 107.5)
PlanetOrbit('Earth', 152.1, 147.1)
PlanetOrbit('Mars', 249.1, 206.7)
PlanetOrbit("Halley's Comet",45900,88)
for i in range(0, 52):
  DrawPlanet('Earth', 152.1, 147.1, i/52 * 365.25 *60*60*24)
for i in range(-2, 3):
  DrawPlanet("Halley's Comet",45900,88, 7*i *60*60*24)
print(60*60*24*365)
  
plt.show()

编辑：感谢@kdubs，sn-p 现在正在使用封闭形式的 theta 解决方案。