¿Cómo puedo trazar la órbita de un satélite en 3D desde un TLE usando Python y Skyfield?

Código actualizado para 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 es el procedimiento estándar con el que se pretende que trabajen los TLE . Todo lo siguiente es extremadamente útil, pero el punto más importante sería usar un paquete SGP4 estándar y reciente que se recomienda, no intente usar los elementos en un TLE usted mismo . Esto se debe a que el paquete TLE y SGP4 están construidos el uno para el otro.

• pagina web de documentacion
• PDF en formato TLE
• Informe de seguimiento espacial n.º 3 PDF
• Revisando la página web del Informe Spacetrack #3
• Revisión del informe Spacetrack n.º 3 PDF

Un punto interesante es que para órbitas de gran altitud con períodos superiores a 225 minutos, una implementación moderna de SGP4 cambiará a un algoritmo llamado SDP4. Véase la pregunta de correcciones de “Espacio profundo” en SGP4; ¿Cómo explica la gravedad del Sol y la Luna? para más sobre eso.

La forma más fácil de acceder a SGP4 que conozco es en el paquete Skyfield de Python . Puede encontrar rutinas SGP4 en muchos idiomas en muchos lugares. Le recomendaría que elija algo que se use ampliamente y esté bien mantenido, y no cualquier código aleatorio en Internet.

El SGP4 de Skyfield se basa en: https://github.com/brandon-rhodes/python-sgp4

Aquí hay un gráfico de la órbita terrestre baja de la etapa superior de SpaceX F9 conocida como Roadster generada usando Skyfield y TLE de Roadster, por supuesto antes de la segunda quemadura que lo puso en órbita heliocéntrica .

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

La respuesta dada antes parece funcionar muy bien, solo estoy aquí para dar una respuesta diferente si está más interesado en aprender sobre el software que se dedica a las trayectorias de propagación numérica y el trazado 3D.

El procedimiento general para solucionar este problema sería:

1. Leer TLE. Te dará inclinación, RAAN, excentricidad, argumento de periapsis, anomalía media, movimiento medio (revoluciones/día), y algunos otros.
2. Calcule el período orbital (T) en segundos (donde n es el movimiento medio)

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

3. Calcular eje semi-mayor

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

4. Calcule la anomalía excéntrica a través del método de resolución de raíces de Newton (no tiene que ser Newton, pero eso es lo que uso). Dado que esta es una ecuación trascendental, no hay una solución algebraica para la anomalía excéntrica E, por lo que se resuelve iterativamente:

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

5. Calcular anomalía verdadera$\theta$:

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

6. Propague la órbita durante un período de tiempo determinado utilizando las perturbaciones que desee
7. Parcela 3D

Muestro cómo hacer los pasos 1-5 en un video sobre TLE, usando un ejemplo de TLE de ISS.

El paso 6 es un poco más complicado y requiere un poco de conocimiento sobre ecuaciones diferenciales ordinarias (EDO) y solucionadores de ODE, que cubro en esta serie: https://www.youtube.com/playlist?list=PLOIRBaljOV8gn074rWFWYP1dCr2dJqWab

Otra extensión de este tipo de análisis sería también calcular y trazar las trayectorias terrestres con respecto a la superficie de la Tierra. Esto es nuevamente más complicado cuando se trata de software, pero aprenderá mucho si sigue el proceso de lo que sucede debajo del capó de los grandes paquetes de Python: https://www.youtube.com/playlist?list=PLOIRBaljOV8ghaN9XcC4ubv- QLPOQByYQ