clf
global gamk r_maa;
g_const=6.67408e-11;
m_maa=5.97219e+24;
r_maa=6371.;
gamk=g_const*m_maa*1e-9;
r0=[0.;6400.;0;0.];
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',1000.,'InitialStep',1.,'Events',@mycheck);
[t,res]=ode45(@func,[0,400000], r0,options);
hold on;
axis equal;
%axis([-2 2. -2. 2.]);
grid on;
z=res(:,1)*0;
sphere
plot3(res(:,1)./mkr,res(:,2)./mkr,z);
comet3(res(:,1)./mkr,res(:,2)./mkr,z,0.01);

function  dy=func(t,r)
global gamk;
dy=zeros(4,1);
dy(1)=r(3);
dy(2)=r(4);
dy(3)=-r(1)*gamk/sqrt((r(1)^2+r(2)^2))^3;
dy(4)=-r(2)*gamk/sqrt((r(1)^2+r(2)^2))^3;
end

function [h,finish,dir]= mycheck(t,y)
global r_maa;
rr=sqrt(y(1)^2+y(2)^2);
h=[rr-r_maa];
finish=[1]; 
dir=[-1];
end