%problema de Kepler - masa M dominante, origen en M - potencial postnewtoniano

global GM Rs
% Parámetros (CGS)
G=6.673e-8;
msol=1.989e33;
%GM = G*6.75*msol; % GM (cgs)
GM = G*msol; % GM (cgs)
c=3e10;
Rs=2*GM/c^2

%Rini=5*Rs; Vini=0.38*c;
%Rini=5*Rs; Vini=sqrt(0.4*GM/Rs);
%Rini=3.2*Rs; Vini=1*sqrt(GM*Rini/(Rini-Rs)^2);
%Rini=2.9*Rs; Vini=1*sqrt(GM*Rini/(Rini-Rs)^2);
%probar con Vini=1.001*sqrt y Vini=0.999*sqrt en los dos casos anteriores

%MERCURIO
UA=1.496e13; % 1 Unidad astronómica = 1.5e13cm (distancia media Sol-Tierra)
a=0.386*UA; %distancia media Mercurio/Sol
ex=0.203;
m=3.302e26;
vp=sqrt(G*(msol+m).*(1+ex)./a./(1-ex)) %velocidad en el perihelio
dp=a.*(1-ex); %distancia al perihelio
Rini=dp; Vini=vp;

disp("unidades cgs")
init=[Rini; 0; 0; Vini]

function fdot = postnewtpot (f,t)
      global GM Rs
      r=(f(1)^2+f(2)^2)^(1/2);
      den=r*(r-Rs)^2;
      fdot(1) = f(3);
      fdot(2) = f(4);
      fdot(3) = -GM*f(1)/den; %potencial postnewtoniano
      fdot(4) = -GM*f(2)/den;
endfunction
%s = lsode ("postnewtpot", init, (t = linspace (0, .03, 30000)')); %Rini=5*Rs; Vini=0.38*c; %variar 1Msol hasta 6.75Msol
%s = lsode ("postnewtpot", init, (t = linspace (0, .05, 1000)')); %Rini=5*Rs; Vini=sqrt(0.4*GM/Rs); %M=6.75msol
%s = lsode ("postnewtpot", init, (t = linspace (0, .003, 30000)')); %estable: Rini=3.2*Rs; Vini=1*sqrt(GM*Rini/(Rini-Rs)^2);
%s = lsode ("postnewtpot", init, (t = linspace (0, .0019, 30000)')); %inestable (PCcrash): Rini=2.9*Rs; Vini=1*sqrt(GM*Rini/(Rini-Rs)^2);
%s = lsode ("postnewtpot", init, (t = linspace (0, .001882, 30000)')); %inestable (no crash): Rini=2.9*Rs; Vini=1*sqrt(GM*Rini/(Rini-Rs)^2);
s = lsode ("postnewtpot", init, (t = linspace (0, 1000*88*86400, 10000)')); %caso de Mercurio


%s = lsode ("postnewtpot", init, (t = linspace (0, 0.0035, 30000)'));
%s = lsode ("postnewtpot", init, (t = linspace (0, 0.00315335, 30000)'));
%s = lsode ("postnewtpot", init, (t = linspace (0, 0.1, 30000)'));
%s = lsode ("postnewtpot", init, (t = linspace (0, 1, 30000)'));

figure(1);
plot(s(:,1),s(:,2)); xlabel("X"); ylabel("Y"); title("postnewtoniano")
axis("equal")
figure(2); 
axis()
plot(t,s(:,4)); xlabel("t"); ylabel("Vy");