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