rckmod.m

航天工程工具箱

function [out1,out2,out3]=rckmod(t,x,flag,thrust,atmos,Ar,Cd,m1,ve)
% x1'=x2
% x2'=a_r(t)-a_g(y1)-a_d(y2)
% x3'=-f(t)
%
% Refer to RCKSIM for details.

% Copyright (c) 2003-06-11, B. Rasmus Anthin.

%Model limitations:
%------------------
%Assumes that the exhaust velocity is constant, even when the thrust isn't.
%Assumes Pe=P0 (no net pressure).
%Constant Cd.
%Constant temperature at h=0 (density profile for a certain time at a certain place).


mu=constant('Me')*constant('G');
R=constant('Ree');

h=x(1);
v=x(2);
m=x(3);
 %gravity
g=mu/(h+R)^2;
 %air density, air temperature
rho=interp1(atmos.h,atmos.rho,x(1),'spline');
 %thrust
t1=thrust.t(end);
if t<=t1
   Fr=interp1(thrust.t,thrust.F,t,'linear');
   f=Fr/ve*(m>m1);
   ar=Fr/m;
else
   f=0;
   ar=0;
end

ag=g;                          %gravitational force
ag=ag;

ad=rho*Ar*Cd*v^2/(2*m);        %drag force
ad=ad*sign(v)*(h<2e6);

%%%%%%%% ODE %%%%%%%%

dx=[v
    ar-ag-ad
    -f];
 
%%%%%%%%%%%%%%%%%%%%%
%Events
 
if strcmp(flag,'events')
   %out1=[t-t1 h];
   %out2=[1 1];
   %out3=[1 -1];
   out1=h;
   out2=1;
   out3=-1;
else
   out1=dx;
end