function [yp]=ecdifmov(t,y)
global mus mut EOTi muj EOJu
r=y(1:3)';%POSICIÓN DEL ASTEROIDE
%calcular el vector de posición del planeta en el instante t:TIERRA
%Usar el modelo de los dos cuerpos
% %
% EOt=pplanetaCP(t,EOTi);
% [coord]=rv(EOt);
% rpT=coord(1:3)';
%Usar el propagador planetario
EOt=pplanetaLPCP(t,EOTi);
[coord]=rv(EOt);
rpT=coord(1:3)';
%calcular el vector de posición del planeta en el instante t:JUPITER
%Usar el modelo de los dos cuerpos
% %
% EOt=pplanetaCP(t,EOJu);
% [coord]=rv(EOt);
% rpJ=coord(1:3)';
%Usar el propagador planetario
EOt=pplanetaLPCP(t,EOJu);
[coord]=rv(EOt);
rpJ=coord(1:3)';
% if norm(rpT-r)>2e8
% pert=(mut*(rpT-r)/norm(rpT-r)^3-mut*rpT/norm(rpT)^3)+(muj*(rpJ-r)/norm(rpJ-r)^3-muj*rpJ/norm(rpJ)^3);%PERTURBACIÓN DEBIDA A LOS PLANETAS
% else
% pert=(muj*(rpJ-r)/norm(rpJ-r)^3-muj*rpJ/norm(rpJ)^3);
% end
% pert=[0 0 0];%PERTURBACIÓN DEBIDA A LOS PLANETAS
pert=(mut*(rpT-r)/norm(rpT-r)^3-mut*rpT/norm(rpT)^3)+(muj*(rpJ-r)/norm(rpJ-r)^3-muj*rpJ/norm(rpJ)^3);%PERTURBACIÓN DEBIDA A LOS PLANETAS
r=sqrt(y(1)^2+y(2)^2+y(3)^2);
yp=[y(4),y(5),y(6),-mus*y(1)/r^3+pert(1),-mus*y(2)/r^3+pert(2),-mus*y(3)/r^3+pert(3)]';
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
