global g_maa r_maa ro_ohk_0 param1 param2;
c_aero=0.47;r_keha=1.1;gamma=6.67408e-11; m_maa=5.97219e24;
p_ohk=102500.;mu_ohk=29;
r_gaas=8314.;t_ohk=300.;
m_keha=2800;r_maa=6371.0; 
g_maa=m_maa*gamma*1e-9;

s_keha=pi*r_keha*r_keha;
ro_ohk_0=p_ohk*mu_ohk/(r_gaas*t_ohk);
param1=mu_ohk*9814./(r_gaas*t_ohk);
param2=c_aero*s_keha*500./m_keha;

angle=0.5;
v0=7.78843;

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

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

rr=sqrt(power(y(:,1),2)+power(y(:,2),2));
v_abs=sqrt(power(y(:,3),2)+power(y(:,4),2));
h=rr-r_maa;
ro_ohk=ro_ohk_0*exp(-param1*h);
ax=-y(:,1)*g_maa./power(rr,3)-param2.*ro_ohk.*v_abs.*y(:,3);
ay=-y(:,2)*g_maa./power(rr,3)-param2.*ro_ohk.*v_abs.*y(:,4);
aabs=sqrt(power(ax,2)+power(ay,2))*1000.;
phi=0:1:360;
xx=sind(phi);
yy=cosd(phi);

subplot(1,3,1)
pbaspect([1 1 1])
hold on
axis([-2 2 -2 2])
plot(xx,yy,'black')
plot(y(:,1)./r_maa,y(:,2)./r_maa,'r')
grid on

subplot(1,3,2)
plot(h,aabs)
xlabel ('korgus,km')
ylabel ('kiirendus, m/s^2')
pbaspect([1 1 1])
grid on

subplot(1,3,3)

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

disp(['max kiirendus=',num2str(max(aabs)),' m/s^2'])
disp(['min kiirus=',num2str(min(v_abs)*1000.),' m/s'])

function dy=keplerf(t,y)
global g_maa r_maa ro_ohk_0 param1 param2 ;
ra=sqrt(y(1)^2+y(2)^2);
v_abs=sqrt(y(3)^2+y(4)^2);
h=ra-r_maa;
ro_ohk=ro_ohk_0*exp(-param1*h);
dy=zeros(4,1);
dy(1)=y(3);
dy(2)=y(4);
dy(3)=-y(1)*g_maa/ra^3-param2*ro_ohk*v_abs*y(3);
dy(4)=-y(2)*g_maa/ra^3-param2*ro_ohk*v_abs*y(4);

end
