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
vs=sqrt(r(:,3).^2+r(:,4).^2);
%rs=(sqrt((r(:,1)./mkr-x(:)).^2+(r(:,2)./mkr-y(:)).^2)-kkr/mkr)
plot(t,vs)
return
axis equal
hold on
nt=size(t)
for i=1:nt
plot(r(i,1)/mkr,r(i,2)/mkr);plot(x(i),y(i));
F(i)=getframe;
end
movie(F,1)
pause 
return
%hold off;
%comet(r(:,1)./mkr,r(:,2)./mkr);
