global r_maa ro_atm_0 p1 p2 p3;
c=0.47;r_keha=2.;gamma=6.67408e-11; m_maa=5.97219e24;
p_atm=101000.;mu_atm=29;
R=8314.;t_atm=300.;
m_keha=2800;r_maa=6378.0; 

s_keha=pi*r_keha*r_keha;
ro_atm_0=p_atm*mu_atm/(R*t_atm);
p1=m_maa*gamma*1e-9;
p2=c*s_keha*500./m_keha;
p3=mu_atm*9814./(R*t_atm);
axis equal
angle=0.5;
v0=7.78843;

vx=-v0*sind(angle);
vy= v0*cosd(angle);
y0=[6578.;0.;vx;vy];

time=0:1:200000;
options = odeset('RelTol',1e-9,'AbsTol',1e-9,'MaxStep',1.,'InitialStep',0.01,'Events',@mytest);
[t,y]=ode45(@keplerf,time,y0,options);

rr=sqrt(power(y(:,1),2)+power(y(:,2),2));axis equal
v_abs=sqrt(power(y(:,3),2)+power(y(:,4),2));
h=rr-r_maa;
ro=ro_atm_0*exp(-p3*h);
ax=-y(:,1)*p1./power(rr,3)-p2.*ro.*v_abs.*y(:,3);
ay=-y(:,2)*p1./power(rr,3)-p2.*ro.*v_abs.*y(:,4);
aabs=sqrt(power(ax,2)+power(ay,2))*1000.;

subplot(1,3,1)

plot(h,v_abs*1000)
xlabel ('korgus,km')
ylabel ('kiirus, m/s')
grid on

subplot(1,3,2)
hold on
plot(h,aabs)
xlabel ('korgus,km')
ylabel ('kiirendus,m/s2')
grid on

subplot(1,3,3)
plot(t/60,h)
xlabel ('aeg')
ylabel ('korgus,km')
grid on


function dy=keplerf(t,y)
global r_maa ro_atm_0 p1 p2 p3;;
ra=sqrt(y(1)^2+y(2)^2);
v_abs=sqrt(y(3)^2+y(4)^2);
h=ra-r_maa;
ro=ro_atm_0*exp(-p3*h);
dy=zeros(4,1);
dy(1)=y(3);
dy(2)=y(4);
dy(3)=-y(1)*p1/ra^3-p2*ro*v_abs*y(3);
dy(4)=-y(2)*p1/ra^3-p2*ro*v_abs*y(4);
end

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