global gm gk omegak kor fi0;

%parameetri
gm=3.98604418e+5;
gk=4.9048695e+3;
mkr=6371.;
kkr=1738.
kor=384399.;
kok=1.022;
omegak=kok/kor;
%fi0=-pi*52.7145/180.;
fi0=-pi*60./180.;
%fi0=-pi/3.62;

%alg tingimused
r0=[-6471.;0.;0.;-11.01];

options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',10000.,'InitialStep',1.);
[t,r]=ode45(@kplf2,[0,400000], r0,options);

grid on;

%z=r(:,1)*0;
%sphere

x=kor*cos(omegak*t+fi0)./mkr;
y=kor*sin(omegak*t+fi0)./mkr;

%comet3(r(:,1)./mkr,r(:,2)./mkr,z);
%for i=1:10:n
%quiver3(r(i,1)./mkr,r(i,2)./mkr,z(i),r(i,3),r(i,4),z(i),0.5,'g')
%end
%plot(x,y);
%comet3(r(:,1)./mkr,r(:,2)./mkr,z);
%plot(x(i),y(i));
%drawnow;
%hold on
%plot(t,vs)
%return
axis equal
hold on
nt=size(t)
h=animatedline();
h1=animatedline();
for i=1:nt
addpoints(h,r(i,1)/mkr,r(i,2)/mkr);addpoints(h1,x(i),y(i));
F(i)=getframe;
end
%movie(F,1)
 
return
%hold off;
%comet(r(:,1)./mkr,r(:,2)./mkr);
function  dy=kplf2(t,y)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
global gm gk omegak kor fi0;
dy=zeros(4,1);
xk=kor*cos(omegak*t+fi0);
yk=kor*sin(omegak*t+fi0);
xx=(y(1)-xk);
yy=(y(2)-yk);
rr=sqrt(xx^2+yy^2);
dy(1)=y(3);
dy(2)=y(4);
dy(3)=-y(1)*gm/sqrt((y(1)^2+y(2)^2))^3-xx*gk/rr^3;
dy(4)=-y(2)*gm/sqrt((y(1)^2+y(2)^2))^3-yy*gk/rr^3;
end