python - 使用极坐标法追踪相位、轴和地球方向的椭圆轨道

Question

我想表示两颗恒星组成的双星系统的椭圆轨道。我的目标是这样的:

我有一个沿轴尺寸的网格，焦点处有一个按比例的恒星，以及第二颗恒星的轨道。沿轨道的十进制数字是轨道相位。底部的箭头是地球方向，轨道上较粗的部分与特定情况的观测有关 - 我不需要它。我想从这个情节中改变的是:

1. 轨道相位:我想要从焦点到轨道的“虚线”，以及它们上方的轨道相位，而不是沿轨道的数字:

2. 我不想沿着 (0, 0) 画十字；

3. 我想重新调整轨道方向，以便 0.0 相位位于图的左上方，地球方向是向上的直箭头(我的系统的参数不同来自此处绘制的)。

我尝试寻找 python 示例，但我得到的唯一结果( from here )是极坐标图:

这并不能真正代表我想要的，但仍然是一个开始:

import numpy as np
import matplotlib.pyplot as plt
cos = np.cos
pi = np.pi

a = 10
e = 0.1
theta = np.linspace(0,2*pi, 360)
r = (a*(1-e**2))/(1+e*cos(theta))

fig = plt.figure()
ax = fig.add_subplot(111, polar=True)
ax.set_yticklabels([])
ax.plot(theta,r)

print(np.c_[r,theta])
plt.show()

Answer 1

这里有一些东西可以让你非常接近。您不需要极坐标来绘制像样的椭圆。您可以轻松利用所谓的艺术家。
您可能必须自定义轴标签，并且如果需要，还可以插入一两个箭头:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

# initializing the figure:
fig = plt.figure()

# the (carthesian) axis:
ax = fig.add_subplot(111,aspect='equal')
ax.grid(True)

# parameters of the ellipse:
a = 5.0
e = 4.0
b = np.sqrt(a**2.0 - e**2.0)

# the center of the ellipse:
x = 6.0
y = 6.0

# the angle by which the ellipse is rotated:
angle = -45.0
#angle = 0.0

# plotting the ellipse, using an artist:
ax.add_artist(Ellipse(xy=[x,y], width=2.0*a, height=2.0*b, \
                                angle=angle, facecolor='none'))
ax.set_xlim(0,2.0*x)
ax.set_ylim(0,2.0*y)

# marking the focus (actually, both)
# and accounting for the rotation of the ellipse by angle
xf = [x - e*np.cos(angle * np.pi/180.0),
      x + e*np.cos(angle * np.pi/180.0)]

yf = [y - e*np.sin(angle * np.pi/180.0),
      y + e*np.sin(angle * np.pi/180.0)]

ax.plot(xf,yf,'xr')

# plotting lines from the focus to the ellipse:
# these should be your "rays"
t = np.arange(np.pi,3.0*np.pi,np.pi/5.0)
p = b**2.0 / a
E = e / a
r = [p/(1-E*np.cos(ti)) for ti in t]

# converting the radius based on the focus 
# into x,y coordinates on the ellipse:
xr = [ri*np.cos(ti) for ri,ti in zip(r,t)]
yr = [ri*np.sin(ti) for ri,ti in zip(r,t)]

# accounting for the rotation by anlge:
xrp = [xi*np.cos(angle * np.pi/180.0) - \
       yi*np.sin(angle * np.pi/180.0) for xi,yi in zip(xr,yr)]
yrp = [xi*np.sin(angle * np.pi/180.0) + \
       yi*np.cos(angle * np.pi/180.0) for xi,yi in zip(xr,yr)]

for q in range(0,len(t)):
    ax.plot([xf[0], xf[0]+xrp[q]],[yf[0], yf[0]+yrp[q]],'--b')

# put labels outside the "rays"
offset = 0.75
rLabel = [ri+offset for ri in r]
xrl = [ri*np.cos(ti) for ri,ti in zip(rLabel,t)]
yrl = [ri*np.sin(ti) for ri,ti in zip(rLabel,t)]

xrpl = [xi*np.cos(angle * np.pi/180.0) - \
        yi*np.sin(angle * np.pi/180.0) for xi,yi in zip(xrl,yrl)]
yrpl = [xi*np.sin(angle * np.pi/180.0) + \
        yi*np.cos(angle * np.pi/180.0) for xi,yi in zip(xrl,yrl)]

# for fancy label rotation reduce the range of the angle t:
tlabel = [(ti -np.pi)*180.0/np.pi for ti in t]
for q in range(0,len(tlabel)):
    if tlabel[q] >= 180.0:
        tlabel[q] -= 180.0

# convert the angle t from radians into degrees:
tl = [(ti-np.pi)*180.0/np.pi for ti in t]

for q in range(0,len(t)):
    rotate_label = angle + tlabel[q]
    label_text = '%.1f' % tl[q]
    ax.text(xf[0]+xrpl[q],yf[0]+yrpl[q],label_text,\
            va='center', ha='center',rotation=rotate_label)

plt.show()

上面的例子将产生下图:

说明:

• 您可以使用 artist to plot the ellipse ，而不是使用极坐标
• 命名法基于 Wikipedia 上的定义
• 艺术家设置中的角度可旋转椭圆。该角度稍后用于旋转光线和标签的坐标(这只是数学)
• 光线源自polar form of the ellipse relative to a focus 。
• 角度 t 从 pi 到 3.0*pi 因为我假设这对应于您的光线应该从哪里开始的想法。对于 0 到 2.0*pi，您会得到相同的光线。我用过np.arange而不是 linspace，因为在此示例中我想要一个定义的增量(pi/5.0，或 36 度)。
• 光线末端的标签放置为 text ，变量 offset 控制椭圆和标签之间的距离。根据需要调整此项。
• 为了使标签文本方向与光线对齐，我将角度 t 减小到 0 到 180 度范围。与 0 到 360 度的整个范围相比，这具有更好的可读性。
• 为了简单起见，对于标签文本，我仅使用角度 t。将其替换为更适合您目的的任何信息。
• 在放置标签的循环之前，角度 t 从弧度转换为度数。在循环内，tl 的每个元素都会转换为字符串。这允许更多的格式控制(例如，如果您需要 3 位小数，则为 %.3f)