%符號說明
%b1,b2,b3,b4為橫航向小擾動運動方程的特征方程的系數
%m：飛機質量（kg）；V：飛行速度（m/s）；a：迎角(deg)；
%CY_b：側力對側滑角的導數
%Cl_b：滾轉力矩系數對側滑角的導數；Cl_R：滾轉力矩系數對偏航角速度的導數
%Cl_P：滾轉力矩系數對滾轉角速度的導數；Cn_b：偏航力矩系數對側滑角的導數
%Cn_P：偏航力矩系數對滾轉角速度的導數；Cn_R：偏航力矩系數對偏航角速度的導數
%CY_b1、Cl_b1、Cl_P1、Cl_R1、Cn_b1、Cn_P1、Cn_R1只是中間符號，其含義如程序中表達式所示
%Ix：飛機相對于x軸的慣性矩；Iz：飛機相對于z軸的慣性矩；Ixz：飛機相對于xz軸的慣性積
%機體坐標系：x軸順機身向后為正，y指向右翼為正，z軸與xy軸成右手系，原點為全機質心

%參照各參數的含義，輸入對應值運行即可

m=2;V=15;Ix=;Iz=;Ixz=;g=9.81;a=3.39;
CY_b=-0.10948;Cl_b=0.034154;Cl_P=-0.54149;Cl_R=0.044063;
Cn_b=-0.037591;Cn_P=-0.050911;Cn_R=-0.034918;
CY_b1=CY_b/(m*V);
Cl_P=Cl_P/Ix;
Cl_R=Cl_R/Ix;
Cl_b=Cl_b/Ix;
Cn_R=Cn_R/Iz;
Cn_P=Cn_P/Iz;
Cn_b=Cn_b/Iz;
I=1-Ixz^2/(Ix*Iz);
Cl_P1=(Cl_P+Cn_P*Ixz/Ix)/I;
Cl_R1=(Cl_R+Cn_R*Ixz/Ix)/I;
Cl_b1=(Cl_b+Cn_b*Ixz/Ix)/I;
Cn_P1=(Cn_P+Cl_P*Ixz/Iz)/I;
Cn_R1=(Cn_R+Cl_R*Ixz/Iz)/I;
Cn_b1=(Cn_b+Cl_b*Ixz/Iz)/I;
b1=-(CY_b1+Cl_P1+Cn_R1)
b2=CY_b1*(Cn_R1+Cl_P1)+(Cn_R1*Cl_P1-Cn_P1*Cl_R1)-Cn_b1-(a*pi/180)*Cl_b1
b3=Cn_b1*Cl_P1-Cl_b1*(g/V+Cn_P1)+CY_b1*(Cn_P1*Cl_R1-Cl_P1*Cn_R1)+(a*pi/180)*(Cl_b1*Cn_R1-Cn_b1*Cl_R1)
b4=(g/V)*(Cl_b1*Cn_R1-Cn_b1*Cl_R1)
R=b1*b2*b3-b1^2*b4-b3^2
if b4>0&&R>0&&Cn_b<0&&Cl_b>0
    disp('橫航向靜安定且橫航向動安定');
else if b4<=0&&R>0&&Cn_b<0&&Cl_b>0
            disp('橫航向靜安定且振蕩安定（航向動安定），但螺旋不安定（橫向動不安定）');
else if b4>0&&R<=0&&Cn_b<0&&Cl_b>0
            disp('橫航向靜安定且螺旋安定（橫向動安定），但振蕩不安定（航向動不安定）');
else if b4>0&&R>0&&Cn_b>=0&&Cl_b>0
            disp('橫航向動安定且橫向靜安定，但航向靜不安定');
else if b4>0&&R<0&&Cn_b<0&&Cl_b<=0
            disp('橫航向動安定且航向靜安定，但橫向靜不安定');
else if b4<=0&&R<=0&&Cn_b<0&&Cl_b>0
            disp('橫航向靜安定，但橫航向動不安定');
else if b4<=0&&R>0&&Cn_b>=0&&Cl_b>0
            disp('橫向靜安定，但航向靜不安定；振蕩安定（航向動安定），但螺旋不安定（橫向動不安定）');
else if b4<=0&&R>0&&Cn_b<0&&Cl_b<=0
            disp('航向靜安定，但橫向靜不安定；振蕩安定（航向動安定），但螺旋不安定（橫向動不安定）');
else if b4>0&&R<=0&&Cn_b>=0&&Cl_b>0
            disp('橫向靜安定，但航向靜不安定；螺旋安定（橫向動安定），但振蕩不安定（航向動不安定）');
else if b4>0&&R<=0&&Cn_b<0&&Cl_b<=0
            disp('航向靜安定，但橫向靜不安定；螺旋安定（橫向動安定），但振蕩不安定（航向動不安定）');
else if b4>0&&R>0&&Cn_b>=0&&Cl_b<=0
            disp('橫航向動安定，但橫航向靜不安定');
else if b4<=0&&R<=0&&Cn_b>=0&&Cl_b>0
            disp('橫向靜安定，但航向靜不安定，且橫航向動不安定');
else if b4<=0&&R<=0&&Cn_b<0&&Cl_b<=0
            disp('航向靜安定，但橫向靜不安定；且橫航向動不安定');
else if b4<=0&&R>0&&Cn_b>=0&&Cl_b<=0
            disp('橫航向靜不安定；振蕩安定（航向動安定），但螺旋不安定（橫向動不安定）');
else if b4>0&&R<=0&&Cn_b>=0&&Cl_b<=0
            disp('橫航向靜不安定；螺旋安定（橫向動安定），但振蕩不安定（航向動不安定）');
else if b4<=0&&R<=0&&Cn_b>=0&&Cl_b<=0
            disp('橫航向靜不安定且橫航向動不安定');
        end
    end
end
end
end
end
end
end
end
end
end
end
end
end
end
end