from skyfield.api import load, EarthSatellite
from skyfield.elementslib import osculating_elements_of
import numpy as np
import plotly.express as px
from IPython.display import Image
Sistema de coordenadas¶
Utilizaré, como es habitual en astrodinámica, el Sistema de Coordenadas Geocéntrico Ecuatorial (GECS), también conocido como $IJK$ ya que es así como se denominan los ejes de coordenadas en este sistema. Su origen se situa en el centro de la Tierra, el plano fundamental (plano $IJ$) es el del ecuador terrestre. La dirección fundamental (eje $I$) es la del equinoccio vernal, y el eje $K$ es perpendicular al plano fundamental y por lo tanto coincide con el eje de rotación terrestre. Se trata de un sistema de referencia fijo que no gira con la Tierra.
En realidad tanto el ecuador terrestre como el equinoccio de primavera van variando con el tiempo y por ello estas posiciones se suelen referir a una época determinada, en la actualidad el año 2000, en lo que se conoce como época J2000.0. En tal caso el sistema $IJK$ se suele utilizar de forma intercambiable con el Earth Centered Equatorial (ECI), el cual a su vez coincide con gran aproximación con los planos y direcciones del sistema ICRF (International Celestial Reference Frame), sistema que se establece tomando como referencia objetos celestes muy lejanos y que por lo tanto podemos suponer fijos en la práctica.
Elementos orbitales clásicos¶
Por elementos orbitales entenderemos un conjunto de valores numéricos que permiten describir el tamaño, forma y orientación de la órbita de un satélite, así como la posición del satélite en la órbita en un instante determinado. A continuación introduciremos los que en su momento definió Kepler y que actualmente se conocen como elementos orbitales clásicos. También veremos como se pueden obtener a partir de los TLE de un satélite utilizando la librería Skyfield.
Comencemos por obtener los TLE de un satélite que nos sirva de ejemplo, siguiendo el método ya visto en una entrada anterior de este blog.
ISS_TLE = '''ISS (ZARYA)
1 25544U 98067A 21103.84943184 .00000176 00000-0 11381-4 0 9990
2 25544 51.6434 300.9481 0002858 223.8443 263.8789 15.48881793278621
'''
L0, L1, L2 = ISS_TLE.splitlines()
ISS = EarthSatellite(L1, L2, L0) # La línea conteniendo el nombre va en último lugar
ISS
<EarthSatellite ISS (ZARYA) catalog #25544 epoch 2021-04-13 20:23:11 UTC>
A continuación sigue una breve descripción de cada elemento orbital, con un ejemplo de su obtención con la librería Skyfield. El esquema siguiente puede servir de ayuda para situar geométricamente algunos de ellos.
Image(filename='Datos/elementos_orbitales_V2.png')
pos = ISS.at(ISS.epoch) # Objeto posición de Skyfield para el satélite en la fecha de los TLE
pos
<Geocentric ICRS position and velocity at date t center=399 target=-125544>
A partir de la posición se puede generar el objeto siguiente cuyos atributos serán precisamente los elementos orbitales.
elementos = osculating_elements_of(pos)
Semieje mayor $a$¶
Es la mitad del diametro mayor de la órbita elíptica del satélite
a = elementos.semi_major_axis.km # Semieje mayor en Km
a
6797.148581523243
Excentricidad $\epsilon$¶
Es el cociente entre la distancia entre los focos de la elipse y su diámetro mayor. Si designamos por $c$ la distancia del centro de la elipse a uno de sus focos, y por $a$ su semieje mayor, la excentricidad viene dada por $\epsilon = c/a$. Es un valor comprendido entre $0$ y $1$, siendo $0$ el valor correspondiente a una circunferencia. La excentricidad indica la forma de la órbita, tanto más ovalada cuanto más próxima su excentricidad al valor $1$.
e = elementos.eccentricity
e
0.0007840496617501508
Se trata de un valor muy próximo a $0$ por lo cual la órbita es prácticamente circular.
Inclinación $i$¶
Es el ángulo formado entre el plano fundamental (plano ecuatorial) y el plano de la órbita.
i = elementos.inclination # Objeto ángulo
i.degrees # en grados, .radians devuelve el ángulo en radianes
51.53787056146399
Ascensión recta del nodo ascendente: $\Omega$¶
El plano de la órbita corta al plano del ecuador terrestre a lo largo de una línea conocida como línea de nodos. Los dos puntos de corte de la órbita con la línea de nodos son los nodos de la órbita. En particular el nodo ascendente es aquel en el cual el satélite, en su movimiento, pasa de sur a norte.
Este parámetro, representado por $\Omega$, es la ascensión recta (es decir, el ángulo a lo largo del ecuador, desde la dirección principal, medida hacia el este), con el nodo ascendente. Su rango varía entre 0 y 360º. Su valor permite determinar la línea de nodos y por lo tanto, junto con la inclinación, establece la posición en el espacio del plano inclinado de la órbita.
Omega = elementos.longitude_of_ascending_node
Omega.degrees
300.6975553670492
Argumento del perigeo: $\omega$¶
Representado por $\omega$, es el ángulo, medido a lo largo de la órbita en la dirección del movimiento, entre el nodo ascendente y el perigeo (punto de la órbita más cercano a la Tierra). $0º \leq \omega < 360º$. Este elemento permite situar la órbita (cuya forma ya es conocida gracias a los elementos $a$ e $i$) dentro de su plano orbital.
omega = elementos.argument_of_periapsis # Argumento del perigeo
omega.degrees
85.40540448984416
Anomalía verdadera: $\nu$¶
$\nu$ es el parámetro que nos permite localizar el satélite dentro de la órbita, de modo que es el único parámetro que depende del tiempo. Es el ángulo que forma el vector de posición del satélite con el vector de posición del perigeo en un momento particular $t_0$ conocido como epoch. recordemos que el objeto posición lo hemos creado para el instante "epoch" del TLE utilizado. Este instante es el tiempo de referencia del conjunto TLE del que hemos partido y representa el instante de tiempo en el cual es mayor su exactitud. En este caso corresponde a la fecha y hora siguiente:
elementos.time.utc_iso()
'2021-04-13T20:23:11Z'
la anomalía verdadera en ese instante será:
nu = elementos.true_anomaly
nu.degrees
42.281676731490144
elementos.time.utc_iso()
'2021-04-13T20:23:11Z'
Situaciones excepcionales¶
En el caso de órbitas circulares, la periapsis no está definida, y se suelen emplear: en lugar de la anomalía verdadera el argumento de latitud $u$, siendo en tal caso innecesario el otro elemento que depende de la posiión del perigeo, el argumento del perigeo.
u = elementos.argument_of_latitude
u.degrees
127.6870812213343
Otra situación especial es el caso de una órbita en el plano ecuatorial. En este caso no está definida la línea de nodos y por tanto no es posible definir la ascensión recta del nodo ascendente $\Omega$ ni el argumento del perigeo $\omega$. Sustituiremos ambos por la longitud del perigeo: $\Pi$
Pi = elementos.longitude_of_periapsis
Pi.degrees
26.102959856893364
Y, en el caso de una órbita circular ecuatorial no es posible definir la ascensión recta del nodo ascendente $\Omega$ (no hay línea de nodos), ni tampoco el argumento del perigeo $\omega$ o la anomalía verdadera $\nu$ ya que el perigeo no está definido. En tales casos utilizaremos como elemento en sustitución de los tres la longitud verdadera $l$.
l = elementos.true_longitude
l.degrees
68.38463658838351
Periodo del satélite¶
Para concluir, mencionaré otro atributo del objeto "elementos" que puede ser útil en muchas situaciones, el cual proporciona el periodo del satélite expresado en días:
elementos.period_in_days
0.0645486831407206
No hay comentarios:
Publicar un comentario