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

2021-04-11 14:54:30 发布

您现在位置:Python中文网/ 问答频道 /正文

我想代表双星双星系统的椭圆轨道。我的目标是这样的:

enter image description here

在这里,我有一个沿着轴线大小的网格,焦点处有一个按比例排列的恒星,还有次恒星的轨道。轨道上的十进制数是轨道相位。底部的箭头是地球的方向,轨道上较厚的部分与特定情况下的观测有关——我不需要它。我想从这个情节中改变的是:

  1. 轨道相位:我想要的是从焦点到轨道的“虚线”,而不是轨道上的数字:

  2. 我不想让十字架沿着(0,0);

  3. 我想重新定位轨道,以便0.0相位在图的左上角,地球方向是一个向上的直箭头(我的系统参数与这里绘制的不同)。

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

enter image description 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()
1条回答
网友
1楼 ·

这里有一些东西可以让你离得很近。你不需要极坐标来绘制一个像样的椭圆。有一种所谓的artist你可以很容易地利用。
您可能需要自定义轴标签,如果需要,可以插入一两个箭头:

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()

上面的示例将生成此图:

example

说明:

  • 您可以使用^{} to plot the ellipse,而不是使用极坐标
  • 该术语基于Wikipedia上可用的定义
  • 艺术家设置中的angle旋转椭圆。该角度稍后用于旋转光线和标签的坐标(这是只是数学)
  • 射线来自polar form of the ellipse relative to a focus。在
  • 角度从pi到{},因为我假设这与你关于光线应该从哪里开始的想法相对应。对02.0*pi得到相同的光线。我使用^{}而不是linspace,因为我希望在这个例子中定义一个增量(pi/5.0,或36度)。在
  • 射线末端的标签被放置为^{},变量offset控制椭圆和标签之间的距离。根据需要调整。在
  • 为了使标签文本方向与光线对齐,我将角度t减小到0180度。与0360度的整个范围相比,这使得可读性更好。在
  • 对于标签文本,为了简单起见,我使用了角度t。用更适合你目的的信息来代替它。在
  • 角度t,在放置标签的循环之前从弧度转换为度。在循环中,tl的每个元素都被转换为一个字符串。这允许更多的格式控制(例如%.3f如果需要3个小数)

相关问题