clear
global g_maa r_maa ro_ohk_0 param1 param2 c_aero m_keha;
c_aero=0.47;r_keha=2.;
p_ohk=102500.;mu_ohk=29;
r_gaas=8314.;t_ohk=300.;
m_keha=2800;r_maa=6378.0; 
g_maa=3.986004418e+5;

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.69;
v0=7.78843;

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

time=0:1:7000;
options = odeset('RelTol',1e-6,'AbsTol',1e-6,'MaxStep',0.1,'InitialStep',0.001,'Event',@mycheck);
[t,y,te,ye,ie]=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;

for i=1:length(h)
if h(i) >=5    
param(i,1)=c_aero*s_keha*500./m_keha;
else
param(i,1)=c_aero*s_keha*500.*100./m_keha;  
end
end

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./rr.^3-param.*ro_ohk.*v_abs.*y(:,3);
ay=-y(:,2)*g_maa./rr.^3-param.*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);

figure ('Position',[100 100 1000 500])

subplot('Position',[0.1 0.1 0.2 0.6])
plot(xx,yy,'black')
axis([-1.5 1.5 -1.5 1.5])
axis equal

hold on
plot(y(:,1)./r_maa,y(:,2)./r_maa,'Color','red')
grid on
itime=100:10:length(ax);
quiver(y(itime,1)./r_maa,y(itime,2)./r_maa,ax(itime),ay(itime),0.5,'Color','black')

subplot('Position',[0.4 0.1 0.2 0.8])
%plot(h,aabs)
%xlabel ('korgus,km')
%ylabel ('kiirendus, m/s^2')
plot(t,v_abs)
grid on

subplot('Position',[0.7 0.1 0.2 0.8])
plot(h,v_abs)
xlabel ('korgus,km')
ylabel ('kiirus, km/s')
grid on

disp(['max kiirendus=',num2str(max(aabs(find(h>5.5))) ),' m/s^2=',num2str(max(aabs(find(h>5.5)))/9.814),'g'])
disp(['min kiirus=',num2str(min(v_abs)*1000.),' m/s'])
disp(['aeg=',num2str(te/60),' min'])


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

function [h,finish,dir]= mycheck(t,y)
global r_maa param2 c_aero m_keha
distance=sqrt(y(1)^2+y(2)^2)-r_maa;

h=[distance];
finish=[1];
dir=[0];

if distance < 5.
r_keha=20.;
s_keha=pi*r_keha^2;
param2=c_aero*s_keha*500./m_keha;
end

end
