【气动学】基于龙格库塔算法实现外弹道仿真含Matlab源码

function comet3m(varargin)%COMET3 3-D Comet-like trajectories.%   COMET3(Z) displays an animated three dimensional plot of the vector Z.%   COMET3(X,Y,Z) displays an animated comet plot of the curve through the%   points [X(i),Y(i),Z(i)].%   COMET3(X,Y,Z,p) uses a comet of length p*length(Z). Default is p = 0.1.%%   COMET3(AX,...) plots into AX instead of GCA.%%   Example:%       t = -pi:pi/500:pi;%       comet3(sin(5*t),cos(3*t),t)%%   See also COMET.%   Charles R. Denham, MathWorks, 1989.%   Revised 2-9-92, LS and DTP; 8-18-92, 11-30-92 CBM.%   Copyright 1984-2006 The MathWorks, Inc. %   $Revision: 5.11.4.4 $  $Date: 2011/03/09 07:03:36 $% Parse possible Axes input[ax,args,nargs] = axescheck(varargin{:});error(nargchk(1,10,nargs,'struct'));% Parse the rest of the inputsif nargs < 2, x = args{1}; endif nargs == 2, y = args{2}; endif nargs < 3, z = x; x = 1:length(z); y = 1:length(z); endif nargs == 3, [x,y,z] = deal(args{:}); endif nargs < 4, p = 0.10; endif nargs == 4, [x,y,z,p] = deal(args{:}); endif nargs ==10, [x,y,z,vx,vy,vz,ph,fy,hg,p] = deal(args{:}); end if ~isscalar(p) || ~isreal(p) || p < 0 || p >= 1    error(message('MATLAB:comet3:InvalidP'));endax = newplot(ax);if ~ishold(ax),  [minx,maxx] = minmax(x);  [miny,maxy] = minmax(y);  [minz,maxz] = minmax(z);  axis(ax,[minx maxx miny maxy minz maxz])endco = get(ax,'colororder');grid onif size(co,1)>=3,  % Choose first three colors for head, body, and tail  head = line('parent',ax,'color',co(1,:),'marker','.','MarkerSize',20,'erase','xor', ...              'xdata',x(1),'ydata',y(1),'zdata',z(1));  body = line('parent',ax,'color',co(2,:),'linestyle','-','erase','none', ...              'xdata',[],'ydata',[],'zdata',[]);  tail = line('parent',ax,'color',co(3,:),'linestyle','-','erase','none', ...              'xdata',[],'ydata',[],'zdata',[]);else  % Choose first three colors for head, body, and tail  head = line('parent',ax,'color',co(1,:),'marker','.','MarkerSize',20,'erase','xor', ...              'xdata',x(1),'ydata',y(1),'zdata',z(1));  body = line('parent',ax,'color',co(1,:),'linestyle','--','erase','none', ...              'xdata',[],'ydata',[],'zdata',[]);  tail = line('parent',ax,'color',co(1,:),'linestyle','-','erase','none', ...              'xdata',[],'ydata',[],'zdata',[]);endmm = length(z);k = round(p*mm);TIME=0;try% Grow the body% for i = 2:k+1%    j = i-1:i;%    set(head,'xdata',x(i),'ydata',y(i),'zdata',z(i))%    set(body,'xdata',x(j),'ydata',y(j),'zdata',z(j))%    drawnow%        % end% Primary loop% h =waitbar(0,'Please wait...');% set(h,'name','仿真进度');tic;m = length(x);for i = 2:m   j = i-1:i;      set(head,'xdata',x(i),'ydata',y(i),'zdata',z(i))   set(body,'xdata',x(j),'ydata',y(j),'zdata',z(j))   %set(tail,'xdata',x(j-k),'ydata',y(j-k),'zdata',z(j-k))   drawnow        t = toc;         str = format_time(TIME + t);        set(findobj('tag','txt_shijian'), 'String', str); %toc 停止计数，设置显示格式，在display 显示结果        set(findobj('tag','txt_feixinggaodu'), 'String', z(i));        set(findobj('tag','txt_vx'), 'String', vx(i));        set(findobj('tag','txt_vy'), 'String', vy(i));        set(findobj('tag','txt_vz'), 'String', vz(i));        set(findobj('tag','txt_hesudu'), 'String', sqrt(vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)));        set(findobj('tag','txt_pianhang'), 'String', ph(i));        set(findobj('tag','txt_fuyang'), 'String', fy(i));        set(findobj('tag','txt_henggun'), 'String', hg(i));        set(findobj('tag','txt_weizhizuobiao'), 'String', ['(' num2str(x(i)) ',' num2str(z(i)) ',' num2str(y(i)) ')']);%         waitbar(i/m,h,[num2str(i*100/m) '%']);end% Clean up the tailfor i = m+1:m+k   j = i-1:i;   set(tail,'xdata',x(j-k),'ydata',y(j-k),'zdata',z(j-k))   drawnowendcatch E    if ~strcmp(E.identifier, 'MATLAB:class:InvalidHandle')        rethrow(E);    endend    % same subfunction as in cometfunction [minx,maxx] = minmax(x)minx = min(x(isfinite(x)));maxx = max(x(isfinite(x)));if minx == maxx  minx = maxx-1;  maxx = maxx+1;end  function str = format_time(t)  hrs = floor(t/3600); min = floor(t/60 - 60*hrs); sec = t - 60*(min + 60*hrs); if hrs < 10   h = sprintf('0%1.0f:', hrs); else   h = sprintf('%1.0f:', hrs); end if min < 10   m = sprintf('0%1.0f:', min); else   m = sprintf('%1.0f:', min); end if sec < 9.9995   s = sprintf('0%1.3f', sec); else   s = sprintf('%1.3f', sec); end str = [h m s];