2016-04-13 48 views
4

我一直在編寫一個Python腳本(GitHub LINK),用於可視化小行星/彗星/流星體軌道。劇本還繪製了行星及其軌道的位置。Matplotlib 3D plot - 從某些角度看「參數化曲線」

它適用於小半長軸(即「較小」軌道)的軌道。但是當我有一個超越海王星的軌道時(例如哈雷式彗星),並且從某些角度來看,有一個奇怪的「環繞」(因爲沒有更好的詞)效應。

讓我告訴你我是什麼意思:

圖像製作:http://i.imgur.com/onSZG8s.png enter image description here

  1. 此圖顯示了從一個角度它不會破壞陰謀。

  2. 當您向右旋轉相同的繪圖時,就好像該軌道摺疊了一半並且顛倒了它的方向!

  3. 如果你從遠距離看情節,你可以看到elipse被繪製爲應該的樣子。

這裏是可用於複製問題的最小版本的代碼。只有當攝像機的視角與大軌道平行時,纔會出現「環繞」。

from mpl_toolkits.mplot3d import Axes3D 
import numpy as np 
import matplotlib.pyplot as plt 

def orbitalElements2Cartesian(a, e, I, peri, node, E): 
    """ Convert orbital elements to Cartesian coordinates in the Solar System. 

    Args: 
     a (float): semi-major axis (AU) 
     e (float): eccentricity 
     I (float): inclination (degrees) 
     peri (float): longitude of perihelion (degrees) 
     node (float): longitude of ascending node (degrees) 
     E (float): eccentric anomaly (radians) 

    """ 

    # The source of equations used: 
    # http://farside.ph.utexas.edu/teaching/celestial/Celestialhtml/node34.html 

    # Check if the orbit is parabolic or hyperbolic 
    if e >=1: 
     e = 0.99999999 

    # Convert degrees to radians 
    I, peri, node = map(np.radians, [I, peri, node]) 

    # True anomaly 
    theta = 2*np.arctan(np.sqrt((1.0 + e)/(1.0 - e))*np.tan(E/2.0)) 

    # Distance from the Sun to the poin on orbit 
    r = a*(1.0 - e*np.cos(E)) 

    # Cartesian coordinates 
    x = r*(np.cos(node)*np.cos(peri + theta) - np.sin(node)*np.sin(peri + theta)*np.cos(I)) 
    y = r*(np.sin(node)*np.cos(peri + theta) + np.cos(node)*np.sin(peri + theta)*np.cos(I)) 
    z = r*np.sin(peri + theta)*np.sin(I) 

    return x, y, z 


if __name__ == '__main__': 

    # Example orbital elements 
    # a, e, incl, peri, node 
    orb_elements = np.array([ 
     [2.363, 0.515, 4.0, 205.0, 346.1], 
     [0.989, 0.089, 3.1, 55.6, 21.2], 
     [0.898, 0.460, 1.3, 77.1, 331.2], 
     [104.585332285, 0.994914, 89.3950, 130.8767, 282.4633] 
     ]) 

    # Setup the plot 
    fig = plt.figure() 
    ax = fig.gca(projection='3d') 


    # Eccentric anomaly (full range) 
    E = np.linspace(-np.pi, np.pi, 100) 

    # Plot the given orbits 
    for i, orbit in enumerate(orb_elements): 
     a, e, I, peri, node = orbit 

     # Take extra steps in E if the orbit is very large 
     if a > 50: 
      E = np.linspace(-np.pi, np.pi, (a/20.0)*100) 

     # Get the orbit in the cartesian space 
     x, y, z = orbitalElements2Cartesian(a, e, I, peri, node, E) 

     # Plot orbits 
     ax.plot(x, y, z, c='#32CD32') 

    # Add limits (in AU) 
    ax.set_xlim3d(-5,5) 
    ax.set_ylim3d(-5,5) 
    ax.set_zlim3d(-5,5) 

    plt.tight_layout() 
    plt.show() 

我有點傻,這似乎找不到合適的解決方案。我將非常感謝一些幫助!

+0

請提供一個最小的例子(內聯)來重現這一點。 – tacaswell

+0

@tcaswell感謝您的建議!我編輯了這篇文章。 –

+0

看起來像渲染前錯誤地剔除了曲線。穿過整個屏幕的那兩條大線很可能是與投影平面相交的兩段曲線。一點是可見的或非常接近可見區域,另一點不是。由於該段沒有被剔除掉,即使它最有可能位於相機後面(這就是爲什麼被鏡像),它也會在屏幕上進行渲染和內插。至少這是我的賭注。嘗試僅渲染完全可見的段。 – Spektre

回答

0

matplotlib對我的經驗來說,對於複雜的三維圖很不理想(我有類似奇怪的行爲和軸外值)。像mayavi可能是值得考慮的,因爲它是專爲三維圖...

blog中給出了一種可能的解決方法,基本上只是針對您所需的軸爲np.NaN設置軸值。如果我在您的示例中添加以下內容,

for r in [x,y,z]: 
    for i in np.arange(len(r)): 
     if r[i] < -5: 
      x[i] = np.NaN 
      y[i] = np.NaN 
      z[i] = np.NaN 
     elif r[i] > 5: 
      x[i] = np.NaN 
      y[i] = np.NaN 
      z[i] = np.NaN 
     else: 
      pass 

它刪除環繞。

+0

Mayavi看起來很有前途,謝謝!現在,將任何事情都限制在一定範圍之外的問題是,當您從角度角度來看軌道時,您會看到它被剪輯,這對於軌道可視化來說並不是真的可以接受的。這意味着軌道可能是拋物線或雙曲線的,而不是橢圓的。儘管如此,謝謝你的回答! –

0

我有類似的問題,並希望做一些更友好的用戶。我將這個庫中的所有函數都移到了javascript中,並在Three.js中創建了一個webGL接口,它可以讓你在這裏做你想做的事情,而且還可以通過時間函數繪製小行星/彗星的位置。只需要一個Web瀏覽器來使用它。檢查它:)

http://rankinstudio.com/asteroids/asteroids.html