利用MATLAB进行弹道仿真源代码

%---------------------------------------------% % % % 工作室提供代做matlab仿真 % % % % 详情请访问：http://cn.mikecrm.com/5k6v1DP % % % %---------------------------------------------% clc clear global a b c m_0 d L p g y_2 y_og t_1 a_0 i=1; n=1; prompt = {' a :------(o=at+a_0)-- ';' a_0 :--';' b:------(F=bt)-- ';' c :------(m=m_0-ct)-- ';' d :------(M=dm_0)-- ';' m_0 :--';' p :-----(p=m/l)-- ';' L :------';' g:--';' pause(?) :';' hold off ( y / n ):';' y_2 :---'}; dlg_title = 'missile......'; num_lines= 1; def = {'-.015','.001','30','.2','6','60','1','1100','9.81','.05','y','5000'}; inp = inputdlg(prompt,dlg_title,num_lines,def); while(strcmp(inp,'')==0) hh = strcmp(inp(11,1),'y'); inp = str2mat( inp ); a = str2num(inp(1,:)); a_0 = str2num(inp(2,:)); b = str2num(inp(3,:)); c = str2num(inp(4,:)); d = str2num(inp(5,:)); m_0 = str2num(inp(6,:)); p = str2num(inp(7,:)); L = str2num(inp(8,:)); g = str2num(inp(9,:)); pp = str2num(inp(10,:)); y_2 = str2num(inp(12,:)); if(L < (m_0)/p) msgbox('L < l') break end % #0# x_0=0 ; y_0 = L/2-((m_0)*(L-(m_0)/p))/(2*((1+d)*m_0)); fi_0=90; t_0=0; m_0; % #1# y_og1=((m_0-c*t_1)*(L-(m_0-c*t_1)/p))/(2*((1+d)*m_0-c*t_1)); x_1=0; y_1=L/2-y_og1; fi_1=90; t_1=m_0*(1+d)/(b/g+c); m_1=m_0-c*t_1; v_22=-13/10*(1175/2+13/40*t_1)/(600-13/5*t_1)+13/40*(50-13/10*t_1)/(600-13/5*t_1)+13/5*(50-13/10*t_1)*(1175/2+13/40*t_1)/(600-13/5*t_1)^2; % #2# tt=1; t_2 = t_1 + tt; [T_2,Y_2] = ode45(@state2,[t_1 t_2],[L/2-y_og1 v_22]); while(Y_2(length(Y_2(:,1)),1) < (y_2 + y_1)) tt = tt + 1;t_2 = t_1 + tt; [T_2,Y_2] = ode45(@state2,[t_1 t_2],[L/2-y_og1 v_22]); end x_2 = 0; fi_2 = 90; m_2 = m_0-c*t_2; y_2 = Y_2(length(Y_2(:,1)),1); v_y2 = Y_2(length(Y_2(:,2)),2); if(m_2 < 0) msgbox('fule is low (1)') break end % #3# t_3 = m_2/c;mm=0; tt_3=1; [T_3,Y_3]=ode45(@state3,[t_2 t_2+tt_3],[0 y_2 90 0 v_y2 0]); while((Y_3(length(Y_3(:,5)),1) > 0)) tt_3 = tt_3 + 1;t_2 = t_1 + tt_3; [T_3,Y_3]=ode45(@state3,[t_2 t_2+tt_3],[0 y_2 90 0 v_y2 0]); mm=mm+1; end if(t_3 < tt_3) msgbox('fule is low (2)') break end x_3 = Y_3(length(Y_3(:,1)),1); y_3 = Y_3(length(Y_3(:,2)),2); fi_3 = Y_3(length(Y_3(:,3)),3); v_x3 = Y_3(length(Y_3(:,4)),4); v_y3 = Y_3(length(Y_3(:,5)),5); w_fi3 = Y_3(length(Y_3(:,6)),6); % #4# t_4 = m_0/c; [T_4,Y_4]=ode45(@state4,[t_2+tt_3 t_4],[x_3 y_3 fi_3 v_x3 v_y3 w_fi3]); x_4 = Y_4(length(Y_4(:,1)),1); y_4 = Y_4(length(Y_4(:,2)),2); fi_4 = Y_4(length(Y_4(:,3)),3); v_x4 = Y_4(length(Y_4(:,4)),4); v_y4 = Y_4(length(Y_4(:,5)),5); w_fi4 = Y_4(length(Y_4(:,6)),6); % #5# tt=1; [T_5,Y_5]=ode45(@state5,[t_4 t_4+tt],[x_4 y_4 fi_4 v_x4 v_y4 w_fi4]); while(Y_5(length(Y_5(:,2)),2) > 0) [T_5,Y_5]=ode45(@state5,[t_4 t_4+tt],[x_4 y_4 fi_4 v_x4 v_y4 w_fi4]); tt = tt +1; end plot(zeros(length(Y_2),1),Y_2(:,1) ,Y_3(:,1),Y_3(:,2) ,'-' ,Y_4(:,1),Y_4(:,2),'--' ,Y_5(:,1),Y_5(:,2),'b--' ),xlabel('x'),ylabel('y'),grid pause T_1=[0:1:t_1]'; y_og2=L/2-((m_0-c*T_1).*(L-(m_0-c*T_1)/p))./(2*((1+d)*m_0-c*T_1)) ; v_2=-13./10*(1175./2+13./40.*T_1)./(600-13./5.*T_1)+13./40.*(50-13./10.*T_1)./(600-13./5.*T_1)+13./5.*(50-13./10.*T_1).*(1175./2+13./40.*T_1)./(600-13./5.*T_1).^2; GG=zeros(length(T_1),1); subplot(2,3,1),plot(T_1,GG,T_2,zeros(length(Y_2),1),T_3,Y_3(:,1),'-g',T_4,Y_4(:,1),'--'),ylabel('x'),xlabel('t'),grid subplot(2,3,2),plot(T_1, y_og2 ,T_2,Y_2(:,1),T_3,Y_3(:,2) ,'-g',T_4,Y_4(:,2),'--'),ylabel('y'),xlabel('t'),grid subplot(2,3,3),plot(T_1,GG,T_2,90*ones(length(Y_2),1),T_3,Y_3(:,3),'-g',T_4,Y_4(:,3),'--'),ylabel('fi'),xlabel('t'),grid subplot(2,3,4),plot(T_1,GG,T_2,zeros(length(Y_2),1),T_3,Y_3(:,4),'-g',T_4,Y_4(:,4),'--'),ylabel('dx'),xlabel('t'),grid subplot(2,3,5),plot(T_1,v_2,T_2,Y_2(:,2),T_3,Y_3(:,5),'-g',T_4,Y_4(:,5),'--'),ylabel('dy'),xlabel('t'),grid subplot(2,3,6),plot(T_1,GG,T_2,zeros(length(Y_2),1),T_3,Y_3(:,6),'-g',T_4,Y_4(:,6),'--'),ylabel('dfi'),xlabel('t'),grid pause cc=Y_3(1,2); i=1; while(abs((Y_3(i,2)-cc)) < 2) i=i+1; end X=[GG;zeros(length(Y_2(:,1)),1);Y_3(i:length(Y_3(:,2)),1);Y_4(:,1);Y_5(:,1)]; Y=[y_og2;Y_2(:,1);Y_3(i:length(Y_3(:,2)),2);Y_4(:,2);Y_5(:,2)]; FI=[90*ones(length(Y_2(:,1))+length(T_1),1) ; Y_3(i:length(Y_3(:,2)),3) ; Y_4(:,3);Y_5(:,3)]; V_x=[GG;zeros(length(Y_2(:,1)),1);Y_3(i:length(Y_3(:,2)),4);Y_4(:,4);Y_5(:,4)]; V_y=[v_2;Y_2(:,2);Y_3(i:length(Y_3(:,2)),5);Y_4(:,5);Y_5(:,5)]; T=[T_1;T_2;T_3(i:length(T_3),1);T_4;T_5]; X_max=X(1,1); Y_max=Y(1,1); X_min=X(1,1); Y_min=Y(1,1); for i=1:1:length(X) if(X(i,1) > X_max) X_max=X(i,1); end if(X(i,1) < X_min) X_min=X(i,1); end if(Y(i,1) > Y_max) Y_max=Y(i,1); end if(Y(i,1) < Y_min) Y_min=Y(i,1); end end %-------------------------1--------------------------- subplot(1,1,1) for i=1:1:length(Y) if(T(i,1) <= (t_2 + t_3)) y_og = ((m_0-c*T(i,1))*(L-(m_0-c*T(i,1))/p))/(2*((1+d)*m_0-c*T(i,1))) ; L_m = (m_0-c*T(i,1))/p; aa = 7/200; L_w = (1-c*T(i,1)/(m_0*(1+d)))*L/5; L_F = L*b*T(i,1)/(5*m_0*(1+d)*g); else y_og=0; L_m=0;aa=0; L_w=L/(5*(1+1/d)); L_F=0; end if((T(i,1) <= (t_2))|(T(i,1) >= (t_2+t_3))) cc=0; else cc=a*(T(i,1)-t_2); end plot([X(i,1)-(L/2-y_og)*cos( FI(i,1)*pi/180 ),X(i,1)+(L/2+y_og)*cos( FI(i,1)*pi/180 )],[Y(i,1)-(L/2-y_og)*sin(FI(i,1)*pi/180),Y(i,1)+(L/2+y_og)*sin(FI(i,1)*pi/180)]); xlim([X_min - L , X_max + L]) ylim([Y_min - L , Y_max + L]) hold on XX_0 = X(i,1) - (L/2-y_og)*cos(FI(i,1)*pi/180); YY_0 = Y(i,1) - (L/2-y_og)*sin(FI(i,1)*pi/180); XX_1 = XX_0 - aa*L*sin( FI(i,1)*pi/180 ); XX_2 = XX_1 + L_m*cos( FI(i,1)*pi/180 ); XX_3 = XX_2 + aa*2*L*sin( FI(i,1)*pi/180 ); XX_4 = XX_3 - L_m*cos( FI(i,1)*pi/180 ); YY_1 = YY_0 + aa*L*cos( FI(i,1)*pi/180 ); YY_2 = YY_1 + L_m*sin( FI(i,1)*pi/180 ); YY_3 = YY_2 - aa*2*L*cos( FI(i,1)*pi/180 ); YY_4 = YY_3 - L_m*sin( FI(i,1)*pi/180 ); plot(X(i,1),Y(i,1),'+g') plot([XX_1,XX_2,XX_3,XX_4,XX_1],[YY_1,YY_2,YY_3,YY_4,YY_1],'r'); %fill([XX_1,XX_2,XX_3,XX_4,XX_1],[YY_1,YY_2,YY_3,YY_4,YY_1],'r'); bb=2/15; %----------w plot([X(i,1),X(i,1) , X(i,1)-bb*L_w/((2)^.5) , X(i,1)-bb*L_w/((2)^.5)+(2)^.5*bb*L_w , X(i,1)],... [Y(i,1),Y(i,1)-L_w, Y(i,1)-L_w+bb*L_w/((2)^.5) , Y(i,1)-L_w+bb*L_w/((2)^.5) , Y(i,1)-L_w],'k') %----------F plot([XX_0 , XX_0-bb*L_F*cos((cc+FI(i,1)-45)*pi/180) , XX_0+bb*L_F*sin((cc+FI(i,1)-45)*pi/180) , XX_0 , XX_0-L_F*cos((cc+FI(i,1))*pi/180)],... [YY_0 , YY_0-bb*L_F*sin((cc+FI(i,1)-45)*pi/180) , YY_0-bb*L_F*cos((cc+FI(i,1)-45)*pi/180) , YY_0 , YY_0-L_F*sin((cc+FI(i,1))*pi/180)],'k') %----------f_x L_vx = L/5*(sign(V_x(i,1)).*0.1.*V_x(i,1).^2)/(m_0*g*(1+d)); plot([X(i,1) + y_og*cos( FI(i,1)*pi/180 ) , X(i,1) + y_og*cos( FI(i,1)*pi/180 ) + L_vx*bb*((2)^(.5))/2 , X(i,1) + y_og*cos( FI(i,1)*pi/180 ) + L_vx*bb*((2)^(.5))/2 , X(i,1) + y_og*cos( FI(i,1)*pi/180 ) , X(i,1) + y_og*cos( FI(i,1)*pi/180 ) + L_vx],... [Y(i,1) + y_og*sin( FI(i,1)*pi/180 ) , Y(i,1) + y_og*sin( FI(i,