Indice | Previo: Problema2Cuerpos.ProblemaRelativo | Siguiente: Problema2Cuerpos.SolucionTiempo.EcuacionKepler

7.9. La órbita en el espacio

En las secciones precedentes encontramos, usando las constantes de movimiento del problema relativo de dos cuerpos, la ecuación de la curva que describe el vector relativo (Ec. doscuerpos_trayectoria) y de allí las ecuaciones de las trayectorias descritas por los cuerpos individuales (Ecs. doscuerpos_trayectoria_m1 y doscuerpos_trayectoria_m2).

Estos resultados implican que una vez se específica un único parámetro $f$, es posible calcular la posición y velocidad relativa de los dos cuerpos en el sistema de coordenadas natural de la cónica que describen (origen en uno de los focos, semieje $x+$ que coincide con el eje de simetría y apuntando hacia el periapsis.)

En las siguientes secciones intentaremos resolver dos problemas clásicos:

  1. Determinación de órbita. Si se conoce la posición y velocidad relativa instantanea de un sistema de dos cuerpos, referida a un sistema de coordenadas arbitrario, ¿cuál es son los elementos orbitales $(p,e,i,\Omega,\omega,f_0)$ que describen la cónica?

  2. Predicción del vector estado. Una vez se conocen los elementos orbitales de un sistema de dos cuerpos $(p,e,i,\Omega,\omega,f_0)$ ¿cuál es el vector de estado del sistema (y de las partículas constituyentes) para cualquier valor de $f$?

7.9.1. Determinación de la órbita

En la Figura (determinacion_orbita) se muestra la posición y velocidad relativa instantánea de un sistema de dos cuerpos en un sistema de coordenadas arbitrario $x-y-z$. El reto consiste es encontrar a partir únicamente de $\vec r$ y $\dot{\vec r}$, los parámetros geométricos de la cónica $p$ y $e$, su orientación en el espacio específicada por los ángulos $i$, $\Omega$ y $\omega$, además del valor de la anomalía $f_0$ correspondiente a la posición relativa del sistema.

Construcción geométrica requerida para resolver el problema de la determinación de los elementos orbitales de la trayectoria del vector relativo a partir del vector de estado $

Figura 7.84. Construcción geométrica requerida para resolver el problema de la determinación de los elementos orbitales de la trayectoria del vector relativo a partir del vector de estado $

La reconstrucción comienza por encontrar los vectores constantes:

\begin{eqnarray} \label{eq:det_hvec} \vec{h} & = & \vec{r}\times\dot{\vec{r}}\\ \label{eq:det_evec} \vec{e} & = & \frac{\dot{\vec{r}} \times \vec{h}}{\mu} - \frac{\vec{r}}{r}\\ \label{eq:det_nvec} \vec{n} & = & {\hat e}_z\times {\vec h} \\ \end{eqnarray}

De aquí se pueden obtener los parámetros de tamaño y forma de la cónica:

\begin{eqnarray} \label{eq:det_p} p & = & \frac{h^2}{\mu}\\ \label{eq:det_e} e & = & |\vec{e}| \end{eqnarray}

La orientación de la cónica puede obtenerse usando el producto escalar entre los vectores claves:

\begin{eqnarray} \label{eq:det_i} \cos i & = & \frac{\vec{h}\cdot \hat{a_z}}{h} = \frac{h_z}{h} \\ \label{eq:det_W} \cos \Omega & = & \frac{\vec{n}\cdot \hat{a_x}}{n} = \frac{n_x}{n}\\ \label{eq:det_w} \cos \omega & = & \frac{\vec{n}\cdot \vec{e} }{ne} \\ \end{eqnarray}

Numéricamente, la función $\theta=\cos^{-1} x$ devuelve un ángulo en el rango $0\leq\theta\leq\pi$. Esto es suficiente para determinar la inclinación orbital $i$, pero no lo es en el caso de la longitud del nodo ascendente $\Omega$ que es un ángulo en el rango $0\leq\Omega\leq 2\pi$, ni tampoco en el caso del argumento del periapsis que tiene un rango similar. En estos dos casos podemos usar la siguientes reglas para determinar el cuadrante correcto para los ángulos respectivos:

  • Longitud del nodo ascendente $\Omega$. El criterio para determinar el cuadrante correcto de este elemento, parte de calcular la proyección del vector nodal $\vec{n}$ sobre el eje $y$. Si el signo es positivo, el valor de $\Omega$ será el valor principal devuelto por $\cos^{-1}$:

    $$ \Omega_p\equiv\cos^{-1}\left(\frac{n_x}{n}\right) $$

    En caso contrario, usaremos el ángulo complementario a $2\pi$. En síntesis:

    \begin{equation} \label{eq:W_cuadrante} \Omega = \left\{ \begin{array}{ll} \Omega_p & \mathrm{Si}\;n_y \geq 0\\ 2\pi-\Omega_p & \mathrm{Si}\;n_y < 0\\ \end{array} \right. \end{equation}

  • Argumento del periapsis $\omega$. En este caso el valor principal del elemento es:

    $$ \omega_p\equiv\cos^{-1}\left(\frac{\vec{n}\cdot\vec{e}}{en}\right) $$

    para determinar el cuadrante correcto, proyectamos el vector de excentricidad sobre el eje $z$. Con esto $\omega$ queda:

    \begin{equation} \label{eq:w_cuadrante} \omega = \left\{ \begin{array}{ll} \omega_p & \mathrm{Si}\;e_z \geq 0\\ 2\pi-\omega_p & \mathrm{Si}\;e_z < 0\\ \end{array} \right. \end{equation}

Finalmente la anomalía verdadera instantánea se puede determinar de la proyección del vector posición sobre el vector excentricidad:

\begin{equation} \label{eq:det_f} f_p \equiv \cos^{-1}\left(\frac{\vec{e}\cdot \vec{r}}{er}\right) \end{equation}

El cuadrante del ángulo en este caso se encuentra determinando si el vector relativo se aleja ($0<f<\pi$) o se acerca al periapsis ($\pi<f<2\pi$). Esta condición puede evaluarse si se calcula la proyección del vector velocidad sobre el vector posición:

$$ v_r = \frac{\vec{r}\cdot \vec{v}}{r} $$

Finalmente, el valor $f$ en el cuadrante correcto:

\begin{equation} \label{eq:f_cuadrante} f = \left\{ \begin{array}{ll} f_p & \mathrm{Si}\;v_r \geq 0\\ 2\pi-f_p & \mathrm{Si}\;v_r < 0\\ \end{array} \right. \end{equation}

7.9.2. Predicción del vector de estado

Una vez hemos calculado los elementos orbitales $(p,e,i,\Omega,\omega,f_0)$ a partir de una posición específica o nos son provistos para una trayectoria particular, el siguiente problema consiste en determinar la posición y velocidad relativa (vecto de estado) del sistema para un valor cualquiera de la anomalía verdadera $f$.

En la Sección Elementos orbitales habíamos visto cómo calcular las coordenadas de la partícula a partir de los elementos orbitales, para una anomalía verdadera $f$ (Ecs. elementos_estado_f):

$$ \begin{array}{rcl} x & = & r[\cos \Omega \cos(\omega+f) - \cos i \sin \Omega \sin(\omega+f)]\\ y & = & r[\sin \Omega \cos(\omega+f) + \cos i \cos \Omega \sin(\omega+f)]\\ z & = & r \sin i \sin(\omega+f)\\ \end{array} $$

donde $r=p/(1+e\cos f)$

Las componentes cartesianas de la velocidad en el espacio se pueden obtener partiendo de sus componentes en el sistema de coordenadas natural de la cónica (Ecs. v_hodografo):

\begin{eqnarray} \nonumber \dot x''' & = & -\frac{\mu}{h} \sin f\\ \nonumber \dot y''' & = & \frac{\mu}{h} (e+\cos f)\\ \nonumber \dot z''' & = & 0 \end{eqnarray}

que pueden, usando las expresiones explícitas para la rotación en tres dimensiones dadas por las Ecs. (rotacion3d_natural_M) y (matrizM_explicita_transpuesta), rotarse al sistema del observador:

\begin{eqnarray} \label{eq:elementos_dotx} \dot x = & & \frac{\mu}{h}[-\cos \Omega \sin(\omega+f) - \cos i \sin \Omega \cos(\omega+f)]\\ \nonumber & - & \frac{\mu e}{h}(\cos \Omega \sin\omega + \cos\omega\cos i\sin\Omega)\\ \label{eq:elementos_doty} \dot y = & & \frac{\mu}{h}[-\sin \Omega \sin(\omega+f) + \cos i \cos \Omega \cos(\omega+f)]\\ \nonumber & + & \frac{\mu e}{h}(-\sin \Omega \sin\omega + \cos\omega\cos i\cos\Omega)\\ \label{eq:elementos_dotz} \dot z = & & \frac{\mu}{h}[\sin i\cos(\omega+f) + e\cos \omega\sin i]\\ \end{eqnarray}

7.9.3. La órbita osculatriz

La existencia de una función que asigna a cada vector de estado $\vec x:(x,y,z,v_x,v_y,v_z)$ un conjunto de elementos orbitales $p(\vec x)$, $e(\vec x)$, $i(\vec x)$, $\Omega(\vec x)$, $\omega(\vec x)$, $f(\vec x)$ y viceversa (Ecs. det_p-elementos_dotz) es una interesante propiedad que conduce a una importante definición en mecánica celeste.

Ilutración del concepto de órbita osculatriz. La trayectoria del cuerpo (curva continua) no es una cónica. Sin embargo por cada punto de la curva (p.e. los puntos P y Q) podemos encontrar una cónica que sea tangente a la curva (curvas rayadas) y que tenga como foco ($F_P$ o $F_Q$) el origen de coordenadas.

Figura 7.85. Ilutración del concepto de órbita osculatriz. La trayectoria del cuerpo (curva continua) no es una cónica. Sin embargo por cada punto de la curva (p.e. los puntos P y Q) podemos encontrar una cónica que sea tangente a la curva (curvas rayadas) y que tenga como foco ($F_P$ o $F_Q$) el origen de coordenadas.

Independientemente de cuál es la trayectoria que siga una partícula en el espacio, sea esta una cónica o no, siempre es posible encontrar para cada punto de la trayectoria, una curva cónica con elementos $p(\vec x)$, $e(\vec x)$, $i(\vec x)$, $\Omega(\vec x)$, $\omega(\vec x)$, $f(\vec x)$ que tiene como foco el origen de coordenadas y es tangente a la trayectoria de la partícula (ver Figura (osculatriz).) Por su naturaleza "rasante", llamamos a esta cónica la órbita osculatriz o cónica osculatriz (la palabra osculatriz viene del latin osculo que significa "beso".)

Naturalmente, en el problema relativo de los dos cuerpos, la cónica osculatriz asociada a cada punto de la trayectoria descrita por el vector relativo será siempre la misma. En otros términos, podemos decir que los elementos orbitales $p(\vec x)$, $e(\vec x)$, $i(\vec x)$, $\Omega(\vec x)$, $\omega(\vec x)$, $f(\vec x)$ son en sí mismos constantes de movimiento en el problema de los dos cuerpos (un interesante resultado sobre el que volveremos en el Capítulo Formalismo de Hamilton-Jacobi).

En situaciones más generales (como veremos en la Sección Aproximación de dos cuerpos a sistemas jerárquicos), los elementos de la órbita osculatriz cambian punto a punto y el estudio de su dinámica (conocido en mecánica celeste como teoría de perturbaciones) se vuelve en sí mismo muy interesante. Si bien la teoría de perturbaciones está más allá del objetivo de este libro (es de esos temas desarrollados en detalle en la mayoría de los textos avanzados de mecánica celeste), por su importancia abordaremos algunos de sus rudimentos más adelante en este capítulo.

7.9.4. Un ejemplo numérico

Para poner en práctica lo visto en esta sección partamos de un sistema de dos cuerpos con condiciones iniciales arbitrarias:

In [1]:
from numpy import array
sistema=[
    dict(m=1.0,
         r=array([0.0,0.0,+0.3]),
         v=array([+1.0,0.0,0.5])),
    dict(m=0.5,
         r=array([+1.0,0.0,0.0]),
         v=array([0.0,+1.0,0.0])),
]

A partir de estas condiciones iniciales, podemos calcular los parámetros iniciales del sistema relativo:

In [2]:
from numpy.linalg import norm
mu=sistema[0]["m"]+sistema[1]["m"]
rvec=sistema[0]["r"]-sistema[1]["r"]
r=norm(rvec)
vvec=sistema[0]["v"]-sistema[1]["v"]
v=norm(vvec)

mu = 1.5
r_vec = [-1.   0.   0.3]
r = 1.044030650891055
v_vec = [ 1.  -1.   0.5]
v = 1.5

A partir de los vectores posición y velocidad podemos obtener los vectores clave (Ecs. det_h-det_n):

In [4]:
from numpy import cross
from numpy.linalg import norm

#Momento angular relativo específico
hvec=cross(rvec,vvec)
h=norm(hvec)
#Vector excentricidad
evec=cross(vvec,hvec)/mu-rvec/r
e=norm(evec)
#Vector nodo ascendente
nvec=cross([0,0,1],hvec)
n=norm(nvec)

hvec = [0.3 0.8 1. ]
h = 1.3152946437965904
evec = [ 0.02449295 -0.56666667  0.44598545]
e = 0.7215358864417007
nvec = [-0.8  0.3  0. ]
n = 0.8544003745317532

Los parámetros de tamaño y forma se obtienen usando:

In [6]:
#Semilatus rectum
p=h**2/mu
#Semieje mayor
a=p/(1-e**2)
#Velocidad angular promedio
from numpy import sqrt
nmed=sqrt(mu/abs(a)**3)

p = 1.153333333333333 u.c.
a = 2.405855445416549 u.c.
nmed = 0.3282020847560834 u.c.

La orientación de la órbita se obtiene usando las Ecs. (det_W)-(det_f):

In [8]:
from numpy import dot,arccos,pi
i=arccos(hvec[2]/h)

Wp=arccos(nvec[0]/n)
W=Wp if nvec[1]>0 else 2*pi-Wp

wp=arccos(dot(nvec,evec)/(e*n))
w=wp if evec[2]>0 else 2*pi-wp

fp=arccos(dot(rvec,evec)/(r*e))
f0=fp if dot(rvec,vvec)>0 else 2*pi-fp

i = 40.510589437332754 grados
W = 159.44395478041653 grados
w = 107.91123121198778 grados
f_0 = 278.34291953929824 grados

Hemos sintetizado este mismo procedimiento como la rutina estado_a_elementos cuyo algoritmo se presenta en el [Apéndice Algoritmos y rutinas útiles].

Con este resultado podemos ahora determinar el vector posición y velocidad en un valor arbitario de la anomalía verdadera:

In [10]:
#Anomalía verdadera
f=pi/2

#Distancia al punto
from numpy import cos
r=p/(1+e*cos(f))

Las coordenadas del punto son:

In [11]:
from numpy import cos,sin
x=r*(cos(W)*cos(w+f)-cos(i)*sin(W)*sin(w+f))
y=r*(sin(W)*cos(w+f)+cos(i)*cos(W)*sin(w+f))
z=r*sin(i)*sin(w+f)

from numpy import array
r_nuevo=array([x,y,z])

Y la velocidad:

In [12]:
#Parametro mu/h
muh=mu/h

vx=muh*(-cos(W)*sin(w+f)-cos(i)*sin(W)*cos(w+f))\
   -muh*e*(cos(W)*sin(w)+cos(w)*cos(i)*sin(W))
vy=muh*(-sin(W)*sin(w+f)+cos(i)*cos(W)*cos(w+f))\
   +muh*e*(-sin(W)*sin(w)+cos(w)*cos(i)*cos(W))
vz=muh*(sin(i)*cos(w+f)+e*cos(w)*sin(i))

from numpy import array
v_nuevo=array([vx,vy,vz])

Una posible comprobación de que nuestra predicción es la correcta (aunque es una condición necesaria y no suficiente de su validez), es verificar si el momento angular específico relativo con el vector de estado sigue siendo el mismo:

In [13]:
from numpy import cross
hvec_nuevo=cross(r_nuevo,v_nuevo)

Inicial, hvec = [0.3 0.8 1. ]
Nueva posición, hvec = [0.3 0.8 1. ]

Y vemos que coinciden.

Una prueba aún más estricta para nuestra teoría sería la de comparar la órbita predicha con una integración numérica de las ecuaciones de movimiento. Para ello nos valdremos de la rutina ncuerpos_solucion definida en el Alg. (ncuerpos_solucion):

In [15]:
#Tiempo característico de la cónica
#Si la cónica es una elipse este es el período
from numpy import pi
T=2*pi/nmed

#Tiempos 
from numpy import linspace
ts=linspace(0,T,20)

#Solución a las e.d.m. del sistema
from pymcel.export import ncuerpos_solucion
rs,vs,rps,vps,constantes=ncuerpos_solucion(sistema,ts)

Para comparar con la solución necesitamos calcular las posiciones y velocidades relativas. Por relativas, no es importante si usamos el vector de estado referido al sistema de referencia original (rs, vs) o el referido al centro de masa (rps,vps):

In [16]:
rs_num=rs[0,:,:]-rs[1,:,:]
vs_num=vs[0,:,:]-vs[1,:,:]

Para visualizar la órbita predicha en tres dimensiones, usaremos la rutina conica_de_elementos:

In [17]:
#Visualización de la cónica 
from numpy import pi
from pymcel.export import conica_de_elementos
fig=conica_de_elementos(p,e,i*180/pi,W*180/pi,w*180/pi,figreturn=True)
ax=fig.gca()
#Posición y velocidad inicial
ax.quiver(0,0,0,
          rvec[0],rvec[1],rvec[2],
          color='k');
ax.quiver(rvec[0],rvec[1],rvec[2],
          vvec[0],vvec[1],vvec[2],
          color='k');

#Posición y velocidad nueva
ax.quiver(0,0,0,
          r_nuevo[0],r_nuevo[1],r_nuevo[2],
          color='b');
ax.quiver(r_nuevo[0],r_nuevo[1],r_nuevo[2],
          v_nuevo[0],v_nuevo[1],v_nuevo[2],
          color='b');

#Posiciones calculadas numéricamente
ax.plot(rs_num[:,0],rs_num[:,1],rs_num[:,2],'rx');

fig.tight_layout()

Figura 7.86.

Este resultado nos permite comprobar plenamente los resultados teóricos de las últimas sesiones:

  • La posición y velocidad inicial coincide con uno de los puntos de la curva cónica.

  • La nueva posición y velocidad coincide también con uno de los puntos de la cónica.

  • Las posiciones relativas calculadas numéricamente coinciden perfectamente con los puntos sobre la cónica.

Indice | Previo: Problema2Cuerpos.ProblemaRelativo | Siguiente: Problema2Cuerpos.SolucionTiempo.EcuacionKepler