clear

global  gm g_const r_maa m_maa ro_atm_0 p1 p2 p3
c_keha=0.47;r_keha=2.0;
m_keha=2800.;p_atm=101000.;
mu_atm=29.;r_gaas=8314.;
t_atm=300.;
g_const=6.67408e-11;
m_maa=5.9722e+24;
gm=g_const*m_maa*1.e-9;
r_maa=6371.;

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);
g_rask=rask(g_const,r_maa,m_maa,0.);
p1=mu_atm*g_rask*1.e3/(r_gaas*t_atm);
p2=s_keha*c_keha*500./m_keha;
p3=v_keha*1.e-3/m_keha;

%v0=7.7884079;
v0=10.;
alpha=5.;
h0=150.;
      
y0=[r_maa+h0;0;-v0*sin(pi*alpha/180.);v0*cos(pi*alpha/180.)];

options = odeset('RelTol',1e-9,'AbsTol',1e-9,'MaxStep',1000.,'InitialStep',1.,'Events',@mycheck);
[t,res]=ode45(@func,[0,400000], y0,options);

vv=sqrt(res(:,3).^2+res(:,4).^2);
rr=sqrt(res(:,1).^2+res(:,2).^2);
h=rr-r_maa;

ro=ro_atm_0*exp(-p1*h./(1+h./r_maa));
g_raskus=rask(g_const,r_maa,m_maa,h);
a(:,1)=-res(:,1)*gm./rr.^3-p2*ro.*vv.*res(:,3)+p3*g_raskus.*ro.*res(:,1)./rr;
a(:,2)=-res(:,2)*gm./rr.^3-p2*ro.*vv.*res(:,4)+p3*g_raskus.*ro.*res(:,2)./rr;
aa=sqrt(a(:,1).^2+a(:,2).^2);

f1=figure;
f2=figure;
clf
figure(f1)
subplot(2,2,1)
plot(t,h)
legend("h(t)")
grid on
subplot(2,2,2)
plot(t,aa*1000)
legend("a(t)")
grid on
subplot(2,2,3)
plot(t,vv)
legend("v(t)")
grid on
subplot(2,2,4)
plot(t,res(:,1))
legend("x(t)")
grid on
%--------------
figure(f2)
axis equal
hold on
th = 0:pi/50:2*pi;
xmaa =  cos(th);
ymaa = sin(th);
plot(xmaa, ymaa);

nt=size(t);
p1=animatedline();
for i=1:nt
addpoints(p1,res(i,1)/r_maa,res(i,2)/r_maa);
pause(0.02)
end
%------------------

t(length(t))/60
t(length(t))/3600


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

function g_rask=rask(g_const,r_maa,m_maa,h)
g_rask=g_const*m_maa*1.e-6./(r_maa+h).^2;
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
