grkt2f.asv

垂直发射弹道仿真程序(采用四元数法).matlab编写。

%龙格库塔右端函数，dy:导数值；y:积分值

function dy=grkt2f(y)

global AJZ;
global AJZT;
global AMa;
global Aa;
global ACX;
global ACY;
global AP;
global APT;
global AXG;
global AXGT;
global AMZWZ1;
global AMZWZ2;
global AMZAF;

RHO=1.225;
SONIC=343.13;

s=0.0227;             %特征面积
q=0.5*RHO*y(1)*y(1);          %动压头
Ma=y(1)/SONIC;          %计算马赫数
L=1.8;                  %特征长度
B=[2.1126 44.0832];      %两级发动机结束工作时间      
g=9.81;


alpha=asin(2*(y(9)*y(10)+y(8)*y(11)))-y(2);
ABS_alpha=abs(alpha);
%插值部分开始
P=interp1(APT,AP,y(7))*g;
CX=twointerp(Aa,AMa,ACX,ABS_alpha,Ma);

CY=twointerp(Aa,AMa,ACY,ABS_alpha,Ma);

if alpha<0
    CY=-CY;
end
XF=CX*q*s;
YF=CY*q*s;
wzt=y(3)*L/y(1);

if y(7)<=B(1) 
    MC=2.362;
else  if y(8)>B(2) 
         MC=0;
         P=0;
      else 
         MC=0.21059;
      end
  end
  
if y(7)<B(2)
XG=interp1(AXGT,AXG,y(7));  %插值重心
end

if y(7)<B(2)
  aa=twointerp(Aa,AMa,AMZWZ1,ABS_alpha,Ma);
  bb=twointerp(Aa,AMa,AMZWZ2,ABS_alpha,Ma);
  ab=[aa,bb];
  a=[.9831,.8896];
  MZWZ=interp1(a,ab,XG);
else
    MZWZ=twointerp(Aa,AMa,AMZWZ2,ABS_alpha,Ma);
end

MZAF=twointerp(Aa,AMa,AMZAF,ABS_alpha,Ma); 
if ABS_alpha>0
    MZAF=MZAF*57.3+CY*(XG-0.9381)/(L*alpha); %由于重心的改变，对MZAF进行修正
end
%插值部分结束

JZ=interp1(AJZT,AJZ,y(7));
dy(1)=(P*cos(alpha)-XF-9.81*y(6)*sin(y(2)))/y(6);
dy(2)=(P*sin(alpha)+YF-9.81*y(6)*cos(y(2)))/(y(6)*y(1));
dy(3)=(MZAF*alpha+MZWZ*wzt)*q*s*L/JZ;
dy(4)=y(1)*cos(y(2));
dy(5)=y(1)*sin(y(2));
dy(6)=-MC;
dy(7)=1;
dy(8)=-0.5*y(11)*y(3);
dy(9)=0.5*y(10)*y(3);
dy(10)=-0.5*y(9)*y(3);
dy(11)=0.5*y(8)*y(3);