global  g_rask ro_atm_0 p1 p2 p3
r_keha=0.05;
C=0.47;
g_rask=9.814;
m_keha=0.5;
p_atm=100400.;
mu_atm=29.;
R=8314.;
t_atm=280.;

s_keha=pi*r_keha^2;
v_keha=4.*pi*r_keha^3/3.;
ro_atm_0=p_atm*mu_atm/(R*t_atm);
p1=mu_atm*g_rask/(R*t_atm);
p2=0.5*C*s_keha/m_keha;
p3=g_rask*v_keha/m_keha;

y0=[0.;0.;10.;10.];
options=odeset('RelTol',1.e-8,'AbsTol',1.e-8,'MaxStep',0.001,'InitialStep',0.0011,'Event',@mycheck) ;
[t,y,te,ye,ie]=ode45(@dim2_rof,[0 5],y0,options);
te
ye
ie
subplot (2,1,1)
plot(y(:,1),y(:,2))
grid on
subplot(2,1,2)
plot(t(:),sqrt(y(:,3).^2+y(:,4).^2))

grid on

function dy = dim2_rof(t,y)
global g_rask ro_atm_0 p1 p2 p3;
dy=zeros(4,1);
v_abs=sqrt(y(3)*y(3)+y(4)*y(4));
ro=ro_atm_0*exp(-p1*y(2));
dy(1)=y(3);
dy(2)=y(4);
dy(3)=-p2*v_abs*y(3);
dy(4)=-g_rask-p2*v_abs*y(4)+p3*ro;
end

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

