global gm gk omegak kor fi0 kkr mkr;

gm=3.98604418e+5;
gk=4.9048695e+3;
mkr=6378.;
kkr=1738.;
kor=384399.;
kok=1.022;
omegak=kok/kor;
%fi0=-pi*52.7145/180.;
%fi0=-pi*60./180.;
%fi0=-pi/3.62;
fi0=-pi*49.145/180.;
%alg tingimused
r0=[-6478.;0.;0.;-11.01];

options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',100.,'InitialStep',1.,'Events',@mycheck);
[t,r]=ode45(@func,[0,400000], r0,options);
hold on
x=kor*cos(omegak*t+fi0);
y=kor*sin(omegak*t+fi0);
v=sqrt(r(:,3).^2+r(:,4).^2);
h=sqrt((r(:,1)-x).^2+(r(:,2)-y).^2)-kkr;
min(h)

subplot(1,3,1)
plot(t,v)
grid on;

subplot(1,3,2)
plot(t,h)
grid on;

subplot(1,3,3)
grid on
axis equal
nt=size(t);
p1=animatedline('Color','blue');
p2=animatedline();
for i=1:100:nt
addpoints(p1,r(i,1)/mkr,r(i,2)/mkr);
addpoints(p2,x(i)/mkr,y(i)/mkr);
drawnow
end

function  dy=func(t,y)
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

function  [h,finish,dir]=mycheck(t,y)
global kkr kor omegak fi0
xk=kor*cos(omegak*t+fi0);
yk=kor*sin(omegak*t+fi0);
h=sqrt((y(1)-xk)^2+(y(2)-yk)^2)-kkr;
finish=1;
dir=0;
end



