% JBG
% problema de Kepler - masa M dominante (Sol), origen de coordenadas en Sol
% TODOS LOS PERIHELIOS ALINEADOS (en y=0) Y TODOS LOS PLANETAS SOBRE LA ECLÍPTICA

global GM  %constante que quiero que comparta la función keplerecus

% PARÁMETROS GENERALES (en unidades SI)
G=6.673e-11; %constante de gravitación universal
M=1.989e30; %masa del Sol
GM = G*M; % GM
UA=1.496e11; % 1 Unidad astronómica = 1.5e11m (distancia media Sol-Tierra)
JD= 2455764.375770; %Julian date 20110721 23h00 hora peninsular

%DATOS DE LOS PLANETAS VISIBLES A SIMPLE VISTA: Mercurio, Venus, Tierra, Marte, Júpiter, Saturno
m=[3.302e23 4.87e24 5.98e24 6.42e23 1.899e27 5.688e26]; %masas
a=[0.386 0.723 1 1.522 5.201 9.516]*UA; %semiejes mayores 20110721
ex=[0.203 0.018 0.021 0.089 0.049 0.055]; %excentricidades 20110721
inc=[7.004 3.395 0 1.849 1.304 2.488]; %inclinación de la órbita (grados) respecto eclíptica
w=[22.37 98.08 -94.64 -73.09 -85.18 -23.69]; %argumento del perihelio (grados) 20110721
OM=[48.32 76.65 175.38 49.52 100.51 113.61]; %longitud nodo ascendente (grados) 20110721
%vm=[47.87e3 35.02e3 30.29e3 24.13e3 13.07e3 9672.4] %velocidades medias
vp=sqrt(G*(M+m).*(1+ex)./a./(1-ex)); %velocidad en el perihelio
dp=a.*(1-ex); %distancia al perihelio
T=[0.24 0.62 1 1.88 11.86 29.46]*365.25*86400; %período orbital aprox en segundos

%ELECCIÓN DE PLANETAS A VISUALIZAR
ip=input("introduce planetas a visualizar (ej: [1:3 5] para Me, V, T, J): ");

%PLANTEAMIENTO DE LAS ECUACIONES DIFERENCIALES
function fdot = keplerecus (f,t)
      global GM
      r3=(f(1)^2+f(2)^2)^(3/2);
      fdot(1) = f(3);
      fdot(2) = f(4);
      fdot(3) = -GM*f(1)/(r3);
      fdot(4) = -GM*f(2)/(r3);
endfunction

%Órbitas elípticas: Vini<sqrt(2*GM/Rini) según 3a ley de Kepler
disp("caso elíptico, origen de coordenadas en foco (masa dominante)")
plot(0,0,'k*'); xlabel("X"); ylabel("Y"); title("trayectorias"); axis("equal");

%% TIERRA
if(!all(ip-3)) %si está dentro del menú de planetas elegido
 hold on
 disp(""); disp("Tierra:")
 Xini=dp(3); Yini=0;
 VXini=0; VYini=0.3*vp(3);
 %VYini=2*pi*a(3)/(365*24*60*60); % velocidad media aproximada de la Tierra
 %VYini=1*sqrt(2*GM/Xini); %probar 1, 0.99, 0.71~=sqrt(0.5) en que vp~=va; por debajo de aquí la velocidad inicial es tan baja que no es la del perihelio sino la del afelio: 0.5, 0.3 
 init=[Xini Yini VXini VYini]
 t = linspace (0, 10*T(3), 1800);
 q = lsode ("keplerecus", init, t);
 plot(0,0,'k*',q(:,1),q(:,2),'b');
 %estimación de a y ecc (válida sólo si se completa al menos una órbita)
 a=(max(q(:,1))-min(q(:,1)))/2 %(apastro+periastro)/2
 ae=(q(find(q(:,2)==max(q(:,2))),1)+q(find(q(:,2)==min(q(:,2))),1))/2; %abscisa (promediada) de los puntos de max y min ordenada
 ecc=ae/a
end %if Tierra

%% MERCURIO
if(!all(ip-1))
 hold on
 disp(""); disp("Mercurio:")
 Xini=dp(1); Yini=0;
 VXini=0; VYini=vp(1);
 %VYini=2*pi*a(1)/(88*24*60*60); % velocidad media aproximada de Mercurio
 %VYini=0.3*sqrt(2*GM/Xini); %probar 1, 0.99, 0.71~=sqrt(0.5) en que vp~=va; por debajo de aquí la velocidad inicial es tan baja que no es la del perihelio sino la del afelio: 0.5, 0.3 
 init=[Xini Yini VXini VYini]
 t = linspace (0, T(1), 800);
 q = lsode ("keplerecus", init, t);
 plot(q(:,1),q(:,2),'r');
 a=(max(q(:,1))-min(q(:,1)))/2 %(apastro+periastro)/2
 ae=(q(find(q(:,2)==max(q(:,2))),1)+q(find(q(:,2)==min(q(:,2))),1))/2; %abscisa (promediada) de los puntos de max y min ordenada
 ecc=ae/a
end %if Mercurio

%% VENUS
if(!all(ip-2))
 hold on
 disp(""); disp("Venus:")
 Xini=dp(2); Yini=0;
 VXini=0; VYini=vp(2);
 init=[Xini Yini VXini VYini]
 t = linspace (0, T(2), 800);
 q = lsode ("keplerecus", init, t);
 plot(q(:,1),q(:,2),'g');
 a=(max(q(:,1))-min(q(:,1)))/2 %(apastro+periastro)/2
 ae=(q(find(q(:,2)==max(q(:,2))),1)+q(find(q(:,2)==min(q(:,2))),1))/2; %abscisa (promediada) de los puntos de max y min ordenada
 ecc=ae/a
end %if Venus

%% MARTE
if(!all(ip-4))
 hold on
 disp(""); disp("Marte:")
 Xini=dp(4); Yini=0;
 VXini=0; VYini=vp(4);
 init=[Xini Yini VXini VYini]
 t = linspace (0, T(4), 800);
 q = lsode ("keplerecus", init, t);
 plot(q(:,1),q(:,2),'c');
 a=(max(q(:,1))-min(q(:,1)))/2 %(apastro+periastro)/2
 ae=(q(find(q(:,2)==max(q(:,2))),1)+q(find(q(:,2)==min(q(:,2))),1))/2; %abscisa (promediada) de los puntos de max y min ordenada
 ecc=ae/a
end %if Marte

%% JUPITER
if(!all(ip-5))
 hold on
 disp(""); disp("Jupiter:")
 Xini=dp(5); Yini=0;
 VXini=0; VYini=vp(5);
 init=[Xini Yini VXini VYini]
 t = linspace (0, T(5), 800);
 q = lsode ("keplerecus", init, t);
 plot(q(:,1),q(:,2),'m');
 a=(max(q(:,1))-min(q(:,1)))/2 %(apastro+periastro)/2
 ae=(q(find(q(:,2)==max(q(:,2))),1)+q(find(q(:,2)==min(q(:,2))),1))/2; %abscisa (promediada) de los puntos de max y min ordenada
 ecc=ae/a
end %if Jupiter

%% SATURNO
if(!all(ip-6))
 hold on
 disp(""); disp("Saturno:")
 Xini=dp(6); Yini=0;
 VXini=0; VYini=vp(6);
 init=[Xini Yini VXini VYini]
 t = linspace (0, T(6), 800);
 q = lsode ("keplerecus", init, t);
 plot(q(:,1),q(:,2),'k');
 a=(max(q(:,1))-min(q(:,1)))/2 %(apastro+periastro)/2
 ae=(q(find(q(:,2)==max(q(:,2))),1)+q(find(q(:,2)==min(q(:,2))),1))/2; %abscisa (promediada) de los puntos de max y min ordenada
 ecc=ae/a
end %if Saturno

hold off
return

