%Función que propaga la posición de un planeta(TIERRA o JUPITER) a un instante t desde la
%época dados los elementos orbitales en la época
function EOt=pplanetaLPCP(t,EO)
EOt=zeros(1,6);
% T=2*pi*sqrt(EO(2)^3/(mus));
% n=2*pi/T;
if t==0
T0=0;
else
DeltaT=t/86400;
T0=DeltaT/36525;
end
EO(1:6)=EO(1:6)+EO(7:12)*T0;
M=EO(6)-EO(5);
M=2*pi*(M/2/pi-floor(M/2/pi));
EO(5)=EO(5)-EO(4);
EO(5)=2*pi*(EO(5)/2/pi-floor(EO(5)/2/pi));
myfun = @(E,M) M-E+EO(1)*sin(E);
OPTIONS=optimset('TolFun',1e-10,'TolX',1e-10,'MaxIter',5000,'MaxFunEvals',5000);
X = fzero(@(E) myfun(E,M),pi,OPTIONS);
TAepoch=2*atan(sqrt((1+EO(1))/(1-EO(1)))*tan(X/2));
EO(6)=2*pi*(TAepoch/2/pi-floor(TAepoch/2/pi));
EOt=EO(1:6);
% EOt(1:12)=EO(1:12);
% M=EOt(6)-EOt(5);
% EOt(5)=EOt(5)-EOt(4);
%
% myfun = @(E,M) M-E+EO(1)*sin(E);
% X = fzero(@(E) myfun(E,M),0.1);
% TAepoch=2*atan(sqrt((1+EO(1))/(1-EO(1)))*tan(X/2));
%
% EOt(6)=TAepoch;
% if floor(DeltaT/T)==0 %Hay q porpagar un tiempo menor al del periodo del planeta
% E=atan(tan(EO(6)/2)*sqrt((1-EO(1))/(1+EO(1))))*2;
% M=E-EO(1)*sin(E);
% TPOS=M/n;
% if (TPOS+DeltaT)>T
% incT=(TPOS+DeltaT)-T;
% M=n*incT;
% myfun = @(E,M) M-E+EO(1)*sin(E);
% X = fzero(@(E) myfun(E,M),0.1);
% TAepoch=2*atan(sqrt((1+EO(1))/(1-EO(1)))*tan(X/2));
% else
% incT=(TPOS+DeltaT);
% M=n*incT;
% myfun = @(E,M) M-E+EO(1)*sin(E);
% X = fzero(@(E) myfun(E,M),0.1);
% TAepoch=2*atan(sqrt((1+EO(1))/(1-EO(1)))*tan(X/2));
% end
% else %Hay que propagar mas de una órbita
% fractT=DeltaT/T-floor(DeltaT/T);
% DeltaT=fractT*T;
%
% E=atan(tan(EO(6)/2)*sqrt((1-EO(1))/(1+EO(1))))*2;
% M=E-EO(1)*sin(E);
% TPOS=M/n;
%
% if (TPOS+DeltaT)>T
% incT=(TPOS+DeltaT)-T;
% M=n*incT;
% myfun = @(E,M) M-E+EO(1)*sin(E);
% X = fzero(@(E) myfun(E,M),0.1);
% TAepoch=2*atan(sqrt((1+EO(1))/(1-EO(1)))*tan(X/2));
% else
% incT=(TPOS+DeltaT);
% M=n*incT;
% myfun = @(E,M) M-E+EO(1)*sin(E);
% X = fzero(@(E) myfun(E,M),0.1);
% TAepoch=2*atan(sqrt((1+EO(1))/(1-EO(1)))*tan(X/2));
% end
% end
end
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
