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

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; %período orbital aprox en días
thrad=(w+OM)*pi/180; %ángulo perihelio desde Aries (caso de planetas en plano eclíptico)

%Fecha de inicio (20091204 por defecto). Luego introduciré un contador !!!!!!!!!!!!!
ti=[61 62 334 227 3866 2384]; %días desde perihelio hasta 20091204

%Fecha de fin (20110721 por defecto). Luego introduciré un contador !!!!!!!!!!!!!
te=27+365+31+28+31+30+31+30+21; % días extra a partir de 20091204 ... hasta el 20110721
%te=27+31+28+31+3; % días extra a partir de 20091204 ... hasta el 20100403
tf=ti+te;

%MODO DE USO
ip=input("introduce planeta para ver trayectoria geocentrica: ");

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)")
axis("equal");

%% TIERRA (siempre hay que calcularla)
 Xini=dp(3)*cos(thrad(3)); Yini=dp(3)*sin(thrad(3));
 VXini=-vp(3)*sin(thrad(3)); VYini=vp(3)*cos(thrad(3));
 init=[Xini Yini VXini VYini]
 t = 0:86400:tf(3)*86400; %tf+1 valores, de los cuales los te últimos representan el intervalo entre ti y tf
 qt = lsode ("keplerecus", init, t);
 kt=ti(3)+1:tf(3)+1;
 figure(1); plot(0,0,'k*',qt(:,1),qt(:,2),'b',qt(length(qt),1),qt(length(qt),2),'b*'); hold on

%% MERCURIO
if(ip==1)
 Xini=dp(1)*cos(thrad(1)); Yini=dp(1)*sin(thrad(1));
 VXini=-vp(1)*sin(thrad(1)); VYini=vp(1)*cos(thrad(1));
 init=[Xini Yini VXini VYini]
 t = 0:86400:tf(1)*86400;
 qp = lsode ("keplerecus", init, t);
 kp=ti(1)+1:tf(1)+1;
 plot(qp(:,1),qp(:,2),'r',qp(length(qp),1),qp(length(qp),2),'r*');hold off
 figure(2); plot(qp(kp,1)-qt(kt,1),qp(kp,2)-qt(kt,2)); xlabel("Xp-Xt metros"); ylabel("Yp-Yt metros"); title ("trayectoria geocentrica");
end %if Mercurio

%% VENUS
if(ip==2)
 Xini=dp(2)*cos(thrad(2)); Yini=dp(2)*sin(thrad(2));
 VXini=-vp(2)*sin(thrad(2)); VYini=vp(2)*cos(thrad(2));
 init=[Xini Yini VXini VYini]
 t = 0:86400:tf(2)*86400;
 qp = lsode ("keplerecus", init, t);
 kp=ti(2)+1:tf(2)+1;
 plot(qp(:,1),qp(:,2),'r',qp(length(qp),1),qp(length(qp),2),'r*');hold off
 figure(2); plot(qp(kp,1)-qt(kt,1),qp(kp,2)-qt(kt,2)); xlabel("Xp-Xt metros"); ylabel("Yp-Yt metros"); title ("trayectoria geocentrica");
end %if Venus

%% MARTE
if(ip==4)
 Xini=dp(4)*cos(thrad(4)); Yini=dp(4)*sin(thrad(4));
 VXini=-vp(4)*sin(thrad(4)); VYini=vp(4)*cos(thrad(4));
 init=[Xini Yini VXini VYini]
 t = 0:86400:tf(4)*86400;
 qp = lsode ("keplerecus", init, t);
 kp=ti(4)+1:tf(4)+1;
 plot(qp(:,1),qp(:,2),'r',qp(length(qp),1),qp(length(qp),2),'r*');hold off
 figure(2); plot(qp(kp,1)-qt(kt,1),qp(kp,2)-qt(kt,2)); xlabel("Xp-Xt metros"); ylabel("Yp-Yt metros"); title ("trayectoria geocentrica");
end %if Marte

%% JUPITER
if(ip==5)
 Xini=dp(5)*cos(thrad(5)); Yini=dp(5)*sin(thrad(5));
 VXini=-vp(5)*sin(thrad(5)); VYini=vp(5)*cos(thrad(5));
 init=[Xini Yini VXini VYini]
 t = 0:86400:tf(5)*86400;
 qp = lsode ("keplerecus", init, t);
 kp=ti(5)+1:tf(5)+1;
 plot(qp(:,1),qp(:,2),'r',qp(length(qp),1),qp(length(qp),2),'r*');hold off
 figure(2); plot(qp(kp,1)-qt(kt,1),qp(kp,2)-qt(kt,2)); xlabel("Xp-Xt metros"); ylabel("Yp-Yt metros"); title ("trayectoria geocentrica");
end %if Jupiter

%% SATURNO
if(ip==6)
 Xini=dp(6)*cos(thrad(6)); Yini=dp(6)*sin(thrad(6));
 VXini=-vp(6)*sin(thrad(6)); VYini=vp(6)*cos(thrad(6));
 init=[Xini Yini VXini VYini]
 t = 0:86400:tf(6)*86400;
 qp = lsode ("keplerecus", init, t);
 kp=ti(6)+1:tf(6)+1;
 plot(qp(:,1),qp(:,2),'r',qp(length(qp),1),qp(length(qp),2),'r*');hold off
 figure(2); plot(qp(kp,1)-qt(kt,1),qp(kp,2)-qt(kt,2)); xlabel("Xp-Xt metros"); ylabel("Yp-Yt metros"); title ("trayectoria geocentrica");
end %if Saturno

return