Órbita de un Planeta

[Graphics:Kepler/Saturno.GIF]

Cálculo de la posición

A partir de las tres leyes de Kepler vamos a calcular la ecuación que determina la posición de un planeta. El resultado está contenido en el paquete Kepler.m que se describe al final del documento. La idea de este trabajo es mostrar cómo se pueden encontrar las órbitas planetarias, utilizando solamente técnicas de Cálculo Diferencial e Integral, sin necesidad de resolver  las ecuaciones del movimiento determinadas por la Ley de Gravitación de Newton.

Leyes de Kepler

La primera ley asegura que la órbita de un planeta es una elipse, encontrándose el Sol en uno de sus focos.  La segunda ley  nos dice que el área barrida por el vector de posición con centro el Sol, es proporcional al tiempo transcurrido. La tercera ley determina  que el cuadrado del periodo (P) es proporcional al cubo del semieje de la elipse (a).

[Graphics:Kepler/Kepler_gr_2.gif]

Parametrización

Calculamos primero la ecuación de la elipse, en coordenadas polares, con centro el Sol:
Si en coordenadas rectangulares, la ecuación viene dada por

[Graphics:Kepler/Kepler_gr_2.gif]

sabiendo que si  r y s son las distancias del punto (x,y) a cada uno de los dos focos, entonces s+r = 2a, y utilizando la fórmula del coseno

s 2 = r 2+ 4 c 2 - 4 r c Cos[z],

substituyendo obtenemos:

4 a 2 + r2 - 4 a r = r2 + 4 c2 - 4 r c Cos[z].

De aquí podemos despejar r como función del ángulo z:

[Graphics:Kepler/Kepler_gr_2.gif]

lo que nos da la ecuación en coordenadas polares.                                               

[Graphics:Kepler/Kepler_gr_3.gif]

Solución del problema

Nuestro objetivo es encontrar el ángulo que determina la posición del planeta sobre la elipse, en función del tiempo: z[t].  Sea A[z0,z1] el área del sector de la elipse, determinada por los radio-vectores que tienen como ángulos los valores z0 y z1.

[Graphics:Kepler/Kepler_gr_4.gif]


Es fácil probar que A[z0,z1] = Integrate[r[z]2 /2,{z,z0,z1}]. Así, la segunda ley de Kepler, con las condiciones iniciales  z[0] = 0  y  z[P] = 2π,  se traduce en la igualdad

                                       A[0,z[t]] = A[z[s],z[t+s]],
                                       
para cualquier valor de s y t. Por lo tanto, si derivamos con respecto a s, y utilizando el Teorema Fundamental del Cálculo, obtenemos

                                r2 [z[t+s]] z'[t+s] = r2 [z[s]] z'[s].

Haciendo s = 0 y observando que r[z[0]] = r[0] = a + c, concluimos que

[Graphics:Kepler/Kepler_gr_4.gif]

Por lo tanto, integrando entre 0 y t, y usando el cambio de variables  u = z[s],  obtenemos

[Graphics:Kepler/Kepler_gr_4.gif]

Hagamos ahora el cálculo de la integral que aparece a la izquierda de la ecuación anterior.


Integrate[1/(a-c Cos[u])^2,u]

y obtenemosla solución  que buscamos:

[Graphics:Kepler/Kepler_gr_2.gif]

Ejemplos


Orbita[a,c,P,n] dibuja la órbita (con n puntos), de un planeta con semieje mayor=a, excentricidad=c  y periodo=P. Por defecto, a=1, P=1 y n=30. La función Orbita se define al final del documento.


Orbita[1,.7,1,30]

[Graphics:Kepler/Kepler_gr_5.gif]

     -Graphics-


BeginPackage["Kepler`"]

Orbita::usage="Orbita[a,c,P,n] dibuja la orbita (con
n puntos), de un planeta con semieje mayor=a,  
excentricidad=c y periodo=P. Por defecto, a=1, P=1 y
n=30."    

Begin["`Private`"]

(* La funcion F corresponde al valor:     *)
(* F=Integrate[1/(A-C Cos[s])^2,{s,0,x}] *)

F=(-2*A*ArcTanh[((A + C)^(1/2)*Tan[x/2])/
       (-A + C)^(1/2)])/
   ((-A + C)^(1/2)*(A + C)^(1/2)*(A^2 - C^2)) +
  (C*Sin[x])/((-A^2 + C^2)*(-A + C*Cos[x]));
  
Off[FindRoot::precw,FindRoot::cvnwt,Plot::plnr];

sol=Graphics[{PointSize[.02],Point[{0,0}]}];
  

Orbita[a_:1,c_,P_:1, n_:30]:=Module[
{b,K,G,Theta,Angulo,lista,radios,datos,ultimo,orb},

b=Sqrt[a^2-c^2];

K=2 Pi a/(b^3 P);

G=Evaluate[Chop[F/.{A->a,C->c}]];

Theta[t_]:=FindRoot[G==K t,{x,Pi/2}];

Angulo[t_]:=N[Which[0<=t <=P/2,x/.Theta[t],
    P/2<t<P,2Pi-(x/.Theta[P-t])]];
    
lista= Table[Angulo[P k/n],{k,0,n-1}];

radios=Map[(b^2/(a-c Cos[#]))&,lista];

datos=
Transpose[{Map[Cos,lista], Map[Sin,lista]}];

ultimo=radios datos;

orb=ListPlot[ultimo,AspectRatio->Automatic,
DisplayFunction-> Identity,Axes->None];

Show[{orb,sol},DisplayFunction->$DisplayFunction]
];

End[];

Protect[Orbita];
EndPackage[]    ;             
    


Converted by Mathematica      March 1, 2002