MATLAB六自由度弹道仿真程序

适用于弹箭飞行控制原理、导弹飞行力学等多门课程 课程设计的需求 在MATLAB平台进行弹道仿真。 适用于弹箭飞行控制原理、导弹飞行力学等多门课程 课程设计的需求 在MATLAB平台进行弹道仿真。

close all;clear;clc
tspan=（0:0.01:55）	;       %t区间、步长
[ty]=ode45（@fun05tspan[260 0 0 7100 225 0 500]）;  %45龙格算法积分参数
%选取特征点，求动力学系数和传递函数
i=find（t==6）    %t在矩阵中位置    0 6  20  34.55 V最大

ma0=[0.60.80.91.01.2];                               
Cx00=[-0.2000-0.2210-0.2315-0.2420-0.2500];
Cxa20=[-0.0005-0.0005-0.0005-0.0005-0.0005];
Cya0=[0.25050.2530.2540.2550.256];
Cydetaz0=[0.040.0470.0430.0450.045];
mza0=[-0.0290-0.0298-0.0300-0.0302-0.0302];
mzdetaz0=[-0.0220-0.0240-0.0232-0.0230-0.0234];
mzwz0=[-0.002-0.0015-0.0015-0.0012-0.0016];
%空军标准气象条件
cc=20.046*（288.34-5.86*10^（-3）*y（i4））^0.5;               %声速
ma=y（i1）/cc;                                             %马赫数
g=9.806*（1.0-2*y（i4）/6371000）;                              %重力加速度
pp=1.225*（1.0-2.0323*10^（-5）*y（i4））^4.83;                  %密度
q= 0.5*pp*y（i1）*y（i1）;                                     %动压
%三次样条插值求各系数
Cx0=interp1（ma0Cx00ma‘spline‘）;
Cxalpha2=interp1（ma0Cxa20ma‘spline‘）;
Cyalpha=interp1（ma0Cya0ma‘spline‘）;
Cydetaz=interp1（ma0Cydetaz0ma‘spline‘）;
Czbeita=-Cyalpha;
mza=interp1（ma0mza0ma‘spline‘）;
mzdetaz=interp1（ma0mzdetaz0ma‘spline‘）;
mzwz=interp1（ma0mzwz0ma‘spline‘）;

L=2.358;     %弹长
Jz=60.85;    %转动惯量
S=0.456;                        %弹翼特征面积
P=3000;                         %推力为3000
Xt=12000;Yt=0;Zt=0;    %目标位置
K=4;
cita=y（i2）;
v=y（i1）;
h=y（i4）;
if i  alpha=5*sin（0.2*pi*t+pi）-y（i2）*57.3  ;                   %攻角  
elseif i  alpha=（y（i5）*g）/（P/57.3+（Cyalpha）*q*S）;                  %攻角
else             %末端攻击段
  r=（（y（i4）-Yt）^2+（y（i3）-Xt）^2）;
  dq=y（i1）*（（y（i3）-Xt）*sin（y（i2））-（y（i4）-Yt）*cos（y（i2）））/r;     %dq/dt  y
 alpha=（y（i5）*y（i1）*K*dq+y（i5）*g*cos（y（i2）））/（P/57.3+（Cyalpha）*q*S）;
end
m=y（i5）;
detaz=-mza/mzdetaz*alpha;

%计算各参数的偏导数
Xv=pp*v*S*（Cx0+Cxalpha2*alpha*alpha）;
Xa=2*Cxalpha2*alpha*q*S;
Mzwz=mzwz*q*S*L*L/v;
Mza=mza*q*S*L;
 Mzdetaz=mzdetaz*q*S*L;
 Yv=v*pp*S*（Cyalpha*alpha+Cydetaz*detaz）;
 Ya=Cyalpha*q*S;
 Ydetaz=Cydetaz*q*S;

%特征点各动力系数，忽略下洗延迟
A11=-Xv/m;A12=0;A13=-g*cos（cita）;A14=-（Xa+P*alpha）/m;A15=0;A16=1/m;
A21=0;A22=Mzwz/Jz;A23=0;A24=Mza/Jz;A25=Mzdetaz/Jz;
A31=Yv/m/v;A32=0;A33=g*sin（cita）/v;A34=（Ya+P）/m/v;A35=Ydetaz/m/v;A36=1/m/v;

%动力系数
a1=-A22
a2=-A24
a3=-A25
a4=A34
a5=A35
%计算弹体传递系数
Km=（-A25*A34+A35*A24）/（A22*A34+A24）     %Km
Tm=1/sqrt（-A24-A22*A34）                 %Tm                         
Em=（-A22+A34）/sqrt（-A24-A22*A34）/2      %ξm 
T1=A25/（A25*A34-A35*A24）                %T1  

Wdz_jieta=tf（Km*[T1 1][Tm^2 2*Tm*Em 1 0]）;  %俯仰角对舵偏的传递函数

Kwz=-0.35;Kp=8;Ki=4;                %设计的参数
%Kwz=-0.08;Kp=20;Ki=4;       
G=tf（1[0.01 1]）;                   %舵机一阶惯性环节
G_Kwz=tf（Kwz*Km*[T1 1][Tm^2 2*Tm*Em 1]）    
figure（‘numbertitle‘‘off‘‘name‘‘内环伯德图‘）
margin（series（G_KwzG））;grid on;                      %内环伯德图
figure（‘numbertitle‘‘off‘‘name‘‘内环闭环传递函数阶跃响应‘）
step（feedbac

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3679  2016-06-21 21:29  913108320128张昊\Bode01.m

     文件       1378  2016-06-21 20:37  913108320128张昊\dandao01.m

     文件       3558  2016-06-21 20:21  913108320128张昊\fun05.m

    I.A.SH.    311296  2016-06-21 21:30  913108320128张昊\Thumbs.db

     文件     391401  2016-06-21 22:08  913108320128张昊\【91310820128张昊】弹箭飞行控制原理.docx

     目录          0  2016-06-21 22:12  913108320128张昊

     文件     676232  2016-06-21 22:12  913108320128张昊\【91310820128张昊】弹箭飞行控制原理.pdf

----------- ---------  ---------- -----  ----

              1387544                    7