Построение эллиптических орбит

Я пытаюсь написать код, который отображает эллиптические пути объекта, используя уравнение для эллипса r=a(1-e^2)/(1+e*cos(theta)). Я также хотел бы, чтобы эти данные были помещены в массив для другого использования.

from numpy import *#Imports Python mathematical functions library
import matplotlib.pyplot as plt #Imports plot library
from pylab import *

a = 5
e = 0.3
theta = 0
while theta <= 2*pi:
    r = (a*(1-e**2))/(1+e*cos(theta))
    print("r = ",r,"theta = ",theta)
    plt.polar(theta, r)
    theta += pi/180

plt.show()

Код выдает правильные значения для r и тета, но график пустой. Появляется окно полярного графика, но ничего не отображается.

Пожалуйста помоги. Заранее спасибо.


python matplotlib astronomy orbital-mechanics
person maelstromscientist    schedule 11.09.2012    source источник
comment
Один быстро очевидный недостаток заключается в том, что тета увеличивается на 180 градусов (в радианах) — разве вам не нужны меньшие шаги, скажем, 1 градус?   -  person DarenW    schedule 13.10.2012

Ответы (2)


arrow_upward
9
arrow_downward

Не вызывайте plt.polar один раз для каждой точки. Вместо этого вызовите его один раз со всеми данными в качестве входных данных:

import numpy as np #Imports Python mathematical functions library
import matplotlib.pyplot as plt #Imports plot library
cos = np.cos
pi = np.pi

a = 5
e = 0.3
theta = np.linspace(0,2*pi, 360)
r = (a*(1-e**2))/(1+e*cos(theta))
plt.polar(theta, r)

print(np.c_[r,theta])

plt.show()

введите здесь описание изображения

Кстати, numpy может выполнять расчет как двухстрочный, вместо использования цикла while:

theta = np.linspace(0,2*pi, 360)   # 360 equally spaced values between 0 and 2*pi
r = (a*(1-e**2))/(1+e*cos(theta))

Это определяет theta и r как массивы numpy (а не отдельные значения).

person unutbu    schedule 11.09.2012
comment
+1 тот же ответ, что и я :) в основном ... плюс лучшая информация о linspace - person Joran Beasley; 11.09.2012

arrow_upward
2
arrow_downward

Я думаю, вам нужно сделать points.append([theta,r]), а затем в конце plt.polar(points) ... это тоже делает довольно аккуратный дизайн

from numpy import *#Imports Python mathematical functions library
import matplotlib.pyplot as plt #Imports plot library
from pylab import *

a = 5
e = 0.3
theta = 0

points = []
while theta <= 2*pi:
    r = (a*(1-e**2))/(1+e*cos(theta))
    print("r = ",r,"theta = ",theta)
    points.append((theta, r))
    theta += pi/180
#plt.polar(points) #this is cool but probably not what you want
plt.polar(*zip(*points))
plt.show()
person Joran Beasley    schedule 11.09.2012
comment
Я надеялся, что это Юпитер, но эксцентриситет намного слишком велик для этого... - person mgilson; 11.09.2012
comment
Наши ответы в основном одинаковы, но наши сюжеты совершенно разные! :) Я немного в недоумении, как исправить (твое или мое?), но поскольку тета идет от 0 до 2pi, разве график не должен делать одну орбиту? - person unutbu; 11.09.2012
comment
Ах, наши графики были бы такими же, если бы вы использовали plt.polar(*zip(*points)), чтобы plt.polar получал два аргумента вместо одного списка. Я не уверен, что делает plt.polar, когда ему дается список кортежей... - person unutbu; 11.09.2012
comment
лол, о, я думал, что это был просто сюжет ... я думал, что это выглядит опрятно - person Joran Beasley; 11.09.2012