implicit none
real y(4),work(27),mu_atm,m_keha,a(4),ar(2),at(2),aa(2)
real g_const,r_maa,m_maa,p1,p2,p3,ro_atm_0,vtuul,t,tout,dt
real pi,r_keha,c,s_keha,v_keha,p_atm,relerr,abserr
real vabs,aabs,atau,anorm,radius,r_gaas,ekin,epot,etot,v0,alpha
real g_rask,t_atm,ro
integer it,nt,iflag, iwork(5),neqn
common g_const,r_maa,m_maa,p1,p2,p3,ro_atm_0,vtuul
external f
open(10,file="dim2.dat")
open(20,file="kiirendus.dat")
pi=3.1415926536
r_keha=0.5
c=0.47
m_keha=5.
m_maa=5.97219e+24
g_const=6.6743015e-11
r_maa=6371000.
p_atm=100400.
t_atm=280.
mu_atm=29.
r_gaas=8314.
 call rask(g_const,m_maa,r_maa,0.,g_rask)
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)
p1=g_rask*mu_atm/(r_gaas*t_atm)
p2=0.5*c*s_keha/m_keha
p3=g_rask*v_keha/m_keha
neqn=4
nt=10000
dt=0.001
t=0.
relerr=1.e-6
abserr=1.e-6

alpha=45
v0=14.142
vtuul=-12.5

y(1)=0.
y(2)=0.
y(3)=v0*sin(alpha*pi/180.)
y(4)=v0*cos(alpha*pi/180.)
iflag=1

do it=1,nt
tout=t+dt

call rkf45(f,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)
if (iflag /= 2)then
iflag=2
endif
if (y(2) <0.) then
stop
endif

ekin=0.5*m_keha*vabs**2
epot=m_keha*g_rask*y(2)
etot=ekin+epot

call f(t,y,a)
vabs=sqrt(y(3)*y(3)+y(4)*y(4))
aabs=sqrt(a(3)*a(3)+a(4)*a(4))
atau=(y(3)*a(3)+y(4)*a(4))/vabs
anorm=sqrt(aabs**2-atau**2)
radius=vabs**2/anorm

ro=ro_atm_0*exp(-p1*y(2))
ar(1)=0.
ar(2)=-g_rask
at(1)=-p2*ro*vabs*(y(3)-vtuul)
at(2)=-p2*ro*vabs*y(4)
aa(1)=0.
aa(2)=ro*p3

write (10,'(15f15.6)')tout,y,ar(1:2),at(1:2),aa(1:2),a(3:4)
write (20,*) atau,anorm,aabs,radius
enddo

stop
end

subroutine f(t,y,yp)
implicit none
real y(4),yp(4),g_const,r_maa,m_maa,p1,p2,p3,ro_atm_0,vtuul
real ro,vabs,t,g_rask
common g_const,r_maa,m_maa,p1,p2,p3,ro_atm_0,vtuul
yp(1)=y(3)
yp(2)=y(4)
call rask(g_const,m_maa,r_maa,y(2),g_rask)
ro=ro_atm_0*exp(-p1*y(2))
vabs=sqrt((y(3)-vtuul)**2+y(4)**2)
yp(3)=-p2*ro*vabs*(y(3)-vtuul)
yp(4)=-g_rask-p2*ro*vabs*y(4)+p3*ro
return
end

subroutine rask(g_const,m_maa,r_maa,h,g_rask)
implicit none
real g_const,m_maa,r_maa,g_rask,h
g_rask=g_const*m_maa/(r_maa+h)**2
return
end


--------------------------------------------------------------------------
a=load('dim2.dat');
b=load('kiirendus.dat');
hold on
tindex=1:100:length(a(:,1));
clf
max(a(:,3))
min(b(:,4))

subplot(1,2,1)

hold on
plot (a(:,2),a(:,3),'Color','blue')
quiver(a(tindex,2),a(tindex,3),a(tindex,6),a(tindex,7))
quiver(a(tindex,2),a(tindex,3),a(tindex,8),a(tindex,9))
quiver(a(tindex,2),a(tindex,3),a(tindex,10),a(tindex,11))
quiver(a(tindex,2),a(tindex,3),a(tindex,12),a(tindex,13))
legend('tr','ag','at','aa','atot')
axis equal
grid on

subplot(1,2,2)
hold on
plot (a(:,1),b(:,1))
plot (a(:,1),b(:,2))
plot (a(:,1),b(:,3))
plot (a(:,1),b(:,4))
legend('at','an','atot','r')
grid on

---------------------------------------------------------------------------

global  ro_keha ro_ohk_0 g_rask param1 param2
r_keha=0.5;
c_aero=0.47;
g_rask=9.814;
m_keha=5.;
p_ohk=102500.;
mu_ohk=29.;
r_gaas=8314.;
t_ohk=266.;

s_keha=pi*r_keha^2;
v_keha=4.*pi*r_keha^3/3.;
ro_ohk_0=p_ohk*mu_ohk/(r_gaas*t_ohk);
ro_keha=m_keha/v_keha;
param2=mu_ohk*g_rask/(r_gaas*t_ohk);
param1=0.5*c_aero*s_keha/m_keha;

y0=[0.,0.,10.,20.];
options=odeset('RelTol',1.e-8,'AbsTol',1.e-8,'MaxStep',0.01,'InitialStep',0.01,'Event',@mycheck) ;
[t,y,te,ye,ie]=ode45(@dim2_rof,[0 5],y0,options);
%*** kiirenduse arvutused ****
v_abs=sqrt(y(:,3).*y(:,3)+y(:,4).*y(:,4));
ro_ohk=ro_ohk_0*exp(-param2*y(:,2));
ax=-param1*ro_ohk.*v_abs.*y(:,3);
ay=g_rask*(ro_ohk/ro_keha-1.)-param1*ro_ohk.*v_abs.*y(:,4);
nt=length(ax);
time=1:100:nt;
%energy
ekin=0.5*m_keha*(y(:,3).^2+y(:,4).^2);
epot=m_keha*g_rask*y(:,2);
etot=ekin+epot;


subplot (2,1,1)
hold on
plot(y(:,1),y(:,2))
quiver(y(time,1),y(time,2),ax(time),ay(time),'Color','red')
quiver(y(time,1),y(time,2),y(time,3),y(time,4),'Color','black')
legend('traj','a','v')
grid on
subplot(2,1,2)
plot(t,ekin,t,epot,t,etot)
legend('ekin','epot','etot')
grid on

function dy = dim2_rof(t,y)
global ro_keha ro_ohk_0 g_rask param1 param2;
dy=zeros(4,1);
v_abs=sqrt(y(3)*y(3)+y(4)*y(4));
ro_ohk=ro_ohk_0*exp(-param2*y(2));
param3=g_rask*(ro_ohk/ro_keha-1.);
param4=param1*ro_ohk;
dy(1)=y(3);
dy(2)=y(4);
dy(3)=-param4*v_abs*y(3);
dy(4)=param3-param4*v_abs*y(4);
end

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