Como hemos deducido en la página, la ecuación de la trayectoria en coordenadas polares es

${\displaystyle r=\frac{d}{1+\epsilon \cos \theta }\epsilon =\sqrt{1+\frac{2L^{2}E}{m^3{G}^2{M}^2}}d=\frac{L^2}{GM{m}^2}}$

Por ser la fuerza de atracción una fuerza central, el momento angular L es constante, y por ser conservativa la energía E es constante. El momento angular en coordenadas polares se escribe

${\displaystyle L=m{r}^2\frac{d\theta }{dt}}$

Integrando, obtenemos la expresión de la posición angular θ en función del tiempo t.

${\displaystyle {\displaystyle \underset{0}{\overset{\theta }{\int }}\frac{d\theta }{{\left(1+\epsilon \cos \theta \right)}^2}}=\frac{G^2{M}^2{m}^3}{L^3}t}$

Vamos a calcular esta integral para las trayectorias elípticas, parabólicas e hiperbólicas

Trayectorias elípticas

Supongamos un planeta, como la Tierra que se mueve en una órbita elíptica de semieje mayor a y de excentricidad ε=c/a con el Sol situado en el foco F. Sea el periodo, o tiempo que tarda en dar una vuelta completa, P. Al cabo de un cierto tiempo t, el planeta se encuentra en una posición (r, θ). El ángulo θ se denomina anomalía verdadera en el instante t.

Tomamos el tiempo t=0, cuando el planeta pasa por la posición más cercana al centro de fuerzas (perihelio o perigeo), θ=0.

Calculamos la integral haciendo el cambio de variable

${\displaystyle \begin{array}{l}\cos \theta =\frac{1-x^2}{1+x^2}\text{,\hspace{1em}}\sin \theta =\frac{2x}{1+x^2}\text{,\hspace{1em}}d\theta =\frac{2}{1+x^2}dx\text{,\hspace{1em}}x=\tan \left(\frac{\theta }{2}\right)\hfill \\ {{\displaystyle \int \frac{d\theta }{{\left(1+\epsilon \cos \theta \right)}^2}}}^{\text{}}=\frac{2}{{\left(1+\epsilon \right)}^2}{\displaystyle \int \frac{1+x^2}{{\left(1+a{x}^2\right)}^2}dx}=\hfill \\ \frac{2}{{\left(1+\epsilon \right)}^2}\left\{\frac{1}{a}{{\displaystyle \int \frac{dx}{1+a{x}^2}}}^{\text{}}+\frac{a-1}{a}{\displaystyle \int \frac{dx}{{\left(1+a{x}^2\right)}^2}}\right\}\hfill \end{array}}$

con a=(1-ε)/(1+ε), no confundir esta constante con el semieje mayor

La primera integral es inmediata

${\displaystyle {{\displaystyle \int \frac{dx}{1+a{x}^2}}}^{\text{}}=\frac{1}{\sqrt{a}}{\displaystyle \int \frac{du}{1+u^2}}=\frac{1}{\sqrt{a}}\arctan u=\frac{1}{\sqrt{a}}\arctan \left(\sqrt{a}x\right)}$

La segunda integral es algo más laboriosa

${\displaystyle {{\displaystyle \int \frac{dx}{{\left(1+a{x}^2\right)}^2}}}^{\text{}}=\frac{1}{\sqrt{a}}{\displaystyle \int \frac{du}{{\left(1+u^2\right)}^2}}}$

Haciendo el cambio, u=tanv, du=dv/cos²(v)

${\displaystyle {\displaystyle \int \frac{du}{{\left(1+u^2\right)}^2}}={\displaystyle \int {\cos }^{2}v·dv}={\displaystyle \int \frac{1+\cos (2v)}{2}dv}=\frac{1}{2}v+\frac{1}{4}\sin (2v)}$

Deshaciendo los cambios

${\displaystyle \begin{array}{l}{{\displaystyle \int \frac{dx}{{\left(1+a{x}^2\right)}^2}}}^{\text{}}=\frac{1}{2\sqrt{a}}\left\{\arctan \left(\sqrt{a}x\right)+\frac{1}{2}\sin \left(2\arctan (\sqrt{a}x)\right)\right\}=\hfill \\ \frac{1}{2\sqrt{a}}\left(\arctan (\sqrt{a}x)+\frac{\sqrt{a}x}{1+a{x}^2}\right)\hfill \end{array}}$

Se ha tenido en cuenta la relación trigonométrica entre sin(2u) y tan(u)

${\displaystyle \sin \left(2u\right)=\frac{2\tan u}{1+{\tan }^{2}u}}$

El resultado de la suma de las dos integrales es

${\displaystyle \begin{array}{l}{\displaystyle \int \frac{d\theta }{{\left(1+\epsilon \cos \theta \right)}^2}}=\frac{2}{{\left(1+\epsilon \right)}^2}\left\{\frac{1}{a}{{\displaystyle \int \frac{dx}{1+a{x}^2}}}^{\text{}}+\frac{a-1}{a}{{\displaystyle \int \frac{dx}{{\left(1+a{x}^2\right)}^2}}}^{\text{}}\right\}=\hfill \\ \frac{1}{{\left(1+\epsilon \right)}^2}\left\{\frac{a+1}{a\sqrt{a}}\arctan \left(\sqrt{a}x\right)+\frac{a-1}{2a\sqrt{a}}\frac{2\sqrt{a}x}{1+a{x}^2}\right\}=\hfill \end{array}}$

Definimos el ángulo E tal que

${\displaystyle \tan \left(\frac{E}{2}\right)=\sqrt{a}x}$

y teniendo en cuenta la relación trigonométrica entre sin(E) y tan(E/2), la integral queda simplificada

${\displaystyle \frac{1}{{\left(1+\epsilon \right)}^{2}a\sqrt{a}}\left\{\left(a+1\right)\frac{E}{2}+\frac{\left(a-1\right)}{2}\sin E\right\}=\frac{1}{{\left(1-{\epsilon }^2\right)}^{3/2}}\left(E-\epsilon \sin E\right)}$

La ecuación de Kepler se escribe

${\displaystyle \begin{array}{l}{\displaystyle \underset{\theta }{\overset{0}{\int }}\frac{d\theta }{{(1+\epsilon \cos \theta )}^2}}=\frac{G^2{M}^2{m}^3}{L^3}t\hfill \\ E-\epsilon \sin E={\left(1-{\epsilon }^2\right)}^{3/2}\frac{G^2{M}^2{m}^3}{L^3}t\hfill \end{array}}$

Donde el ángulo E se ha definido como

${\displaystyle \tan \left(\frac{E}{2}\right)=\sqrt{\frac{1-\epsilon }{1+\epsilon }}\tan \left(\frac{\theta }{2}\right)}$

El periodo P es el tiempo que tarda en dar una vuelta completa un planeta, que describe una órbita elíptica

${\displaystyle \begin{array}{l}{P}^2=4{\pi }^2\frac{a^3}{GM}\hfill \\ P=2\pi \frac{L^3}{G^2{M}^2{m}^3}\frac{1}{{(1-{\epsilon }^2)}^{3/2}}\hfill \end{array}}$

La forma más simplificada de la ecuación de Kepler para las trayectorias elípticas, es

${\displaystyle E-\epsilon \sin E=2\pi \frac{t}{P}}$

El miembro izquierdo se denota como M_e (anomalía media)

Representamos M_e en función de la posición angular θ, para los siguientes valores de la excentricidad ε=0, 0.2, 0.5 0.8 y 0.9

x=(0:360)*pi/180;
hold on
for e=[0,0.2,0.5,0.8,0.9]
    y=2*atan(sqrt((1-e)/(1+e))*tan(x/2))-e*sqrt(1-e^2)*sin(x).
/(1+e*cos(x));
    yy=y.*(y>0)+(2*pi+y).*(y<0);
    plot(x,yy, 'displayName',num2str(e))
end
legend('-DynamicLegend','Location','northwest')
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','\pi/2','\pi','3\pi/2','2\pi'})
set(gca,'YTick',0:pi/2:2*pi)
set(gca,'YTickLabel',{'0','\pi/2','\pi','3\pi/2','2\pi'})
hold off
grid on
xlabel('\theta')
ylabel('M_e')
title('La ecuación del tiempo')

Ecuación de la trayectoria en términos del ángulo E

La distancia más cercana al Sol r₁ se obtiene en la ecuación de la trayectoria con θ=0 , y la más alejada r₂ con θ=π. Con r₁+r₂=2a. Siendo a el semieje mayor de la elipse

${\displaystyle \begin{array}{l}r=\frac{d}{1+\epsilon \cos \theta }\\ \frac{d}{1+\epsilon }+\frac{d}{1-\epsilon }=2a,d=a\left(1-{\epsilon }^2\right)\end{array}}$

con la excentricidad ε<1. La ecuación de la trayectoria en términos del ángulo θ es

${\displaystyle r=\frac{a\left(1-{\epsilon }^2\right)}{1+\epsilon \cos \theta }}$

Conociendo la relación entre los ángulos θ y E

${\displaystyle \tan \left(\frac{\theta }{2}\right)=\sqrt{\frac{1+\epsilon }{1-\epsilon }}\tan \left(\frac{E}{2}\right)}$

Determinaremos la ecuación de la trayectoria en términos del ángulo E, utilizando las relaciones trigonométricas

${\displaystyle {\cos }^2\theta =\frac{1}{1+{\tan }^2\theta },\tan \theta =\frac{2\tan \left(\frac{\theta }{2}\right)}{1-{\tan }^2\left(\frac{\theta }{2}\right)}}$

El resultado es

${\displaystyle \cos \theta =\frac{1-\frac{1+\epsilon }{1-\epsilon }{\tan }^2\left(\frac{E}{2}\right)}{1+\frac{1+\epsilon }{1-\epsilon }{\tan }^2\left(\frac{E}{2}\right)}=\frac{\left(1-\epsilon \right){\cos }^2\left(\frac{E}{2}\right)-\left(1+\epsilon \right){\sin }^2\left(\frac{E}{2}\right)}{\left(1-\epsilon \right){\cos }^2\left(\frac{E}{2}\right)+\left(1+\epsilon \right){\sin }^2\left(\frac{E}{2}\right)}=\frac{\cos E-\epsilon }{1-\epsilon \cos E}}$

Hemos utilizado las relaciones trigonométricas

${\displaystyle \begin{array}{l}{\cos }^2\left(\frac{\theta }{2}\right)+{\sin }^2\left(\frac{\theta }{2}\right)=1\\ \cos \theta ={\cos }^2\left(\frac{\theta }{2}\right)-{\sin }^2\left(\frac{\theta }{2}\right)\end{array}}$

Finalmente, la ecuación de la trayectoria es

${\displaystyle r=\frac{a\left(1-{\epsilon }^2\right)}{1+\epsilon \frac{\cos E-\epsilon }{1-\epsilon \cos E}}=a\left(1-\epsilon \cos E\right)}$

Dada la posición angular θ del satélite determinar el tiempo t

Teniendo en cuenta las relaciones trigonométricas

${\displaystyle \sin E=\frac{2\tan \left(\frac{E}{2}\right)}{1+{\tan }^2\left(\frac{E}{2}\right)}\cos E=\frac{1-{\tan }^2\left(\frac{E}{2}\right)}{1+{\tan }^2\left(\frac{E}{2}\right)}}$

Dada la posición angular θ, determinamos el ángulo E mediante las relaciones

${\displaystyle \sin E=\frac{\sqrt{1-{\epsilon }^2}\sin \theta }{1+\epsilon \cos \theta }\cos E=\frac{\epsilon +\cos \theta }{1+\epsilon \cos \theta }}$

Dependiendo del signo del seno y del coseno, se determina el cuadrante y se calcula el ángulo E. Sustituyendo en la ecuación de Kepler, obtenemos M_e

${\displaystyle M_e=E-\epsilon \sin E}$

Conocido M_e, obtenemos el tiempo

${\displaystyle t=\frac{M_e}{2\pi }P\text{,\hspace{1em}}{P}^2=\frac{4{\pi }^2}{GM}{a}^3}$

Ejemplo

Un satélite artificial describe una órbita elíptica con cuyas distancias máxima y mínima al centro de la Tierra son, respectivamente r₁=9.6·10⁶ m y r₂=21·10⁶. Determinar el instante t en el que la posición del satélite es θ=2π/3 (120°)

Elaboramos un script para realizar los cálculos

M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G

r1=9.6e6; %distancia mínima
r2=21e6; %distancia máxima
a=(r1+r2)/2; %semieje mayor
c=a-r1; %semidistancia focal
ex=c/a; %excentricidad
th=2*pi/3; %posición

c_E=(ex+cos(th))/(1+ex*cos(th));
s_E=sqrt(1-ex^2)*sin(th)/(1+ex*cos(th));
E=atan2(s_E,c_E);
if E<0
    E=2*pi+E;
end
t=(E-ex*s_E)*a^(3/2)/sqrt(G*M);
disp(t)

  4.0757e+03

El instante en el que se alcanza la posición angular θ=2π/3 es t=4076 s

Dado el tiempo t determinar la posición angular θ del satélite

Sustituyendo M_e en la ecuación de Kepler, obtenemos una ecuación trascendente cuya raíz tenemos que determinar por procedimientos numéricos

${\displaystyle E-\epsilon \sin E=M_e}$

Un procedimiento apropiado para resolver esta ecuación trascendente es el método de Newton. Sea f(x)=0. El proceso iterativo para calcular la raíz es

${\displaystyle \begin{array}{l}{x}_{i+1}=x_i-\frac{f\left(x_i\right)}{f\text{'}\left(x_i\right)}\\ E_{i+1}=E_i-\frac{E_i-\epsilon \sin E_i-M_e}{1-\epsilon \cos E_i}\end{array}}$

Se repite el proceso iterativo hasta que se encuentra la raíz E con el nivel de precisión deseado. Una vez determinado E se calcula la posición angular θ mediante la relación

${\displaystyle \tan \left(\frac{\theta }{2}\right)=\sqrt{\frac{1+\epsilon }{1-\epsilon }}\tan \left(\frac{E}{2}\right)}$

Ejemplo

Encontrar posición angular θ del satélite del ejercicio anterior en el instante t=10800 s (3 horas)

${\displaystyle \begin{array}{l}{M}_e=2\pi \frac{t}{P}\\ M_e=2\pi \frac{10800}{18828}=3.60\end{array}}$

Resolvemos la ecuación trascendente E-ε·sinE=M_e por el procedimiento de Newton tomando como valor inicial de partida E₀=M_e-ε/2=3.418

x0=3.418;
e=0.3725;
Me=3.604;
while(1)
    x=x0-(x0-e*sin(x0)-Me)/(1-e*cos(x0));
    if abs((x-x0)/x0)<1e-6
        break;
    end
    x0=x;
end
fprintf('La raíz buscada es %1.3f\n',x0)

La raíz buscada es 3.480

Si E=3.480, a partir de la relación

${\displaystyle \tan \left(\frac{\theta }{2}\right)=\sqrt{\frac{1+0.37}{1-0.37}}\tan \left(\frac{3.48}{2}\right)\text{,\hspace{1em}}\theta =3.372}$

obtenemos la posición angular θ=3.372 del satélite en el instante t=10800 s (3 horas)

Interpretación geométrica del ángulo E

Trazamos la trayectoria elíptica de semieje mayor a y excentricidad ε en color rojo. c=ε·a es la semidistancia focal entre O y el foco F donde se encuentra el Sol.

Trazamos una circunferencia de radio a en color azul y dibujamos una recta perpendicular al eje mayor que corta a la circunferencia en el punto Q. El ángulo de la recta OQ se denomina E, anomalía excéntrica del planeta en el instante t. Obtenemos la siguiente relación geométrica: acosE=c+rcosθ

${\displaystyle \begin{array}{l}a\cos E=a\epsilon +\frac{a\left(1-{\epsilon }^2\right)}{1+\epsilon \cos \theta }\cos \theta \\ \cos E=\frac{\epsilon +\cos \theta }{1+\epsilon \cos \theta }\\ {\sin }^{2}E+{\cos }^{2}E=1\text{,\hspace{1em}}\sin E=\frac{\sqrt{1-{\epsilon }^2}\sin \theta }{1+\epsilon \cos \theta }\end{array}}$

A partir de la primera relación entre cosE y cosθ obtenemos la ecuación de la trayectoria en función de E

${\displaystyle \begin{array}{l}r=\frac{a\left(1-{\epsilon }^2\right)}{1+\epsilon \cos \theta }\\ r=\frac{a\left(1-{\epsilon }^2\right)}{1+\epsilon \left(\frac{\epsilon -\cos E}{\epsilon \cos E-1}\right)}=a\left(1-\epsilon \cos E\right)\end{array}}$

Que es la ecuación de la trayectoria en términos del ángulo E, que ya hemos obtenido

Utilizando la relación trigonométrica de la tangente del ángulo mitad, llegamos a una relación entre los ángulos E y la posición angular θ

${\displaystyle \begin{array}{l}{\tan }^2\left(\frac{E}{2}\right)=\frac{1-\cos E}{1+\cos E}\\ {\tan }^2\left(\frac{E}{2}\right)=\frac{1-\epsilon }{1+\epsilon }\cdot \frac{1-\cos \theta }{1+\cos \theta }=\frac{1-\epsilon }{1+\epsilon }{\tan }^2\left(\frac{\theta }{2}\right)\\ \tan \left(\frac{E}{2}\right)=\sqrt{\frac{1-\epsilon }{1+\epsilon }}\tan \left(\frac{\theta }{2}\right)\\ E=2{\tan }^{-1}\left(\sqrt{\frac{1-\epsilon }{1+\epsilon }}\tan \left(\frac{\theta }{2}\right)\right)\end{array}}$

Que es la relación entre los ángulos θ y E obteneida anteriormente

Trayectorias parabólicas

Una trayectoria parabólica se caracteriza por que la energía de la partícula es nula y la excentricidad es ε=1. La ecuación de la trayectoria es

${\displaystyle r=\frac{d}{1+\cos \theta }}$

La posición de máximo acercamiento r_p al centro de fuerzas se produce para θ=0, por lo que r_p=d/2

La posición angular θ en función del tiempo t es

${\displaystyle \underset{0}{\overset{\theta }{{\displaystyle \int }}}\frac{d\theta }{{\left(1+\cos \theta \right)}^2}=\frac{G^2{M}^2{m}^3}{L^3}t}$

Haciendo estos cambios la integral es inmediata

${\displaystyle \begin{array}{l}\cos \theta =\frac{1-x^2}{1+x^2}\text{,\hspace{1em}}\sin \theta =\frac{2x}{1+x^2}\text{,\hspace{1em}}d\theta =\frac{2}{1+x^2}dx\text{,\hspace{1em}}x=\tan \left(\frac{\theta }{2}\right)\\ {{\displaystyle \int }}^{\text{}}\frac{d\theta }{{\left(1+\cos \theta \right)}^2}=\frac{1}{2}{{\displaystyle \int }}^{\text{}}(1+x^2)dx=\frac{1}{2}\left(x+\frac{1}{3}{x}^3\right)=\frac{1}{2}\tan \left(\frac{\theta }{2}\right)+\frac{1}{6}{\tan }^3\left(\frac{\theta }{2}\right)\end{array}}$

La ecuación de Kepler para órbitas parabólicas se escribe

${\displaystyle {\tan }^3\left(\frac{\theta }{2}\right)+3\tan \left(\frac{\theta }{2}\right)=6\frac{G^2{M}^2{m}^3}{L^3}t}$

Dada la posición θ del cuerpo celeste se calcula el tiempo t. Sin embargo, dado el tiempo t requiere algo más de esfuerzo calcular la posición angular θ como veremos en este ejemplo

Dado el tiempo t determinar la posición angular θ del satélite

Un satélite artificial describe una trayectoria parabólica, su velocidad en el perigeo es de v_p=10000 m/s. Determinar la distancia del satélite al centro de la Tierra seis horas más tarde.

Dado v_p, aplicamos el principio de conservación de la energía para calcular la distancia más próxima al centro de la Tierra, r_p. Para una trayectoria parabólica, la energía total es 0.

${\displaystyle \begin{array}{l}\frac{1}{2}m{v}_p^2-G\frac{Mm}{r_p}=0\\ r_p=\frac{2GM}{v_p^2}\end{array}}$

El momento angular vale L=m·r_p·v_p

Dado el instante t, el ángulo θ es la raíz real de la ecuación cúbica

${\displaystyle {\tan }^3\left(\frac{\theta }{2}\right)+3\tan \left(\frac{\theta }{2}\right)-6\frac{G^2{M}^2{m}^3}{L^3}t=0}$

Utilizamos la función raices_3, para obtener las tres raíces de la ecuación cúbica, una real y dos complejas conjugadas. Conocido el ángulo θ, obtenemos la distancia r al centro de la Tierra mediante la ecuación de la trayectoria

Elaboramos un script para realizar los cálculos

M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G
vp=10000; %velocidad en el perigeo
rp=2*G*M/vp^2; %perigeo
L=rp*vp; %momento angular
 
t=6*60*60; %instante 6 h
%resolvemos la ecuación cúbica
 
c=6*G^2*M^2*t/L^3;
x=raices_3([1,0,3,-c]);
th=0;
for i=1:3
    if isreal(x(i))
        th=2*atan(x(i));
        break;
    end
end
r=2*rp/(1+cos(th));
disp(r)

La distancia del satálite al centro de la Tierra en el instante t=6 h, en km, es

  8.6993e+04

Trayectorias hiperbólicas

La energía total de una partícula que describe una trayectoria hiperbólica es positiva, la excentricidad es ε>1. La posición angular en función del tiempo es

${\displaystyle \underset{0}{\overset{\theta }{{\displaystyle \int }}}\frac{d\theta }{{(1+\epsilon \cos \theta )}^2}=\frac{G^2{M}^2{m}^3}{L^3}t}$

Calculamos la integral, sabiendo que ε>1, para ello seguimos los mismos pasos que para la trayectoria elíptica

${\displaystyle \begin{array}{l}\cos \theta =\frac{1-x^2}{1+x^2}\text{,\hspace{1em}}\sin \theta =\frac{2x}{1+x^2}\text{,\hspace{1em}}d\theta =\frac{2}{1+x^2}dx\text{,\hspace{1em}}x=\tan \left(\frac{\theta }{2}\right)\\ {{\displaystyle \int }}^{\text{}}\frac{d\theta }{{\left(1+\epsilon \cos \theta \right)}^2}=\frac{2}{{\left(1+\epsilon \right)}^2}{{\displaystyle \int }}^{\text{}}\frac{1+x^2}{{\left(1-a{x}^2\right)}^2}dx=\\ \frac{2}{{\left(1+\epsilon \right)}^2}\left\{-\frac{1}{a}{{\displaystyle \int }}^{\text{}}\frac{dx}{1-a{x}^2}+\frac{a+1}{a}{{\displaystyle \int }}^{\text{}}\frac{dx}{{\left(1-a{x}^2\right)}^2}\right\}\end{array}}$

con a=(ε-1)/(ε+1)

La primera integral es inmediata

${\displaystyle {{\displaystyle \int }}^{\text{}}\frac{dx}{1-a{x}^2}=\frac{1}{\sqrt{a}}{{\displaystyle \int }}^{\text{}}\frac{du}{1-u^2}=\frac{1}{\sqrt{a}}\text{arctanh}u=\frac{1}{\sqrt{a}}\text{arctanh}\left(\sqrt{a}x\right)}$

La segunda integral es algo más laboriosa

${\displaystyle {{\displaystyle \int }}^{\text{}}\frac{dx}{{\left(1-a{x}^2\right)}^2}=\frac{1}{\sqrt{a}}{{\displaystyle \int }}^{\text{}}\frac{du}{{\left(1-u^2\right)}^2}}$

Haciendo el cambio, u=tanhv, du=dv/cosh²(v),

${\displaystyle {{\displaystyle \int }}^{\text{}}\frac{du}{{\left(1-u^2\right)}^2}={\displaystyle \int {\cosh }^{2}v·dv=}{\displaystyle \int \frac{\exp (2v)+\exp (-2v)+2}{4}·dv=}\frac{1}{2}v+\frac{1}{4}\sinh \left(2v\right)}$

Deshaciendo los cambios

${\displaystyle \begin{array}{l}{{\displaystyle \int }}^{\text{}}\frac{dx}{{\left(1-a{x}^2\right)}^2}=\frac{1}{2\sqrt{a}}\left(\arctan \left(\sqrt{a}x\right)+\frac{1}{2}\sinh \left(2\text{arctanh}\left(\sqrt{a}x\right)\right)\right)=\\ \frac{1}{2\sqrt{a}}\left(\text{arctanh}\left(\sqrt{a}x\right)+\frac{\sqrt{a}x}{1-a{x}^2}\right)\end{array}}$

Se ha utilizado la relación

${\displaystyle \sinh \left(2u\right)=\frac{2\tanh (u/2)}{1-{\tanh }^2(u/2)}}$

El resultado de la suma de las dos integrales es

${\displaystyle \begin{array}{l}{{\displaystyle \int }}^{\text{}}\frac{d\theta }{{\left(1+\epsilon \cos \theta \right)}^2}=\frac{2}{{\left(1+\epsilon \right)}^2}\left\{-\frac{1}{a}{{\displaystyle \int }}^{\text{}}\frac{dx}{1-a{x}^2}+\frac{a+1}{a}{{\displaystyle \int }}^{\text{}}\frac{dx}{{\left(1-a{x}^2\right)}^2}\right\}=\\ \frac{1}{{\left(1+\epsilon \right)}^2}\left\{\frac{a-1}{a\sqrt{a}}\text{arctanh}\left(\sqrt{a}x\right)+\frac{a+1}{2a\sqrt{a}}\frac{2\sqrt{a}x}{1-a{x}^2}\right\}\end{array}}$

Definimos el ángulo F

${\displaystyle \tanh \left(\frac{F}{2}\right)=\sqrt{a}x}$

La suma de las dos integrales se convierte en

${\displaystyle \frac{1}{{\left(\epsilon +1\right)}^{2}a\sqrt{a}}\left\{\left(a-1\right)\frac{F}{2}+\frac{\left(a+1\right)}{2}\sinh F\right\}=\frac{1}{{\left({\epsilon }^2-1\right)}^{3/2}}\left(-F+\epsilon \sinh F\right)}$

La ecuación de Kepler para trayectorias hipérbólicas es

${\displaystyle \begin{array}{l}\underset{0}{\overset{\theta }{{\displaystyle \int }}}\frac{d\theta }{{(1+\epsilon \cos \theta )}^2}=\frac{G^2{M}^2{m}^3}{L^3}t\\ \epsilon \sinh F-F={\left({\epsilon }^2-1\right)}^{3/2}\frac{G^2{M}^2{m}^3}{L^3}t\end{array}}$

Donde se ha definido el ángulo F

${\displaystyle \tanh \left(\frac{F}{2}\right)=\sqrt{\frac{\epsilon -1}{\epsilon +1}}\tan \left(\frac{\theta }{2}\right)}$

Para una partícula que viene del infinito con velocidad v₀ y parámetro de impacto b, el momento angular L y la energía E valen

${\displaystyle \begin{array}{l}L=mb{v}_0\\ E=\frac{1}{2}m{v}_0^2\end{array}}$

La ecuación de Kepler se escribe

${\displaystyle t=\frac{b}{p{v}_0}\left\{\sqrt{1+p^2}\sinh (F)-F\right\}\text{,\hspace{1em}}p=\sqrt{{\epsilon }^2-1}=\frac{b{v}_0^2}{GM}}$

Ecuación de la trayectoria en términos del ángulo F

La ecuación de la trayectoria es

${\displaystyle r=\frac{a\left({\epsilon }^2-1\right)}{1+\epsilon \cos \theta }}$

la excentricidad ε>1. Conociendo la relación entre los ángulos θ y F

${\displaystyle \tan \left(\frac{\theta }{2}\right)=\sqrt{\frac{\epsilon +1}{\epsilon -1}}\tanh \left(\frac{F}{2}\right)}$

Determinaremos la ecuación de la trayectoria en términos del ángulo E, utilizando las relaciones trigonométricas

${\displaystyle {\cos }^2\theta =\frac{1}{1+{\tan }^2\theta },\tan \theta =\frac{2\tan \left(\frac{\theta }{2}\right)}{1-{\tan }^2\left(\frac{\theta }{2}\right)}}$

El resultado es

${\displaystyle \cos \theta =\frac{1-\frac{\epsilon +1}{\epsilon -1}{\tanh }^2\left(\frac{F}{2}\right)}{1+\frac{\epsilon +1}{\epsilon -1}{\tanh }^2\left(\frac{F}{2}\right)}=\frac{\left(\epsilon -1\right){\cosh }^2\left(\frac{F}{2}\right)-\left(\epsilon +1\right){\sinh }^2\left(\frac{F}{2}\right)}{\left(\epsilon -1\right){\cosh }^2\left(\frac{F}{2}\right)+\left(\epsilon +1\right){\sinh }^2\left(\frac{F}{2}\right)}=\frac{\epsilon -\cosh F}{\epsilon \cosh F-1}}$

Hemos utilizado las relaciones

${\displaystyle \begin{array}{l}{\cosh }^2\left(\frac{\theta }{2}\right)-{\sinh }^2\left(\frac{\theta }{2}\right)=1\\ \cosh \theta ={\cosh }^2\left(\frac{\theta }{2}\right)+{\sinh }^2\left(\frac{\theta }{2}\right)\end{array}}$

Finalmente, la ecuación de la trayectoria es

${\displaystyle r=\frac{a\left({\epsilon }^2-1\right)}{1+\epsilon \frac{\epsilon -\cosh F}{\epsilon \cosh F-1}}=a\left(\epsilon \cosh F-1\right)}$

Dada la posición angular θ del satélite determinar el tiempo t

Un satélite artificial tiene un perigeo de 300 km de altura sobre la superficie de la Tierra, llevando una velocidad de 15000 m/s. Calcular la distancia al centro de la Tierra cuando su posición angular es θ=100°.

La ecuación de la hipérbola es

${\displaystyle r=\frac{d}{1+\epsilon \cos \theta }}$

La energía y el momento angular valen

${\displaystyle \begin{array}{l}E=\frac{1}{2}m{v}_p^2-G\frac{Mm}{r_p}\\ L=m{r}_p{v}_p\end{array}}$

que nos permiten calcular el parámetro d y la excentricidad ε

${\displaystyle \epsilon =\sqrt{1+\frac{2L^{2}E}{m^3{G}^2{M}^2}}\text{,\hspace{1em}}d=\frac{L^2}{GM{m}^2}}$

M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G
R=6370000; %radio de la Tierra
vp=15000; %velocidad en el perigeo
rp=300000+R; %perigeo

L=rp*vp; %momento angular
E=vp^2/2-G*M/rp; %energía
d=L^2/(G*M);
e=sqrt(1+2*L^2*E/(G*M)^2);

La excentricidad ε>1, la trayectoria es hiperbólica

>> e
e =    2.7625

Calculamos el mayor ángulo θ posible, el ángulo θ_∞ de la asíntota

>> acos(-1/e)*180/pi
ans =  111.2222

que es un ángulo mayor que θ=100°

Calculamos el ángulo F

${\displaystyle \tanh \left(\frac{F}{2}\right)=\sqrt{\frac{\epsilon -1}{\epsilon +1}}\tan \left(\frac{\theta }{2}\right)}$

La ecuación de Kepler para trayectorias hiperbólicas nos proporciona el instante t en el que el satélite alcanza la posición θ

${\displaystyle \epsilon \sinh F-F={\left({\epsilon }^2-1\right)}^{3/2}\frac{G^2{M}^2{m}^3}{L^3}t}$

M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G
R=6370000; %radio d ela Tierra
vp=15000; %velocidad en el perigeo
rp=300000+R; %perigeo

L=rp*vp; %momento angular
E=vp^2/2-G*M/rp; %energía
d=L^2/(G*M);
e=sqrt(1+2*L^2*E/(G*M)^2);

th=100*pi/180; %posición angular
F=2*atanh(sqrt((e-1)/(e+1))*tan(th/2));
t=(e*sinh(F)-F)*L^3/((e^2-1)^(3/2)*(G*M)^2);
disp(t/60)

El instante t en minutos es

   68.6725

Dada el instante t determinar la posición angular θ del satélite

Completamos este ejemplo, calculando la posición y velocidad del satélite 3 h después de haber pasado por la posición θ=100°, es decir en el instante t=14 920 s

Conocido t, la ecuación de Kepler despejamos F en la ecuación transcendente

${\displaystyle \epsilon \sinh F-F={\left({\epsilon }^2-1\right)}^{3/2}\frac{G^2{M}^2{m}^3}{L^3}t}$

Conocido F despejamos la posición angular θ

${\displaystyle \tanh \left(\frac{F}{2}\right)=\sqrt{\frac{\epsilon -1}{\epsilon +1}}\tan \left(\frac{\theta }{2}\right)}$

M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G
R=6370000; %radio de la Tierra
vp=15000; %velocidad en el perigeo
rp=300000+R; %perigeo

L=rp*vp; %momento angular
E=vp^2/2-G*M/rp; %energía
d=L^2/(G*M);
e=sqrt(1+2*L^2*E/(G*M)^2);

th=100*pi/180; %posición angular
F=2*atanh(sqrt((e-1)/(e+1))*tan(th/2));
t=(e*sinh(F)-F)*L^3/((e^2-1)^(3/2)*(G*M)^2);

t=t+3*3600;
k=(e^2-1)^(3/2)*(G*M)^2*t/L^3;
f=@(x) e*sinh(x)-x-k;
F=fzero(f,[0,5]);
th=2*atan(tanh(F/2)*sqrt((e+1)/(e-1))); %posición angular
r=d/(1+e*cos(th)); %distancia al centro de la Tierra
fprintf('Posición angular %3.1f, distancia al centro %3.1f km\n',
th*180/pi,r/1000)

Posición angular 107.8, distancia al centro 162819.7 km

Para conocer aproximadamente la raíz de la ecuación transcendente, representamos la función f(x)=esinh(x)-x-k

>> fplot(f,[0,5])
>> grid on

Aproximadamente, F=3.5 rad

Calculamos las componentes de la velocidad a la distancia r=162819.7 km del centro de la Tierra

${\displaystyle \begin{array}{l}L=r\left(r\frac{d\theta }{dt}\right)=r{v}_{\theta }\\ E=\frac{1}{2}m\left(v_r^2+v_{\theta }^2\right)-G\frac{Mm}{r}\end{array}}$

Cuando la partícula llega al infinito, la energía potencial es nula,

${\displaystyle E=\frac{1}{2}m{v}_{\infty }^2}$

Calculamos las componentes de la velocidad v_r y v_θ, el módulo de la velocidad v y la comparamos con la velocidad del satélite muy alejado de la Tierra, v_∞

>> v_th=L/r
v_th =  614.4836
>> v_r=sqrt(2*(E+G*M/r)-v_th^2)
v_r =   1.0484e+04
>> sqrt(v_r^2+v_th^2)
ans =   1.0502e+04
>> v_inf=sqrt(2*E)
v_inf =   1.0266e+04

Completamos el script, trazando la trayectoria hiperbólica y las posiciones del satélite en dos instantes: t₁=4 120.3 s y t₂=14 920 s

M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G
R=6370000; %radio de la Tierra
vp=15000; %velocidad en el perigeo
rp=300000+R; %perigeo

L=rp*vp; %momento angular
E=vp^2/2-G*M/rp; %energía
d=L^2/(G*M);
e=sqrt(1+2*L^2*E/(G*M)^2);

th_1=100*pi/180; %posición angular
F=2*atanh(sqrt((e-1)/(e+1))*tan(th_1/2));
t1=(e*sinh(F)-F)*L^3/((e^2-1)^(3/2)*(G*M)^2);

t2=t1+3*3600;
k=(e^2-1)^(3/2)*(G*M)^2*t2/L^3;
f=@(x) e*sinh(x)-x-k;
F=fzero(f,[0,5]);
th_2=2*atan(tanh(F/2)*sqrt((e+1)/(e-1))); %posición angular
r=d/(1+e*cos(th_2)); %distancia al centro de la Tierra

hold on
f=@(x) d./(1+e*cos(x));
fplot(@(x) f(x).*cos(x), @(x) f(x).*sin(x),[0,acos(-1/e)])
plot(f(th_1)*cos(th_1),f(th_1)*sin(th_1),'bo','markersize',3,'markeredgecolor'
,'b','markerfacecolor','b')
line([0,f(th_1)*cos(th_1)],[0,f(th_1)*sin(th_1)],'lineStyle','--','color','b')
plot(f(th_2)*cos(th_2),f(th_2)*sin(th_2),'ro','markersize',3,'markeredgecolor'
,'r','markerfacecolor','r')
line([0,f(th_2)*cos(th_2)],[0,f(th_2)*sin(th_2)],'lineStyle','--','color','r')
grid on
xlabel('x')
ylabel('y')
title('Trayectoria hiperbólica')

Dada la posición radial r del satélite determinar el tiempo t

Expresamos la energía del cuerpo celeste de masa m en términos de la distancia radial r al centro de fuerzas

${\displaystyle E=\frac{1}{2}m{\left(\frac{dr}{dt}\right)}^2+\frac{L^2}{2m{r}^2}-\frac{GMm}{r}}$

Para un cuerpo celeste que describe una trayectoria hiperbólica, E>0, es la energía cinética en el infinito (la potencial es nula)

${\displaystyle E=\frac{1}{2}m{v}_{\infty }^2}$

Establecemos el tiempo t=0, cuando el cuerpo celeste pasa por la posición más cercana al centro de fuerzas, a una distancia r_p. La componente radial de la velocidad en esta posición es nula, dr/dt=0

${\displaystyle \begin{array}{l}\frac{1}{2}m{v}_{\infty }^2=\frac{L^2}{2m{r}_p^2}-\frac{GMm}{r_p}\\ L^2=m^2{r}_p^2{v}_{\infty }^2+2GM{m}^2{r}_p\end{array}}$

En la expresión de la energía despejamos dt e integramos

${\displaystyle \begin{array}{l}\frac{dr}{\sqrt{\frac{2}{m}\left(E-\frac{L^2}{2m{r}^2}+\frac{GMm}{r}\right)}}=dt\\ t=\sqrt{\frac{m}{2E}}{\displaystyle \underset{r_p}{\overset{r}{\int }}\frac{r·dr}{\sqrt{r^2+\frac{GMm}{E}r-\frac{L^2}{2mE}}}}\end{array}}$

Resolvemos la integral

${\displaystyle {\displaystyle \int \frac{x·dx}{\sqrt{x^2+bx-c}}},b=\frac{GMm}{E},c=\frac{L^2}{2mE}}$

Descomponemos la integral en suma de dos

${\displaystyle {\displaystyle \int \frac{x·dx}{\sqrt{x^2+bx-c}}}=\frac{1}{2}{\displaystyle \int \frac{(2x+b)·dx}{\sqrt{x^2+bx-c}}}-\frac{b}{2}{\displaystyle \int \frac{dx}{\sqrt{x^2+bx-c}}}}$

La primera es inmediata y en la segunda, completamos cuadrados

${\displaystyle \sqrt{x^2+bx-c}-\frac{b}{2}{\displaystyle \int \frac{dx}{\sqrt{{\left(x+\frac{b}{2}\right)}^2-\left(c+\frac{b^2}{4}\right)}}}=\sqrt{x^2+bx-c}-\frac{b}{2}{\cosh }^{-1}\left(\frac{x+\frac{b}{2}}{\sqrt{c+\frac{b^2}{4}}}\right)}$

Ahora resolvemos la integral definida entre r_p posición más cercana y r

${\displaystyle \begin{array}{l}t=\sqrt{\frac{m}{2E}}{\left(\sqrt{r^2+\frac{GMm}{E}r-\frac{L^2}{2mE}}-\frac{GMm}{2E}{\cosh }^{-1}\left(\frac{r+\frac{GMm}{2E}}{\sqrt{\frac{L^2}{2mE}+\frac{G^2{M}^2{m}^2}{4E^2}}}\right)\right)}_{r_p}^r=\\ \frac{1}{v_{\infty }}{\left(\sqrt{r^2-r_p^2+\frac{2GM}{v_{\infty }^2}(r-r_p)}-\frac{GM}{v_{\infty }^2}{\cosh }^{-1}\left(\frac{r+\frac{GM}{v_{\infty }^2}}{r_p+\frac{GM}{v_{\infty }^2}}\right)\right)}_{r_p}^r=\\ \frac{1}{v_{\infty }}\left(\sqrt{r^2-r_p^2+\frac{2GM}{v_{\infty }^2}(r-r_p)}-\frac{GM}{v_{\infty }^2}{\cosh }^{-1}\left(\frac{r+\frac{GM}{v_{\infty }^2}}{r_p+\frac{GM}{v_{\infty }^2}}\right)\right)\end{array}}$

El resultado final es

${\displaystyle t=\frac{GM}{v_{\infty }^3}\left(\sqrt{{\left(1+\frac{r{v}_{\infty }^2}{GM}\right)}^2-{\left(1+\frac{r_p{v}_{\infty }^2}{GM}\right)}^2}-{\cosh }^{-1}\left(\frac{1+\frac{r{v}_{\infty }^2}{GM}}{1+\frac{r_p{v}_{\infty }^2}{GM}}\right)\right)}$

M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G
R=6370000; %radio de la Tierra
vp=15000; %velocidad en el perigeo
rp=300000+R; %perigeo

L=rp*vp; %momento angular
E=vp^2/2-G*M/rp; %energía
d=L^2/(G*M);
e=sqrt(1+2*L^2*E/(G*M)^2);
th=100*pi/180; %posición angular
r=d/(1+e*cos(th)); %distancia al centro de la Tierra

F=2*atanh(sqrt((e-1)/(e+1))*tan(th/2)); 
t1=(e*sinh(F)-F)*L^3/((e^2-1)^(3/2)*(G*M)^2);

vi=sqrt(2*E); %velocidad en el infinito
t2=(G*M/vi^3)*(sqrt((1+r*vi^2/(G*M))^2-(1+rp*vi^2/(G*M))^2)-
acosh((1+r*vi^2/(G*M))/(1+rp*vi^2/(G*M))));
disp([t1,t2])

Los tiempos calculados para la posición angular θ=100°, o para la distancia radial r=4.8235e·10⁷ m son los mismos 4 120.3 s

   1.0e+03 *    4.1203    4.1203

Referencias

Howard D. Curtis. Orbital Mechanics for Engineering Students. Elsevier, capítulo 3