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.
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).
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 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:
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 /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
Integrate[1/(a-c Cos[u])^2,u]
y obtenemosla solución que buscamos:
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[] ;