rcksim.m

来自「航天工程工具箱」· M 代码 · 共 302 行

function [flight_,t,x]=rcksim(name,eng,Cd,Dr,mr,dt)
%RCKSIM  Simulation of rocket trajectory.
%   [DATA,T,X] = RCKSIM(NAME,ENG,CD,DR,MR,DT) simulates the trajectory
%   of a multi staged rocket with name NAME, having for each stage
%   a motor ENG given by LOADENG, an overall "static" subsonic drag
%   coefficient CD (for current stage and upper stages),
%   a body diameter DR, an empty mass MR (without the motor) and a
%   separation delay (burnout-separation) DT. It is assumed that the delay
%   for separation-ignition equals zero.
%
%   [...] = RCKSIM(ROCKET) does the same as above but for a rocket
%   definition structure ROCKET given by RCKTOOL.
%
%   LENGTH(DT)+1 = LENGTH(MR) = LENGTH(DR) = LENGTH(CD) = LENGTH(ENG)
%                = # stages
%
%   The returned structure DATA contains information about the
%   simulation results. The format for DATA will be:
%
%      DATA.name         : Name of rocket
%      DATA.Cd           : Drag coefficient for each stage
%      DATA.Dr           : Diameter of each stage
%      DATA.mass.mr      : Mass of each stage w/o motor
%      DATA.mass.mf      : Fuel mass for each stage
%      DATA.mass.mm      : Motor mass for each stage
%      DATA.mass.m[0-3]  : Total mass of rocket at event
%      DATA.alt.h[0-3]   : Altitude at event
%      DATA.vel.ve       : Exhaust velocity for each stage
%      DATA.vel.v[0-3]   : Velocity at event
%      DATA.acc.amax     : Maximum acceleration
%      DATA.acc.a[0-3]   : Acceleration at event
%
%   where an event is associated with each index 0-3 as:
%      0 : ignition for each stage
%      1 : burn-out for each stage
%      2 : turning point (maximum altitude)
%      3 : impact
%
%   T is the sampled time points and X is:
%
%      X(:,1) = altitude
%      X(:,2) = velocity
%      X(:,3) = acceleration
%      X(:,4) = mass
%
%   See also RCKMOD, LOADENG, RCKTOOL, MOTORCLASS.

% Copyright (c) 2002-02-09, B. Rasmus Anthin.
% Revisited 2003-06-11,
%           2003-06-14 - 2003-06-15,
%           2003-06-18 - 2003-06-19.

g=constant('g');

if nargin==1 & isempty(name), return,end

if isstruct(name)
   r=name;
   name=r.name;
   slen=r.stages;
   eng=r.eng;
   Cd=r.Cd;
   Dr=r.Dr;
   mr=r.mr;
   dt=r.dt;
else
   slen=length(eng);
end

%name=[];
for i=1:slen
%   name=[name ' - ' eng(i).name];
   mf(i)=eng(i).mf;
   mm(i)=eng(i).mm;
end
%name=name(4:end);

for i=1:slen
   m0(i)=sum(mr(i:end))+sum(mm(i:end))+sum(mf(i:end));   %start
   m1(i)=m0(i)-mf(i);                                    %burnout   
end

t0=0;t1=eng(1).thrust(end,1);
for i=2:slen
   t0(i)=t1(i-1)+dt(i-1);                                %start
   t1(i)=t0(i)+eng(i).thrust(end,1);                     %burnout
end

for i=1:slen
   Ar(i)=pi*max(Dr(i:end))^2/4;

   thrust(i).t=[t0(i);t0(i)+eng(i).thrust(:,1)];
   thrust(i).F=[0;eng(i).thrust(:,2)];

   ve(i)=trapz(thrust(i).t,thrust(i).F)/mf(i);
end

atmmodel='US';    %alter to change to another model.
switch(atmmodel)
case 'MSIS-E-90'
   %MSIS-E-90 (up to 1000 km)
   %http://nssdc.gsfc.nasa.gov/space/model/models/msis_n.html#height
   load atmdata
   atmos.h=atm(:,1)*1e3; %km -> m.
   atmos.rho=atm(:,5)*1e3; %g/cm^3 -> kg/m^3.
   atmos.T=atm(:,6);
case 'US'
   %US Standard Atmosphere (up to 65 km)
   %http://dss.ucar.edu/docs/equations/std.atmos.html
   stdatm
   atmos.h=atm(:,2);
   atmos.rho=atm(:,5);
   atmos.T=convert(atm(:,4),'oC','K'); %K -> oC
end

tic
optn=odeset('refine',2,'reltol',1e-6,'abstol',1e-5,'stats','off','events','on');
tt=[]; xx=[];
tj=NaN; xj=NaN*ones(1,3);
init=[0 0];
for i=1:slen
   %IGNITION #I -> BURN-OUT #I
   [t,x]=ode15s('rckmod',[t0(i) t1(i)],[init m0(i)],optn,thrust(i),atmos,Ar(i),Cd(i),m1(i),ve(i));
   tt=[tt;t]; xx=[xx;x];
   init=x(end,1:2);
   if i<slen
      if t1(i)<t0(i+1)          %BURN-OUT #I -> IGNITION #I+1
         [t,x]=ode15s('rckmod',[t1(i) t0(i+1)],[init m1(i)],optn,thrust(i),atmos,Ar(i),Cd(i),m1(i),ve(i));
         tt=[tt;t]; xx=[xx;x];
         init=x(end,1:2);
      end
      null.t=0; null.F=0;       %BURN-OUT #I -> IMPACT #I
      [t,x]=ode15s('rckmod',[t0(i+1) inf],[init mr(i)+mm(i)],optn,null,atmos,Ar(i),Cd(i),m1(i),ve(i));
      tj=[tj;t;NaN]; xj=[xj;x;NaN*ones(1,3)];
   end
end
%BURN-OUT #I+N -> IMPACT #I+N
[t,x]=ode15s('rckmod',[t(end) inf],x(end,:),optn,thrust(end),atmos,Ar(end),Cd(end),m1(end),ve(end));
tt=[tt;t]; xx=[xx;x];
t=tt;x=xx;
tsim=toc;

for i=1:slen
   i0(i)=min(find(t==t0(i)));  %ignition
   i1(i)=min(find(t==t1(i)));  %burn-out
end
i2=min(find(x(:,1)==max(x(:,1))));      %coast
i3=length(t);

t0=t(i0);
t1=t(i1);
t2=t(i2);
t3=t(i3);
h0=x(i0,1);
h1=x(i1,1);
h2=x(i2,1);
h3=x(i3,1);
v0=x(i0,2);
v1=x(i1,2);
v2=x(i2,2);
v3=x(i3,2);

m=x(:,3);
f=-gradient(x(:,3),t);
a=gradient(x(:,2),t);
x(:,3)=a;
x(:,4)=m;
m=x(:,4);    % mass
a=x(:,3);    % acceleration
v=x(:,2);    % velocity
h=x(:,1);    % height

gh=figure;set(gh,'name',name,'numbertitle','off')
subplot(4,1,1), hold on
plot(tj,xj(:,1),'g')
plot(t,x(:,1))
ax=axis; axis tight
events(t0,t1,t2,ax)
xlabel t,ylabel h
grid, set(gca,'xgrid','off')

subplot(4,1,2), hold on
plot(tj,xj(:,2),'g')
plot(t,x(:,2))
ax=axis; axis tight
events(t0,t1,t2,ax)
xlabel t,ylabel('h''')
grid, set(gca,'xgrid','off')

subplot(4,1,3), hold on
%plot(tj,gradient(xj(:,2)),'g')
plot(t,x(:,3))
ax=axis; axis tight
events(t0,t1,t2,ax)
xlabel t,ylabel('h''''')
grid, set(gca,'xgrid','off')

subplot(4,1,4), hold on
plot(t,x(:,4))
ax=axis; axis tight
events(t0,t1,t2,ax)
xlabel t,ylabel m
grid, set(gca,'xgrid','off')

   disp(' ')
   disp(['Simulation time : ' num2str(tsim) ' s'])
   disp(['# samples       : ' num2str(size(x,1))])

   disp(' ')
   disp( 'MASS---------------------------------')
for i=1:slen
   disp(['  Rocket     ' int2str(i) '    : ' num2mass(mr(i),6)])
   disp(['  Propellant ' int2str(i) '    : ' num2mass(mf(i),6)])
   disp(['  Motor      ' int2str(i) '    : ' num2mass(mm(i),6)])
end
for i=1:slen
   disp(['  Ignition   ' int2str(i) '    : ' num2mass(m0(i),6)])
   disp(['  Burn-out   ' int2str(i) '    : ' num2mass(m1(i),6)])
end
   disp(' ')
   disp( 'ALTITUDE-----------------------------')
for i=1:slen
   disp(['  Ignition   ' int2str(i) '    : ' num2sci(h0(i)) 'm'])
   disp(['  Burn-out   ' int2str(i) '    : ' num2sci(h1(i)) 'm'])
end
   disp(['  Turning         : ' num2sci(h2) 'm'])
   disp(' ')
   disp( 'VELOCITY-----------------------------')
for i=1:slen
   disp(['  Exit       ' int2str(i) '    : ' num2sci(ve(i)) 'm/s'])
end
for i=1:slen
   disp(['  Ignition   ' int2str(i) '    : ' num2sci(v0(i)) 'm/s'])
   disp(['  Burn-out   ' int2str(i) '    : ' num2sci(v1(i)) 'm/s'])
end
   disp(['  Impact          : ' num2sci(v3) 'm/s'])
   disp(' ')
   disp( 'ACCELERATION-------------------------')
   disp(['  Maximum         : ' num2sci(max(abs(a))/g) 'g'])
for i=1:slen
   disp(['  Ignition   ' int2str(i) '    : ' num2sci(a(i0(i))/g) 'g'])
   disp(['  Burn-out   ' int2str(i) '    : ' num2sci(a(i1(i))/g) 'g'])
end
   disp(['  Turning         : ' num2sci(a(i2)/g) 'g'])
   disp(['  Impact          : ' num2sci(a(i3)/g) 'g'])
   disp(' ')
   disp( 'TIME---------------------------------')
for i=1:slen
   tm0=sec2hms(t0(i));
   tm1=sec2hms(t1(i));
   disp(['  Ignition   ' int2str(i) '    : ' num2str(tm0(1)) ' h, ' num2str(tm0(2)) ' m, ' num2str(tm0(3)) ' s'])
   disp(['  Burn-out   ' int2str(i) '    : ' num2str(tm1(1)) ' h, ' num2str(tm1(2)) ' m, ' num2str(tm1(3)) ' s'])
end
   tm2=sec2hms(t2);
   tm3=sec2hms(t3);
   disp(['  Turning         : ' num2str(tm2(1)) ' h, ' num2str(tm2(2)) ' m, ' num2str(tm2(3)) ' s'])
   disp(['  Impact          : ' num2str(tm3(1)) ' h, ' num2str(tm3(2)) ' m, ' num2str(tm3(3)) ' s'])
   disp(' ')
   disp( 'OTHER--------------------------------')
for i=1:slen
   disp(['  Diameter   ' int2str(i) '    : ' num2sci(Dr(i)) 'm'])
end
   disp(['  Impact energy   : ' num2sci(m1(end)*v3^2/2) 'J'])
   disp(' ')
   
flight.name=name;
flight.Cd=Cd;
flight.Dr=Dr;
flight.mass.mr=mr;
flight.mass.mf=mf;
flight.mass.mm=mm;
flight.mass.m0=m0;
flight.mass.m1=m1;
flight.mass.m2=m(i2);
flight.mass.m3=m(i3);
flight.alt.h0=h0;
flight.alt.h1=h1;
flight.alt.h2=h2;
flight.alt.h3=h3;
flight.vel.ve=ve;
flight.vel.v0=v0;
flihgt.vel.v1=v1;
flight.vel.v2=v2;
flight.vel.v3=v3;
flight.acc.max=max(a);
flight.acc.a0=a(i0);
flight.acc.a1=a(i1);
flight.acc.a2=a(i2);
flight.acc.a3=a(i3);


if nargout, flight_=flight;end


%%%%%%%%%%%%%%%%%%%%

function events(t0,t1,t2,ax)
for i=1:length(t0)
   plot([t0(i) t0(i)],ax(3:4),'r:')
   plot([t1(i) t1(i)],ax(3:4),'m:')
end
plot([t2 t2],ax(3:4),'b:')