% JBG
% problema de Kepler - masa M dominante (Sol), origen de coordenadas en Sol
% PERIHELIOS REALES Y TODOS LOS PLANETAS EN EL PLANO ECLÍPTICO

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

% Parámetros (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

%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
%ap=[77.77 131.61 102.97 336.10 14.59 -999]; %argumentos del perihelio??
thrad=(w+OM)*pi/180; %ángulo perihelio desde Aries (caso de planetas en plano eclíptico)
de=27+365+31+28+31+30+31+30+21; % días extra a partir de 20091204 ... hasta el 20110721
%de=27+31+28+31+3; %hasta 20100403
%de=0; %hasta 20091204
%de=27+365+31+28+16; % días extra a partir de 20091204 ... hasta 20110316 (perihelio Júpiter)
dd=[61 62 334 227 3866 2384]; %días desde perihelio hasta 20091204
%dd=[61 -164 338 220 3856 2384]; %días desde perihelio hasta 20091204
tp=(dd+de)*86400; %tiempo desde perihelio


%MODO DE USO
disp("importante mejorar la inicialización: fechas de perihelios !!")
ip=input("introduce planetas a visualizar (ej: [1:3 5] para Me, V, T, J): ");

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

plot(0,0,'k*'); xlabel("X"); ylabel("Y"); title("trayectorias"); axis("equal");

%Tiempo de integración
%tmax=T(max(ip))*ones(1,6);
tmax=tp; tmax/86400

%% TIERRA
if(!all(ip-3))
 hold on
 disp(""); disp("Tierra:")
 Xini=dp(3)*cos(thrad(3)); Yini=dp(3)*sin(thrad(3));
 VXini=-vp(3)*sin(thrad(3)); VYini=vp(3)*cos(thrad(3));
 %VYini=2*pi*UA/(365*24*60*60); % velocidad media aproximada de la Tierra
 %VYini=1*sqrt(2*GM/Rini);
 init=[Xini Yini VXini VYini]
 t = linspace (0, tmax(3), 459);
 q = lsode ("keplerecus", init, t);
 xx=linspace(0,Xini,10); %para dibujar el vector perihelio
 plot(xx,Yini*xx/Xini,'k')
 plot(0,0,'k*',q(:,1),q(:,2),'b',q(459,1),q(459,2),'b*');
 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)*cos(thrad(1)); Yini=dp(1)*sin(thrad(1));
 VXini=-vp(1)*sin(thrad(1)); VYini=vp(1)*cos(thrad(1));
 %VYini=2*pi*UA/(365*24*60*60); % velocidad media aproximada de la Tierra
 %VYini=1*sqrt(2*GM/Rini);
 init=[Xini Yini VXini VYini]
 t = linspace (0, tmax(1), 185);
 q = lsode ("keplerecus", init, t); 
 xx=linspace(0,Xini,10);
 plot(xx,Yini*xx/Xini,'k')
 plot(q(:,1),q(:,2),'r',q(185,1),q(185,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)*cos(thrad(2)); Yini=dp(2)*sin(thrad(2));
 VXini=-vp(2)*sin(thrad(2)); VYini=vp(2)*cos(thrad(2));
 %VYini=2*pi*UA/(365*24*60*60); % velocidad media aproximada de la Tierra
 %VYini=1*sqrt(2*GM/Rini);
 init=[Xini Yini VXini VYini]
 t = linspace (0, tmax(2), 800);
 q = lsode ("keplerecus", init, t);
 %xx=linspace(0,Xini,10); %para dibujar el vector perihelio
 %plot(xx,Yini*xx/Xini,'k'); 
 plot(q(:,1),q(:,2),'g',q(800,1),q(800,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)*cos(thrad(4)); Yini=dp(4)*sin(thrad(4));
 VXini=-vp(4)*sin(thrad(4)); VYini=vp(4)*cos(thrad(4));
 %VYini=2*pi*UA/(365*24*60*60); % velocidad media aproximada de la Tierra
 %VYini=1*sqrt(2*GM/Rini);
 init=[Xini Yini VXini VYini]
 t = linspace (0, tmax(4), 800);
 q = lsode ("keplerecus", init, t);
 %xx=linspace(0,Xini,10);
 %plot(xx,Yini*xx/Xini,'k')
 plot(q(:,1),q(:,2),'c',q(800,1),q(800,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)*cos(thrad(5)); Yini=dp(5)*sin(thrad(5));
 VXini=-vp(5)*sin(thrad(5)); VYini=vp(5)*cos(thrad(5));
 %VYini=2*pi*UA/(365*24*60*60); % velocidad media aproximada de la Tierra
 %VYini=1*sqrt(2*GM/Rini);
 init=[Xini Yini VXini VYini]
 t = linspace (0, tmax(5), 800);
 q = lsode ("keplerecus", init, t);
 xx=linspace(0,Xini,10);
 plot(xx,Yini*xx/Xini,'k')
 plot(q(:,1),q(:,2),'m',q(800,1),q(800,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)*cos(thrad(6)); Yini=dp(6)*sin(thrad(6));
 VXini=-vp(6)*sin(thrad(6)); VYini=vp(6)*cos(thrad(6));
 %VYini=2*pi*UA/(365*24*60*60); % velocidad media aproximada de la Tierra
 %VYini=1*sqrt(2*GM/Rini);
 init=[Xini Yini VXini VYini]
 t = linspace (0, tmax(6), 800);
 q = lsode ("keplerecus", init, t);
 xx=linspace(0,Xini,10);
 plot(xx,Yini*xx/Xini,'k')
 plot(q(:,1),q(:,2),'k',q(800,1),q(800,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

%t=input("días transcurridos desde perihelio: ");
%t=27+28+31+30+31+30+31+31+30+23 %Tierra de 20090104 a 20091023
%t= 19   %Mercurio de 20091004 a 20091023
%t= 19.7 %Venus 20100516-224.701 (=20091003) a 20091023
%t= 9+31+30+31+31+30+23   %Marte de 20090421 a 20091023
%t= 26+30+31+31+30+23+366*3+365*7   %Júpiter 19990505 a 20091023
 
%http://almanako.com/aphelion_perihelion.shtml
%20091004 Mercury at perihelion 0.30750 AU (28,590,278 miles or 46,001,757 km)
%20100516 Venus at perihelion 0.71840 AU (66,794,170 miles or 107,471,819 km)
%20090421 Mars at perihelion
%19990505, 20110316 Jupiter at perihelion
%http://www.astronomycafe.net/qadir/q47.html
%20030526 Saturn at perihelion
%20500301 Uranus at perihelion
%Unless I have made some dumb mistake in adding and converting from elapsed days to a calendar date, these dates should be more or less correct. Jupiter's last Aphelion was on June 12, 1993, and Saturn's was on September 16, 1988. Uranus was at perihelion on June 4, 1966. 

