如何在已绘制的椭圆上绘制作为时间函数的行星轨道？
How do I plot a planet's orbit as a function of time on an already plotted ellipse?
我正在尝试为我的课程创建一个程序，用户在 2016 年 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)
在 this gamedev question. Someone there linked 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，该代码段现在对 theta 使用封闭形式的解决方案。
我正在尝试为我的课程创建一个程序，用户在 2016 年 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)
在 this gamedev question. Someone there linked 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，该代码段现在对 theta 使用封闭形式的解决方案。