Py-ser
Py-ser

Reputation: 2098

Elliptic orbit with polar method tracing phases, axes and Earth direction

I would like to represent the elliptical orbit of a binary system of two stars. What I aim to, is something like this:

enter image description here

Where I have a grid of the sizes along the axes, an in-scale star at the focus, and the orbit of the secondary star. The decimal numbers along the orbit are the orbital phases. The arrow at the bottom is the Earth direction, and the thick part at the orbit is related to the observation for that specific case - I don't need it. What I want to change from this plot is:

  1. Orbital phase: instead of numbers along the orbit, I would like "dashed rays" from the focus to the orbit, and the orbital phase above them:

  2. I don't want the cross along (0, 0);

  3. I would like to re-orient the orbit, in order that the 0.0 phase is around the top left part of the plot, and the Earth direction is an upward pointing straight arrow (the parameters of my system are different from the one plotted here).

I tried to look for python examples, but the only thing I came out with (from here), is a polar plot:

enter image description here

which is not really representative of what I want, but still is a beginning:

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

Upvotes: 1

Views: 1746

Answers (1)

Schorsch
Schorsch

Reputation: 7905

Here's something that gets you very close. You do not need polar coordinates to plot a decent ellipse. There is a so-called artist you can readily utilize.
You probably have to customize the axis labels and maybe insert an arrow or two if you want:

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

The example above will result in this figure:

example

Explanations:

  • You can use an artist to plot the ellipse, instead of using polar coordinates
  • The nomenclature is based on the definitions available on Wikipedia
  • The angle in the artist setup rotates the ellipse. This angle is later used to rotate coordinates for the rays and labels (this is just math)
  • The rays are derived from the polar form of the ellipse relative to a focus.
  • The angle t runs from pi to 3.0*pi because I assumed that this would correspond to your idea of where the rays should start. You get the same rays for 0 to 2.0*pi. I used np.arange instead of linspace because I wanted a defined increment (pi/5.0, or 36 degrees) in this example.
  • The labels at the end of the rays are placed as text, the variable offset controls the distance between the ellipse and the labels. Adjust this as needed.
  • For the alignment of the label text orientation with the rays, I reduced the angle, t, to the range 0 to 180 degrees. This makes for better readability compared to the full range of 0 to 360 degrees.
  • For the label text I just used the angle, t, for simplicity. Replace this with whatever information better suits your purpose.
  • The angle, t, was converted from radians to degrees before the loop that places the label. Inside the loop, each element of tl is converted to a string. This allows for more formatting control (e.g. %.3f if you needed 3 decimals)

Upvotes: 3

Related Questions