El problema de Kepler

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.

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

r= d 1+εcosθ ε= 1+ 2 L 2 E m 3 G 2 M 2 d= L 2 GM m 2

La distancia más cercana al Sol r1 se obtiene en la ecuación de la trayectoria con θ=0 , y la más alejada r2 con θ=π. Como vemos en la figura r1+r2=2a

r 1 =( L 2 GM m 2 ) 1 1+ε r 2 =( L 2 GM m 2 ) 1 1ε 2a= r 1 + r 2 a=( L 2 GM m 2 ) 1 1 ε 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

L=m r 2 dθ dt

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

0 θ dθ ( 1+εcosθ ) 2 = G 2 M 2 m 3 L 3 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 de la izquierda

cosθ= 1 x 2 1+ x 2 sinθ= 2x 1+ x 2 dθ= 2 1+ x 2 dxx=tan( θ 2 ) dθ ( 1+εcosθ ) 2 = 2 ( 1+ε ) 2 1+ x 2 ( 1+a x 2 ) 2 dx= 2 (1+ε) 2 { 1 a dx 1+a x 2 + a1 a dx ( 1+a x 2 ) 2 }

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

La primera integral es inmediata

dx 1+a x 2 = 1 a du 1+ u 2 = 1 a arctanu= 1 a arctan( a x )

La segunda integral es algo más laboriosa

dx ( 1+a x 2 ) 2 = 1 a du ( 1+ u 2 ) 2

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

du ( 1+ u 2 ) 2 = cos 2 v·dv= 1+cos( 2v ) 2 ·dv= 1 2 v+ 1 4 sin( 2v )

Deshaciendo los cambios

dx ( 1+a x 2 ) 2 = 1 2 a ( arctan( a x )+ 1 2 sin( 2arctan( a x ) ) )= 1 2 a ( arctan( a x )+ a x 1+a x 2 )

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

sin( 2u )= 2tanu 1+ tan 2 u

El resultado de la suma de las dos integrales es

dθ ( 1+εcosθ ) 2 = 2 (1+ε) 2 { 1 a dx 1+a x 2 + a1 a dx ( 1+a x 2 ) 2 }= 1 ( 1+ε ) 2 { a+1 a a arctan( a x )+ a1 2a a 2 a x 1+a x 2 }=

Definimos el ángulo E tal que

tan( E 2 )= a x

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

1 ( 1+ε ) 2 a a { ( a+1 ) E 2 + ( a1 ) 2 sinE }= 1 ( 1 ε 2 ) 3/2 ( E+εsinE )

La ecuación de Kepler se escribe

0 θ dθ (1+εcosθ) 2 = G 2 M 2 m 3 L 3 t E+εsinE= ( 1 ε 2 ) 3/2 G 2 M 2 m 3 L 3 t

Donde el ángulo E se ha definido como

tan( E 2 )= 1ε 1+ε tan( θ 2 )

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

P 2 =4 π 2 a 3 GM P=2π L 3 G 2 M 2 m 3 1 (1 ε 2 ) 3/2

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

E+εsinE=2π t P

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

Representamos Me 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')

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

Teniendo en cuenta las relaciones trigonométricas

sinE= 2tan( E 2 ) 1+ tan 2 ( E 2 ) cosE= 1 tan 2 ( E 2 ) 1+ tan 2 ( E 2 )

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

sinE= 1 ε 2 sinθ 1+εcosθ cosE= ε+cosθ 1+εcosθ

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 Me

M e =EεsinE

Conocido Me, obtenemos el tiempo

t= M e 2π P P 2 = 4 π 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 r1=9.6·106 m y r2=21·106. 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 Me en la ecuación de Kepler, obtenemos una ecuación trascendente cuya raíz tenemos que determinar por procedimientos numéricos

EεsinE= 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

x i+1 = x i f( x i ) f'( x i ) E i+1 = E i E i εsin E i M e 1εcos E i

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

tan( θ 2 )= 1+ε 1ε tan( E 2 )

Ejemplo

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

M e =2π t P M e =2π 10800 18828 =3.60

Resolvemos la ecuación trascendente E-ε·sinE=Me por el procedimiento de Newton tomando como valor inicial de partida E0=Me/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

tan( θ 2 )= 1+0.37 10.37 tan( 3.48 2 )θ=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 las siguientes relaciones geométricas: acosE=c+rcosθ

acosE=aε+ a( 1 ε 2 ) 1+εcosθ cosθ cosE= ε+cosθ 1+εcosθ sin 2 E+ cos 2 E=1sinE= 1 ε 2 sinθ 1+εcosθ

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

r= a( 1 ε 2 ) 1+εcosθ r= a( 1 ε 2 ) 1+ε( εcosE εcosE1 ) =a( 1εcosE )

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 θ

tan 2 ( E 2 )= 1cosE 1+cosE tan 2 ( E 2 )= 1ε 1+ε 1cosθ 1+cosθ = 1ε 1+ε tan 2 ( θ 2 ) tan( E 2 )= 1ε 1+ε tan( θ 2 ) E=2 tan 1 ( 1ε 1+ε tan( θ 2 ) )

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

r= d 1+cosθ

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

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

0 θ dθ ( 1+cosθ ) 2 = G 2 M 2 m 3 L 3 t

Haciendo estos cambios la integral es inmediata

cosθ= 1 x 2 1+ x 2 sinθ= 2x 1+ x 2 dθ= 2 1+ x 2 dxx=tan( θ 2 ) dθ ( 1+cosθ ) 2 = 1 2 (1+ x 2 )dx= 1 2 ( x+ 1 3 x 3 )= 1 2 tan( θ 2 )+ 1 6 tan 3 ( θ 2 )

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

tan 3 ( θ 2 )+3tan( θ 2 )=6 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 vp=10000 m/s. Determinar la distancia del satélite al centro de la Tierra seis horas más tarde.

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

1 2 m v p 2 G Mm r p =0 r p = 2GM v p 2

El momento angular vale L=m·rp·vp

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

tan 3 ( θ 2 )+3tan( θ 2 )6 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

0 θ dθ (1+εcosθ) 2 = 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

cosθ= 1 x 2 1+ x 2 sinθ= 2x 1+ x 2 dθ= 2 1+ x 2 dxx=tan( θ 2 ) dθ ( 1+εcosθ ) 2 = 2 ( 1+ε ) 2 1+ x 2 ( 1a x 2 ) 2 dx= 2 ( 1+ε ) 2 { 1 a dx 1a x 2 + a+1 a dx ( 1a x 2 ) 2 }

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

La primera integral es inmediata

dx 1a x 2 = 1 a du 1 u 2 = 1 a arctanhu= 1 a arctanh( a x )

La segunda integral es algo más laboriosa

dx ( 1a x 2 ) 2 = 1 a du ( 1 u 2 ) 2

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

du ( 1 u 2 ) 2 = cosh 2 v·dv= exp(2v)+exp(2v)+2 4 ·dv= 1 2 v+ 1 4 sinh( 2v )

Deshaciendo los cambios

dx ( 1a x 2 ) 2 = 1 2 a ( arctan( a x )+ 1 2 sinh( 2arctanh( a x ) ) )= 1 2 a ( arctanh( a x )+ a x 1a x 2 )

Se ha utilizado la relación

sinh( 2u )= 2tanh(u/ 2) 1 tanh 2 (u/ 2)

El resultado de la suma de las dos integrales es

dθ ( 1+εcosθ ) 2 = 2 ( 1+ε ) 2 { 1 a dx 1a x 2 + a+1 a dx ( 1a x 2 ) 2 }= 1 ( 1+ε ) 2 { a1 a a arctanh( a x )+ a+1 2a a 2 a x 1a x 2 }

Definimos el ángulo F

tanh( F 2 )= a x

La suma de las dos integrales se convierte en

1 ( ε+1 ) 2 a a { ( a1 ) F 2 + ( a+1 ) 2 sinhF }= 1 ( ε 2 1 ) 3/2 ( F+εsinhF )

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

0 θ dθ (1+εcosθ) 2 = G 2 M 2 m 3 L 3 t εsinhFF= ( ε 2 1 ) 3/2 G 2 M 2 m 3 L 3 t

Donde se ha definido el ángulo F

tanh( F 2 )= ε1 ε+1 tan( θ 2 )

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

L=mb v 0 E= 1 2 m v 0 2

La ecuación de Kepler se escribe

t= b p v 0 { 1+ p 2 sinh(E)E }p= ε 2 1 = b v 0 2 GM

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

r= d 1+εcosθ

La energía y el momento angular valen

E= 1 2 m v p 2 G Mm r p L=m r p v p

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

ε= 1+ 2 L 2 E m 3 G 2 M 2 d= 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

tanh( F 2 )= ε1 ε+1 tan( θ 2 )

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

εsinhFF= ( ε 2 1 ) 3/2 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

εsinhFF= ( ε 2 1 ) 3/2 G 2 M 2 m 3 L 3 t

Conocido F despejamos la posición angular θ

tanh( F 2 )= ε1 ε+1 tan( θ 2 )

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

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 velcidad a la distancia r=162819.7 km del centro de la Tierra

L=r( r dθ dt )=r v θ E= 1 2 m( v r 2 + v θ 2 )G Mm r

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

E= 1 2 m v 2

Calculamos las componentes de la velocidad vr 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: t1=4 120.3 s y t2=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')

Referencias

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