% TRAZADO BÁSICO DE LA ÓRBITA TERRESTRE (Curso de Sigüenza 2011)

% PARÁMETROS GENERALES (en unidades SI: m, kg, s)
G=6.673e-11; %constante de gravitación universal
M=1.989e30; %masa del Sol
UA=1.496e11; % 1 Unidad astronómica = 1.5e11m (distancia media Sol-Tierra)

%DATOS DE LA TIERRA (servidor de efemérides del OAN)
m=5.98e24; %masas
a=1*UA; %semieje mayor 20110721
ex=0.021; %excentricidad 20110721
T=365*86400; %período orbital aprox en segundos

%PLANTEAMIENTO DE LAS ECUACIONES DIFERENCIALES
function fdot = keplerecus (f,t)
      G=6.673e-11; M=1.989e30;  %hay que repetirlas porque las funciones son así de independientes 
      r3=(f(1)^2+f(2)^2)^(3/2); %para que el denominador luego quede más recogidito
      fdot(1) = f(3);
      fdot(2) = f(4);
      fdot(3) = -G*M*f(1)/(r3);
      fdot(4) = -G*M*f(2)/(r3);
endfunction

%% TIERRA
 Xini=1*a; Yini=0;
 VXini=0; VYini=1*(2*pi*a/(365*24*60*60));;
 init=[Xini Yini VXini VYini]
 t = linspace (0, 1*T, 800); %he pedido dividir una órbita entera (tiempo T) en 800 posiciones
 q = lsode ("keplerecus", init, t); %las 800 posiciones de la Tierra se guardarán en el vector q
%Dibujo el Sol en el centro y la órbita de la Tierra (800 posiciones) en azul. Pongo etiquetas. 
 plot(0,0,'k*',q(:,1),q(:,2),'b'); xlabel("X"); ylabel("Y"); title("trayectoria"); axis("equal");

 
%COMPROBACIONES 

%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=abs(ae/a)

%velocidades
 %figure(2); plot(t/86400,q(:,3)/1000,t/86400,q(:,4)/1000); xlabel("tiempo"); ylabel("velocidades"); %evolución de las velocidades

%probar diferentes tmax (más o menos de una órbita), Xini, VYini 

return

