%龍格庫塔右端函數，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=y(4)-y(2);
ABS_alpha=abs(alpha);
%插值部分開始
P=interp1(APT,AP,y(8))*g;
CX=twointerp(Aa,AMa,ACX,ABS_alpha,Ma);
CY=twointerp(Aa,AMa,ACX,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(8)<=B(1) 
    MC=2.362;
else  if y(8)>B(2) 
         MC=0;
         P=0;
      else 
         MC=0.21059;
      end
  end
  
if y(8)<B(2)
XG=interp1(AXGT,AXG,y(8));  %插值重心
end

if y(8)<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(8));
dy(1)=(P*cos(alpha)-XF-9.81*y(7)*sin(y(2)))/y(7);
dy(2)=(P*sin(alpha)+YF-9.81*y(7)*cos(y(2)))/(y(7)*y(1));
dy(3)=(MZAF*alpha+MZWZ*wzt)*q*s*L/JZ;
dy(4)=y(3);
dy(5)=y(1)*cos(y(2));
dy(6)=y(1)*sin(y(2));
dy(7)=-MC;
dy(8)=1;