Las leyes de Kepler

Las leyes de Kepler describen la cinemática del movimiento de los planetas en torno al Sol.

Primera ley

Los planetas describen órbitas elípticas estando el Sol en uno de sus focos

Segunda ley

El vector posición de cualquier planeta respecto del Sol, barre áreas iguales de la elipse en tiempos iguales.

La ley de las áreas es equivalente a la constancia del momento angular, es decir, cuando el planeta está más alejado del Sol (afelio) su velocidad es menor que cuando está más cercano al Sol (perihelio). En el afelio y en el perihelio, el momento angular L es el producto de la masa del planeta, por su velocidad y por su distancia al centro del Sol.

L=mr₁·v₁=mr₂·v₂

Tercera ley

Los cuadrados de los periodos P de revolución son proporcionales a los cubos de los semiejes mayores a de la elipse.

P²=k·a³

Como podemos apreciar, el periodo de los planetas depende solamente del eje mayor de la elipse. Los tres planetas de la animación tienen el mismo eje mayor 2a=6 unidades, por tanto, tienen el mismo periodo.

Geometría de la elipse

Como se ha descrito en la página Ecuación de la trayectoria otras páginas, una partícula sometida a una fuerza central atractiva, inversamente proporcional al cuadrado de la distancia al centro de fuerzas, describe una trayectoria elíptica si su energía total E<0 es negativa,

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

Los parámetros d y ε (excentricidad) están relacionados con el momento angular L y la energía E de la partícula

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

Una elipse es una cónica cuya excentricidad ε=c/a es menor que la unidad. Para que una partícula sometida a una fuerza central, atractiva, inversamente proporcional al cuadrado de las distancias al centro de fuerzas, describa dicha trayectoria tiene que tener una energía total negativa (E<0).

La posición más cercana al foco r₁ se obtiene cuando θ=0 y la posición más alejada r₂ se obtiene cuando θ=π. Es decir,

${\displaystyle r_1=\frac{d}{1+\epsilon }{r}_2=\frac{d}{1-\epsilon }}$

Los semiejes a y b de la elipse valen

${\displaystyle \begin{array}{l}2a=r_2+r_{1}a=\frac{d}{1-{\epsilon }^2}\\ a^2=b^2+c^{2}b=a\sqrt{1-{\epsilon }^2}\end{array}}$

El semieje mayor de la elipse a es independiente del momento angular L y solamente, depende de la energía total E. El semieje menor b depende del momento angular L y de la energía E

${\displaystyle a=-\frac{mGM}{2E}b=\frac{L\sqrt{a}}{m\sqrt{GM}}}$

Ecuación de la elipse en coordenadas rectangulares

Dada la ecuación de la elipse en coordendas polares, vamos a deducir la ecuación de la elipse en coordendas rectangulares, estableciendo el origen del sistema de referencia en el centro de la elipse.

Partimos de la ecuación de la trayectoria en coordenadas polares

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

Las coordenadas rectangulares (x,y) del cuerpo celeste son:

x-c=r·cosθ
y=r·sinθ

La ecuación de la trayectoria se transforma en

${\displaystyle \begin{array}{l}d=r+\epsilon r\cos \theta \\ d=\sqrt{{\left(x-c\right)}^2+y^2}+\epsilon (x-c)\end{array}}$

Elevando al cuadrado

${\displaystyle \begin{array}{l}d-\epsilon (x-c)=\sqrt{{\left(x-c\right)}^2+y^2}\\ d^2+{\epsilon }^2{\left(x-c\right)}^2-2d\epsilon (x-c)={\left(x-c\right)}^2+y^2\\ d^2+{\epsilon }^2{x}^2+{\epsilon }^2{c}^2-2{\epsilon }^{2}cx-2d\epsilon x+2d\epsilon c=x^2+c^2-2cx+y^2\end{array}}$

Teniendo en cuenta la relación entre el parámetro d y el semieje mayor de la elipse a, y la definición de excentricidad ε

${\displaystyle \begin{array}{l}d=a(1-{\epsilon }^2)\\ c=a\epsilon \end{array}}$

La expresión anterior se simplifica notablemente y conduce a la ecuación de la elipse en coordenadas rectangulares

${\displaystyle \begin{array}{l}{x}^2(1-{\epsilon }^2)+y^2=a^2(1-{\epsilon }^2)\\ \frac{x^2}{a^2}+\frac{y^2}{a^2(1-{\epsilon }^2)}=1\end{array}}$

Periodo

Se denomina periodo, al tiempo que tarda el móvil en dar una vuelta completa. En la figura vemos que el radio vector que une el Sol con el planeta barre en el intervalo de tiempo comprendido entre t y t+dt el área de color rojo de forma triangular.

El ángulo del vértice de dicho triángulo es dθ y la base del triángulo es un arco de longitud rdθ. El área del triángulo  es (base por altura dividido por dos)

${\displaystyle \frac{r(r·d\theta )}{2}=\frac{r^{2}d\theta }{2}}$

El vector posición de cualquier planeta respecto del Sol, barre áreas iguales de la elipse en tiempos iguales.

La ley de las áreas es equivalente a la constancia del momento angular, es decir, cuando el planeta está más alejado del Sol (afelio) su velocidad es menor que cuando está más cercano al Sol (perihelio). En el afelio y en el perihelio, el momento angular L es el producto de la masa del planeta, por su velocidad y por su distancia al centro del Sol.

Integrando la ecuación del momento angular expresado en coordenadas polares

${\displaystyle {\displaystyle \underset{0}{\overset{2\pi }{\int }}{r}^{2}d\theta }={\displaystyle \underset{0}{\overset{P}{\int }}\frac{L}{m}}dt{\displaystyle \underset{0}{\overset{2\pi }{\int }}\frac{r^{2}d\theta }{2}}=\frac{L}{2m}{\displaystyle \underset{0}{\overset{P}{\int }}dt}}$

La primera integral es el área total de la elipse πab, que es igual a la suma de las áreas de todos triángulos infinitesimales. La integral del segundo miembro es el periodo P del planeta, por tanto

${\displaystyle P=\frac{2m\pi ab}{L}=\frac{2\pi GM}{{\left(-2\frac{E}{m}\right)}^{3/2}}}$

Poniendo el semieje b en función del semieje a, (final del apartado anterior) llegamos a la fórmula que relaciona el periodo de la órbita de un planeta P y el semieje mayor de la elipse a, denominada tercera ley de Kepler.

${\displaystyle P^2=\frac{4{\pi }^2{a}^3}{\left(GM\right)}}$

Datos de los planetas

Escribimos un script para dibujar las órbitas elípticas de los planetas: Venus, Tierra, Marte y Júpiter conocidos sus ejes mayores en U.A. (Unidades Astronómicas) y sus excentricidades.

a=[0.723 1 1.524 5.203];
ex=[0.007 0.017 0.093 0.048];
ang=linspace(0,2*pi,100);
colores=['r' 'b' 'k' 'g'];
planeta=['V' 'T' 'M' 'J'];
hold on
for i=1:length(a)
    d=a(i)*(1-ex(i)^2);
    r=d./(1+ex(i)*cos(ang));
    x=r.*cos(ang);
    y=r.*sin(ang);
    plot(x,y,colores(i),'displayName',planeta(i)) 
end
axis equal
grid on
xlabel('x')
ylabel('y')
title('Sistema solar')
legend('-DynamicLegend','location','SouthEast')
hold off 

Relacionamos el semieje mayor de la trayectoria elíptica a con el periodo P. Representamos en el eje X, lna y en el eje Y, lnP. Los datos se ajustan a una línea recta, cuya pendiente determinamos.

a=[0.387,0.723,1,1.524,5.203,9.539,19.182,30.058]; %semieje mayor (UA)
P=[0.24,0.62,1,1.88,11.86,29.46,84.01,164.8]; %periodo (años)
plot(log(a),log(P), 'ro', 'markersize',4,'markeredgecolor' ,'r'
,'markerfacecolor','r')
xlabel('ln(a)')
ylabel('ln(P)')
grid on
title('Relación semieje mayor, periodo')

En el menú seleccionamos Tools/Basic Fitting, aparece un cuadro de diálogo donde marcamos la casilla linear debajo de Types of fit

La pendiente de la recta es 3/2, lnP=3·lna/2 o bien,

${\displaystyle P=a^{\frac{3}{2}},P^2=a^3}$

Los cuadrados de los periodos P (en años) son iguales a los cubos del semieje mayor a (en unidades UA) de la elipse (Tercera ley de Kepler)

Ejemplo 1

Se dispara un proyectil desde el polo Norte de la Luna con velocidad de v₀=2200 m/s. Calcular el tiempo de vuelo. Datos

• Masa de la Luna, M=7.349·10²² kg
• Radio de la Luna, R=1737. 4 km
• Constante, G=6.67·10^-11 N·m²/kg²

La ecuación de la elipse en coordenadas polares es

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

Para θ=π/2, r=R el radio de la Luna, por que el parámetro d=R.

El momento angular vale

${\displaystyle L=m\sqrt{GMR}}$

La energía en la posición de disparo

${\displaystyle E=\frac{1}{2}m{v}_0^2-G\frac{Mm}{R}}$

Conocida la energía E y el momento angular L, calculamos la excentricidad ε

${\displaystyle \epsilon =\sqrt{\frac{R{v}_0^2}{GM}-1}}$

El semieje mayor a de la elipse es

${\displaystyle \begin{array}{l}2a=\frac{R}{1-\epsilon }+\frac{R}{1+\epsilon }\\ a=\frac{R}{1-{\epsilon }^2}\end{array}}$

El semieje menor b

${\displaystyle \begin{array}{l}\epsilon =\frac{c}{a}\\ b=\sqrt{a^2-c^2}=a\sqrt{1-{\epsilon }^2}\end{array}}$

Para calcular el tiempo empleado por el proyectil desde su disparo hasta que llega al punto de impacto, tenemos que determinar el área barrida por el radio vector, desde la posición inicial, π/2 a la final, 3π/2, el área A sombreada en la figura

${\displaystyle \begin{array}{l}{\displaystyle \underset{\pi /2}{\overset{3\pi /2}{\int }}\frac{r^{2}d\theta }{2}}=\frac{L}{2m}{\displaystyle \underset{0}{\overset{T}{\int }}dt}\hfill \\ A=\frac{L}{2m}T\hfill \end{array}}$

La ecuación de la elipse en coordenadas rectangulares es

${\displaystyle \frac{x^2}{a^2}+\frac{y^2}{b^2}=1}$

La mitad A/2 del área sombreada es

${\displaystyle {\displaystyle \underset{-a}{\overset{c}{\int }}y·dx=}{\displaystyle \underset{-a}{\overset{c}{\int }}b\sqrt{1-\frac{x^2}{a^2}}}dx}$

Hacemos el cambio de variable, x=a·sinz, dx=a·cosz·dz

${\displaystyle \frac{A}{2}=b{\displaystyle \underset{-a}{\overset{c}{\int }}\sqrt{1-\frac{x^2}{a^2}}dx}=ab{\displaystyle \underset{z_1}{\overset{z_2}{\int }}{\cos }^{2}z·dz=}\frac{ab}{2}{\displaystyle \underset{z_1}{\overset{z_2}{\int }}\left(1+\cos (2z)\right)dz=\frac{ab}{2}\left(z_2-z_1+\frac{1}{2}\sin (2z_2)-\frac{1}{2}\sin (2z_1)\right)}}$

z₁=-π/2, z₂=arccos(c/a)=arccos(ε)

El tiempo de viuelo del proyectil entre el polo Norte y el polo Sur de la Luna es

${\displaystyle \begin{array}{l}A=ab\left(\arcsin (\epsilon )+\frac{\pi }{2}+\frac{1}{2}\sin (2\arcsin (\epsilon ))\right)=ab\left(\arcsin (\epsilon )+\frac{\pi }{2}+\epsilon \sqrt{1-{\epsilon }^2}\right)\\ T=\frac{2A}{\sqrt{GMR}}\end{array}}$

R=1737.4*1000; %radio de la Luna
M=7.349e22; %masa de la Luna
G=6.67e-11; %constante G
v0=2200; %velocidad de disparo
E=v0^2/2-G*M/R; %energía
L=sqrt(G*M*R); %momento angular
ex=sqrt(v0^2*R/(G*M)-1); %excentricidad
a=R/(1-ex^2); %semieje a
b=a*sqrt(1-ex^2); %semieje b
A=a*b*(asin(ex)+pi/2+ex*sqrt(1-ex^2)); %área
T=2*A/sqrt(G*M*R); %tiempo
disp(T)

El tiempo de vuelo del proyectil es T=41 309 s

4.1309e+04

Ejemplo 2

Este ejemplo se estudia la trayectoria de un proyectil disparado desde una altura h por encima de la Tierra en dirección perpendicular a la radial. El proyectil describe una órbita elíptica en uno de cuyos focos está el centro de la Tierra y en un plano perpendicular al ecuatorial chocando con la Tierra en un punto de latitud λ

En la página titulada Trayectoria de un proyectil disparado desde una altura h sobre la superficie de la Tierra se estudia este problema en profundidad

Para dibujar parte de esta figura se ha empleado el código

R=1;
lambda=pi/6;
alfa=sqrt((1-cos(lambda))/(2-cos(lambda)));
e=1-alfa^2;
d=2*alfa^2*R;

hold on
fplot(@(t) R*cos(t), @(t) R*sin(t),[0,2*pi])
ang=linspace(0,2*pi,180);
r=d./(1+e*cos(ang));
x=2*d*e/(1-e^2)+r.*cos(ang);
y=r.*sin(ang);
plot(x,y)
plot(R*cos(lambda),R*sin(lambda),'ro','markersize',3,'markeredgecolor','r',
'markerfacecolor','r')
plot(0,0,'ko','markersize',2,'markeredgecolor','k','markerfacecolor','k')
axis equal
axis off

El proyectil de masa m se dispara desde una altura h=R igual al radio de la Tierra, en dirección perpendicular al radio vector, con velocidad.

${\displaystyle v=\alpha \sqrt{\frac{GM}{2R}}}$

Determinaremos el valor del parámetro α para que el proyectil impacte en una localidad del hemisferio norte situada a una latitud λ.

Recuérdese que el término que multiplica al parámetro α es la velocidad de un cuerpo celeste que describe una órbita circular de radio 2R. Aplicando la dinámica del movimiento circular uniforme

${\displaystyle m\frac{v^2}{2R}=G\frac{Mm}{{\left(2R\right)}^2},v=\sqrt{\frac{GM}{2R}}}$

La energía y momento angular del proyectil es

${\displaystyle \begin{array}{l}E=\frac{1}{2}m{v}^2-\frac{GMm}{2R}=-\frac{1}{2}\frac{GMm}{R}\left(1-\frac{1}{2}{\alpha }^2\right)\\ L=mv(2R)=2mR\alpha \sqrt{\frac{GM}{2R}}\end{array}}$

La trayectoria es elíptica (E<0) si $\alpha <\sqrt{2}$

La ecuación de la elipse en coordenadas polares es

${\displaystyle \begin{array}{l}r=\frac{d}{1+\epsilon \cos \theta }\\ d=\frac{L^2}{GM{m}^2}=2{\alpha }^{2}R\\ \epsilon =\sqrt{1+\frac{2L^{2}E}{m^3{G}^2{M}^2}}=1-{\alpha }^2\end{array}}$

El semieje mayor a, menor b y semidistancia focal c son

${\displaystyle \begin{array}{l}a=\frac{d}{1-{\epsilon }^2}=\frac{2R}{2-{\alpha }^2}\\ c=\epsilon a=2R\frac{1-{\alpha }^2}{2-{\alpha }^2}\\ b^2=a^2(1-{\epsilon }^2)=\frac{4R^2{\alpha }^2}{2-{\alpha }^2}\end{array}}$

Impacto del proyectil

Calculamos el punto de intersección de la elipse uno de cuyos focos es el centro de la Tierra y la circunferencia de radio R

${\displaystyle \begin{array}{l}\{\begin{array}{l}{x}^2+y^2=R^2\\ \frac{{\left(x-c\right)}^2}{a^2}+\frac{y^2}{b^2}=1\end{array}\\ x=\frac{Ra-b^2}{c}=\frac{Ra-a^2(1-{\epsilon }^2)}{\epsilon a}=\frac{R-d}{\epsilon }=R\frac{1-2{\alpha }^2}{1-{\alpha }^2}\end{array}}$

Relacionamos el parámetro α con la latitud λ del punto de impacto

${\displaystyle \begin{array}{l}R\cos \lambda =R\frac{1-2{\alpha }^2}{1-{\alpha }^2}\\ \alpha =\sqrt{\frac{1-\cos \lambda }{2-\cos \lambda }}\end{array}}$

Supongamos que queremos que el proyectil impacte en una localidad del círculo polar ártico, λ=66.5°, el parámetro α=0.6128

Hay otro valor del parámetro α que da lugar a un impacto en el círculo polar ártico, la posición del punto de impacto es simétrica respecto del eje Y

${\displaystyle \begin{array}{l}-R\cos \lambda =R\frac{1-2{\alpha }^2}{1-{\alpha }^2}\\ \alpha =\sqrt{\frac{1+\cos \lambda }{2+\cos \lambda }}\end{array}}$

El valor mínimo de α se obtiene para λ=0, α=0, caída libre del cuerpo. El valor máximo de α se obtiene para λ=π (180°), $\alpha =\sqrt{2/3}$

Representamos tras trayectorias del proyectil para tres valores del parámero α=0.6128, 0.7636 y 0.8165

R=1; %radio
hold on
ang=linspace(0,2*pi,180);
fill(R*cos(ang),R*sin(ang),'c')
fplot(@(t) R*cos(t), @(t) R*sin(t),[0,2*pi], 'color','b')
for lambda=[66.5, -66.5, 180]*pi/180 %latitud
    alfa=sqrt((1-sign(lambda)*cos(lambda))/(2-sign(lambda)*cos(lambda))); 
    ex=1-alfa^2; %excentricidad
    d=2*alfa^2*R;
    r=d./(1+ex*cos(ang));
    x=2*d*ex/(1-ex^2)+r.*cos(ang);
    y=r.*sin(ang);
    plot(x,y)
    plot(sign(lambda)*R*cos(lambda),sign(lambda)*R*sin(lambda),'ro',
'markersize',3,'markeredgecolor','r','markerfacecolor','r')
end
plot(0,0,'ko','markersize',2,'markeredgecolor','k','markerfacecolor','k')
hold off
axis equal
axis off

Tiempo que tarda el proyectil en llegar al punto de impacto

Para calcular el tiempo empleado por el proyectil desde su disparo hasta que llega al punto de impacto, tenemos que calcular el área barrida por el radio vector, desde la posición inicial a la final, el área sombreada en gris en la figura

${\displaystyle \begin{array}{l}{\displaystyle \underset{0}{\overset{\lambda }{\int }}\frac{r^{2}d\theta }{2}}=\frac{L}{2m}{\displaystyle \underset{0}{\overset{T}{\int }}dt}\\ A=\frac{L}{2m}T\end{array}}$

Donde A se refiere al área barrida por el radio vector

El área sombreada es la suma del área de un triángulo de base Rcosλ y altura Rsinλ y el área de la parte de elipse comprendida entre -c+Rcosλ y a.

El primer sumando vale

$\frac{R\cos \lambda ·R\sin \lambda }{2}=R^2\frac{\sin (2\lambda )}{4}$

La ecuación de la elipse es

${\displaystyle \frac{x^2}{a^2}+\frac{y^2}{b^2}=1}$

El segundo sumando vale

${\displaystyle {\displaystyle \underset{-c+R\cos \lambda }{\overset{a}{\int }}b\sqrt{1-\frac{x^2}{a^2}}}dx}$

Hacemos el cambio de variable, x=a·sinz, dx=a·cosz·dz

${\displaystyle {\displaystyle \underset{-c+R\cos \lambda }{\overset{a}{\int }}b\sqrt{1-\frac{x^2}{a^2}}}dx=ab{\displaystyle \underset{z_1}{\overset{\pi /2}{\int }}{\cos }^{2}z·dz}=\frac{ab}{2}{\displaystyle \underset{z_1}{\overset{\pi /2}{\int }}\left(1+\cos (2z)\right)·dz}=\frac{ab}{2}\left(\frac{\pi }{2}-z_1-\frac{1}{2}\sin (2z_1)\right)}$

El área total es

${\displaystyle \begin{array}{l}A=\frac{R^2\sin (2\lambda )}{4}+\frac{ab}{2}\left(\frac{\pi }{2}-z_1-\frac{1}{2}\sin (2z_1)\right)\\ z_1=\arcsin \left(\frac{-c+R\cos \lambda }{a}\right)\end{array}}$

R=6.67e6; %radio de la Tierra
M=5.98e24; %masa de la Tierra
G=6.67e-11; %constante G
lambda=66.5*pi/180; %latitud
alfa=sqrt((1-cos(lambda))/(2-cos(lambda))); %parámetro
ex=1-alfa^2; %excentricidad
d=2*alfa^2*R;
a=2*R/(2-alfa^2); %semieje mayor
c=ex*a; %semidistancia focal
b=sqrt(a^2-c^2); %semieje menor

z1=asin((-c+R*cos(lambda))/a);
area=R^2*sin(2*lambda)/4+a*b*(pi/2-z1-sin(2*z1)/2)/2;
v=alfa*sqrt(G*M/(2*R));
T=2*area/(v*2*R);
disp(T/60)

El proyectil tarda 48.5 minutos en llegar al punto de impacto

   48.5392

Referencias

O. L. de Lange, J. Pierrus. Solved Problems in Classical Mechanics. Analytical and numerical solutions with comments. Oxford University Press (2010). Question 8.20, pp. 260-262

2020 Online Physics Olympiad: Open Contest. pp. 48-49. Problem 51