%此程序計算巡航狀態縱向動穩定性
%算法取自熊海泉、劉旭、鄭本武編著《飛機飛行動力學》178～180頁

m=2.4;S_w=0.575;ba=0.268;Iz=50;yp=0.01;%????Iz=50;yp=0.01;
g=9.81;den=1.225;a=340;
alpha=0.26;V0=15;P=4.5;M=0.045;
Cl=0.2966;Cd=0.050338;Cl_alpha=5.0189;Cd_alpha=0.12227;
Cl_M=4*m*g/(den*S_w*a^2*M^3);
Cd_M=4*P/(den*S_w*a^2*M^3);
mz_M=4*P*yp/(den*S_w*a^2*M^3);
P_V=den*S_w*Cd*V0;
mz_alpha=-0.88253;mz_Q=-8.98;mz_alpha_=0.5*mz_Q;
q=0.5*den*S_w*V0^2;

Xc_alpha=q*S_w*Cd_alpha/m
Xc_V=P_V/m-q*S_w*(2*Cd+M*Cd_M)/(m*V0)
Y_alpha=q*S_w*Cl_alpha/(m*V0)
Y_V=q*S_w*(2*Cl+M*Cl_M)/(m*V0^2)
Mz_alpha=q*S_w*ba*mz_alpha/Iz
Mzc_V=(q*S_w*ba*mz_M/alpha-P_V*yp)/Iz
Mz_Q=q*S_w*ba*mz_Q/Iz
Mz_alpha_=q*S_w*ba*mz_alpha_/Iz

A1=Y_alpha-(Xc_V+Mz_Q+Mz_alpha_)
A2=Y_V*(Xc_alpha+g)-Xc_V*(Y_alpha-Mz_alpha_-Mz_Q)-Mz_alpha-Mz_Q*Y_alpha
A3=-Xc_alpha*(Y_V*Mz_Q+Mzc_V)-g*Y_V*(Mz_Q+Mz_alpha_)+Xc_V*(Y_alpha*Mz_Q+Mz_alpha)
A4=g*(Y_alpha*Mzc_V-Y_V*Mz_alpha)
R=A1*A2*A3-A1^2*A4-A3^2

if  A1>0&&A2>0&&A3>0&&A4>0&&R>0
     fprintf('\n縱向動穩定\n');
 else  fprintf('\n縱向動不穩定\n');
 end