global g_const r_maa m_maa p1 p2 p3 ro_atm_0 vtuul

pi=3.1415926536;
r_keha=0.5;
c=0.47;
m_keha=5.0;
m_maa=5.97219e+24;
g_const=6.6743015e-11;
r_maa=6371000.;
p_atm=100400.;
t_atm=280.;
mu_atm=29.;
r_gaas=8314.;
g_rask=rask(g_const,m_maa,r_maa,0.)
s_keha=pi*r_keha^2;
v_keha=4.*pi*r_keha^3/3;
ro_atm_0=p_atm*mu_atm/(r_gaas*t_atm);
p1=g_rask*mu_atm/(r_gaas*t_atm);
p2=0.5*c*s_keha/m_keha;
p3=v_keha/m_keha;
vtuul=-12.5;

nt=40000;
dt=0.001;
time=0:dt:nt*dt;
y=[0;0;10;10];

options=odeset('RelTol',1.e-8,'AbsTol',1.e-8,'MaxStep',0.01,'InitialStep',0.01,'Event',@mycheck) ;
[t,y]=ode45(@func,time,y0,options);

vabs_t=sqrt((y(:,3)-vtuul).^2+y(:,4).^2);
vabs=sqrt(y(:,3).^2+y(:,4).^2);

ro=ro_atm_0*exp(-p1*y(:,2)./(1+y(:,2)/r_maa));
g_rask=rask(g_const,m_maa,r_maa,y(:,2));
ax=-p2.*ro.*vabs_t.*(y(:,3)-vtuul);
ay=-g_rask-p2*ro.*vabs_t.*y(:,4)+p3*g_rask.*ro;

aabs=sqrt(ax.^2+ay.^2);
atau=(y(:,3).*ax+y(:,4).*ay)./vabs;
anorm=sqrt(aabs.^2-atau.^2);
radius=vabs.^2./anorm;

ekin=0.5*m_keha*vabs.^2;
epot=m_keha*g_rask.*y(:,2);
etot=ekin+epot;

subplot(1,3,1)
plot(y(:,1),y(:,2))
grid on


subplot(1,3,2)
plot(t,[ekin,epot,etot])
grid on


subplot(1,3,3)
plot(t,[anorm,atau,radius])
grid on


function [dy]=func(t,y)
global g_const r_maa m_maa p1 p2 p3 ro_atm_0 vtuul;
dy=zeros(4,1);
dy(1)=y(3);
dy(2)=y(4);
g_rask=rask(g_const,m_maa,r_maa,y(2));
ro=ro_atm_0*exp(-p1*y(2)/(1+y(2)/r_maa));
vabs=sqrt((y(3)-vtuul)^2+y(4)^2);
dy(3)=-p2*ro*vabs*(y(3)-vtuul);
dy(4)=-g_rask-p2*ro*vabs*y(4)+p3*g_rask*ro;
end

function g_rask=rask(g_const,m_maa,r_maa,h)
g_rask=g_const*m_maa./(r_maa+h).^2;
end

function [h,finish,dir]= mycheck(t,y)
h=[y(2)];
finish=[1]; 
dir=[0];
end
