Как я могу построить орбиту спутника в 3D из TLE, используя Python и Skyfield?

Код обновлен для python3*

def makecubelimits(axis, centers=None, hw=None):
    lims = ax.get_xlim(), ax.get_ylim(), ax.get_zlim()
    if centers == None:
        centers = [0.5*sum(pair) for pair in lims] 

    if hw == None:
        widths  = [pair[1] - pair[0] for pair in lims]
        hw      = 0.5*max(widths)
        ax.set_xlim(centers[0]-hw, centers[0]+hw)
        ax.set_ylim(centers[1]-hw, centers[1]+hw)
        ax.set_zlim(centers[2]-hw, centers[2]+hw)
        print("hw was None so set to:", hw)
    else:
        try:
            hwx, hwy, hwz = hw
            print("ok hw requested: ", hwx, hwy, hwz)

            ax.set_xlim(centers[0]-hwx, centers[0]+hwx)
            ax.set_ylim(centers[1]-hwy, centers[1]+hwy)
            ax.set_zlim(centers[2]-hwz, centers[2]+hwz)
        except:
            print("nope hw requested: ", hw)
            ax.set_xlim(centers[0]-hw, centers[0]+hw)
            ax.set_ylim(centers[1]-hw, centers[1]+hw)
            ax.set_zlim(centers[2]-hw, centers[2]+hw)

    return centers, hw

TLE = """1 43205U 18017A   18038.05572532 +.00020608 -51169-6 +11058-3 0  9993
2 43205 029.0165 287.1006 3403068 180.4827 179.1544 08.75117793000017"""
L1, L2 = TLE.splitlines()

from skyfield.api import Loader, EarthSatellite
from skyfield.timelib import Time
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

halfpi, pi, twopi = [f*np.pi for f in [0.5, 1, 2]]
degs, rads = 180/pi, pi/180

load = Loader('~/Documents/fishing/SkyData')
data = load('de421.bsp')
ts   = load.timescale()

planets = load('de421.bsp')
earth   = planets['earth']

Roadster = EarthSatellite(L1, L2)

print(Roadster.epoch.tt)
hours = np.arange(0, 3, 0.01)

time = ts.utc(2018, 2, 7, hours)

Rpos    = Roadster.at(time).position.km
Rposecl = Roadster.at(time).ecliptic_position().km

print(Rpos.shape)

re = 6378.

theta = np.linspace(0, twopi, 201)
cth, sth, zth = [f(theta) for f in [np.cos, np.sin, np.zeros_like]]
lon0 = re*np.vstack((cth, zth, sth))
lons = []
for phi in rads*np.arange(0, 180, 15):
    cph, sph = [f(phi) for f in [np.cos, np.sin]]
    lon = np.vstack((lon0[0]*cph - lon0[1]*sph,
                     lon0[1]*cph + lon0[0]*sph,
                     lon0[2]) )
    lons.append(lon)

lat0 = re*np.vstack((cth, sth, zth))
lats = []
for phi in rads*np.arange(-75, 90, 15):
    cph, sph = [f(phi) for f in [np.cos, np.sin]]
    lat = re*np.vstack((cth*cph, sth*cph, zth+sph))
    lats.append(lat)

if True:    
    fig = plt.figure(figsize=[10, 8])  # [12, 10]

    ax  = fig.add_subplot(1, 1, 1, projection='3d')

    x, y, z = Rpos
    ax.plot(x, y, z)
    for x, y, z in lons:
        ax.plot(x, y, z, '-k')
    for x, y, z in lats:
        ax.plot(x, y, z, '-k')

    centers, hw = makecubelimits(ax)

    print("centers are: ", centers)
    print("hw is:       ", hw)

    plt.show()

r_Roadster = np.sqrt((Rpos**2).sum(axis=0))
alt_roadster = r_Roadster - re

if True:
    plt.figure()
    plt.plot(hours, r_Roadster)
    plt.plot(hours, alt_roadster)
    plt.xlabel('hours', fontsize=14)
    plt.ylabel('Geocenter radius or altitude (km)', fontsize=14)
    plt.show()

SGP4 — это стандартная процедура, для которой предназначены TLE . Все нижеследующее чрезвычайно полезно, но наиболее важным моментом будет использование стандартного, последнего пакета SGP4, который рекомендуется, не пытайтесь использовать элементы в TLE самостоятельно . Это связано с тем, что пакеты TLE и SGP4 созданы друг для друга.

• Веб-страница документации
• PDF в формате TLE
• Отчет о космическом треке № 3 PDF
• Повторное посещение веб-страницы отчета Spacetrack # 3
• Пересмотр отчета Spacetrack #3 PDF

Один интересный момент заключается в том, что для высотных орбит с периодом более 225 минут современная реализация SGP4 переключится на алгоритм, который называется SDP4. См. исправления вопроса «Глубокий космос» в SGP4; как это объясняет гравитацию Солнца и Луны? подробнее об этом.

Самый простой доступ к SGP4, который я знаю, находится в Python-пакете Skyfield . Вы можете найти подпрограммы SGP4 на многих языках во многих местах. Я бы порекомендовал вам выбрать то, что широко используется и хорошо поддерживается, а не просто случайный код в Интернете.

Skyfield SGP4 основан на: https://github.com/brandon-rhodes/python-sgp4 .

Вот график низкой околоземной орбиты разгонного блока SpaceX F9, известного как Roadster , созданный с использованием TLE Skyfield и Roadster, конечно, до второго запуска, который вывел его на гелиоцентрическую орбиту .

def makecubelimits(axis, centers=None, hw=None):
    lims = ax.get_xlim(), ax.get_ylim(), ax.get_zlim()
    if centers == None:
        centers = [0.5*sum(pair) for pair in lims] 

    if hw == None:
        widths  = [pair[1] - pair[0] for pair in lims]
        hw      = 0.5*max(widths)
        ax.set_xlim(centers[0]-hw, centers[0]+hw)
        ax.set_ylim(centers[1]-hw, centers[1]+hw)
        ax.set_zlim(centers[2]-hw, centers[2]+hw)
        print("hw was None so set to:", hw)
    else:
        try:
            hwx, hwy, hwz = hw
            print("ok hw requested: ", hwx, hwy, hwz)

            ax.set_xlim(centers[0]-hwx, centers[0]+hwx)
            ax.set_ylim(centers[1]-hwy, centers[1]+hwy)
            ax.set_zlim(centers[2]-hwz, centers[2]+hwz)
        except:
            print("nope hw requested: ", hw)
            ax.set_xlim(centers[0]-hw, centers[0]+hw)
            ax.set_ylim(centers[1]-hw, centers[1]+hw)
            ax.set_zlim(centers[2]-hw, centers[2]+hw)

    return centers, hw

TLE = """1 43205U 18017A   18038.05572532 +.00020608 -51169-6 +11058-3 0  9993
2 43205 029.0165 287.1006 3403068 180.4827 179.1544 08.75117793000017"""
L1, L2 = TLE.splitlines()

from skyfield.api import Loader, EarthSatellite
from skyfield.timelib import Time
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
degs, rads = 180/pi, pi/180

load = Loader('~/Documents/fishing/SkyData')
data = load('de421.bsp')
ts   = load.timescale()

planets = load('de421.bsp')
earth   = planets['earth']

Roadster = EarthSatellite(L1, L2)

print(Roadster.epoch.tt)
hours = np.arange(0, 3, 0.01)

time = ts.utc(2018, 2, 7, hours)

Rpos    = Roadster.at(time).position.km
Rposecl = Roadster.at(time).ecliptic_position().km

print(Rpos.shape)

re = 6378.

theta = np.linspace(0, twopi, 201)
cth, sth, zth = [f(theta) for f in (np.cos, np.sin, np.zeros_like)]
lon0 = re*np.vstack((cth, zth, sth))
lons = []
for phi in rads*np.arange(0, 180, 15):
    cph, sph = [f(phi) for f in (np.cos, np.sin)]
    lon = np.vstack((lon0[0]*cph - lon0[1]*sph,
                     lon0[1]*cph + lon0[0]*sph,
                     lon0[2]) )
    lons.append(lon)

lat0 = re*np.vstack((cth, sth, zth))
lats = []
for phi in rads*np.arange(-75, 90, 15):
    cph, sph = [f(phi) for f in (np.cos, np.sin)]
    lat = re*np.vstack((cth*cph, sth*cph, zth+sph))
    lats.append(lat)

if True:    
    fig = plt.figure(figsize=[10, 8])  # [12, 10]

    ax  = fig.add_subplot(1, 1, 1, projection='3d')

    x, y, z = Rpos
    ax.plot(x, y, z)
    for x, y, z in lons:
        ax.plot(x, y, z, '-k')
    for x, y, z in lats:
        ax.plot(x, y, z, '-k')

    centers, hw = makecubelimits(ax)

    print("centers are: ", centers)
    print("hw is:       ", hw)

    plt.show()

r_Roadster = np.sqrt((Rpos**2).sum(axis=0))
alt_roadster = r_Roadster - re

if True:
    plt.figure()
    plt.plot(hours, r_Roadster)
    plt.plot(hours, alt_roadster)
    plt.xlabel('hours', fontsize=14)
    plt.ylabel('Geocenter radius or altitude (km)', fontsize=14)
    plt.show()

Ответ, данный ранее, кажется, отлично работает, я просто здесь, чтобы дать другой ответ, если вы больше заинтересованы в изучении программного обеспечения, которое используется для численных траекторий распространения и трехмерного построения.

Общая процедура решения этой проблемы будет следующей:

1. Читайте ТЛЭ. Это даст вам наклон, RAAN, эксцентриситет, аргумент перицентра, среднюю аномалию, среднее движение (оборотов в день) и некоторые другие.
2. Рассчитайте период обращения (T) в секундах (где n — среднее движение)

$T=\frac{24.0 * 3600.0}{n}$

3. Вычислить большую полуось

$sma = (\frac{T^2 \mu}{4\pi^2})^\frac{1}{3}$

4. Рассчитайте эксцентрическую аномалию с помощью метода решения корней Ньютона (это не обязательно должен быть метод Ньютона, но я его использую). Поскольку это трансцендентное уравнение, для эксцентрической аномалии E нет алгебраического решения, поэтому оно решается итеративно:

$M_{e} = E - esin(E)$

5. Рассчитать истинную аномалию$\theta$:

$\theta = 2\space arctan(\sqrt{\frac{1+e}{1-e}}tan(\frac{E}{2}))$

6. Распространяйте орбиту в течение заданного времени, используя любые возмущения, которые вы хотите
7. Сюжет 3D

Я покажу, как выполнять шаги 1–5 в видео о TLE на примере ISS TLE.

Шаг 6 немного сложнее и требует некоторых знаний об обыкновенных дифференциальных уравнениях (ОДУ) и решателях ОДУ, о которых я расскажу в этой серии: https://www.youtube.com/playlist?list=PLOIRBaljOV8gn074rWFWYP1dCr2dJqWab .

Другим расширением этого типа анализа может быть также расчет и построение траекторий движения относительно поверхности Земли. Это снова более сложно, когда дело доходит до программного обеспечения, но вы многое узнаете, если пройдете процесс того, что происходит под капотом больших пакетов Python: https://www.youtube.com/playlist?list=PLOIRBaljOV8ghaN9XcC4ubv- QLPOQByYQ