Órbita de un Planeta

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

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

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

substituyendo obtenemos:

4 a ² + r² - 4 a r = r² + 4 c² - 4 r c Cos[z].

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

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.

Es fácil probar que A[z0,z1] = Integrate[r[z]² /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

                                r² [z[t+s]] z'[t+s] = r² [z[s]] z'[s].

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

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

y obtenemosla solución  que buscamos:

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-

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[]    ;