global  rho_atm_0 p1 p2 p3 G m_maa r_maa 
r_keha=1.8;
s_keha=pi*r_keha^2;
c_keha=0.47;
v_keha=4.*pi*r_keha^3/3;
m_keha=15.;
p_atm=100000.;
t_atm=293.;
mu_atm=29.;
r_gaas=8314.;
m_maa=5.97219e+24;
G=6.6743015e-11;
r_maa=6371000.;

g0=rask(G,m_maa,r_maa,0.);

rho_atm_0=p_atm*mu_atm/(r_gaas*t_atm);
p1=mu_atm*g0/(r_gaas*t_atm);
p2=0.5*c_keha*s_keha/m_keha;
p3=v_keha/m_keha;
time=0:0.1:4000;

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

g_rask=rask(G,m_maa,r_maa,y(:,1));
rho=rho_atm_0*exp(-p1*y(:,1)./(1.+y(:,1)/r_maa));
v=abs(y(:,2));
ar=-g_rask;
aa=g_rask.*rho*p3;
at=-p2*rho.*v.*y(:,2);
atot=ar+aa+at;


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

subplot(2,2,3)
plot(t,[ar,aa,at,atot])
legend(["ar","aa","at","atot"]) 
grid on

function yp = func(t,y)
global  rho_atm_0 p1 p2 p3 G m_maa r_maa 
yp=zeros(2,1);
yp(1)=y(2);
g_rask=rask(G,m_maa,r_maa,y(1));
rho=rho_atm_0*exp(-p1*y(1)/(1.+y(1)/r_maa));
v=abs(y(2));
ar=-g_rask;
aa=g_rask*rho*p3;
at=-p2*rho*v*y(2);
yp(2)=ar+aa+at;
end

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

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