%此程序計算全機升力、力矩特性曲線和極曲線
%符號說明
% AR_w：機翼展弦比
% Cd_w：機翼翼型阻力系數（取翼型零升阻力系數）；Cd_fus：機身阻力系數；Cd_tl：尾翼阻力系數
% Cd_bld：螺旋槳阻力系數；Cd_lg：起落架鋼絲阻力系數；Cd_wh：機輪阻力系數
% Cd_w：全機廢阻力系數；Cd_t：全機總阻力系數
%Cl_alpha：升力線斜率（1/deg）；Cl：升力系數；
%Cm_Cl：Cm對Cl的導數；Cm：力矩系數；Cm0：零升力矩系數
% S_w：機翼面積（m^2）；S_fus：機身最大截面積（m^2）；S_lg：起落架鋼絲迎風截面積（m^2）
% S_bld：螺旋槳正面投影面積（m^2）；S_tl：尾翼面積（m^2）；S_wh：機輪最大截面積（m^2）
% Kmax:全機最大升阻比
%alpha：迎角（deg）；alpha_0：零升迎角（deg）

%所需參數值均取自VLM的計算結果
clc;
clear;
Cl_alpha=0.08772;
alpha_0=-3.12166;
alpha=-5:0.1:10;
Cl=Cl_alpha*(alpha-alpha_0);
figure(1)
plot(alpha,Cl);
title('全機升力特性曲線');
xlabel('alpha （deg）');
ylabel('Cl');
Cm_Cl=-0.12106;
Cm0=0.00534;
Cm=Cm0+Cm_Cl*Cl;
figure(2)
plot(Cl,Cm);
title('Cm-Cl曲線');
xlabel('Cl');
ylabel('Cm');
AR_w=8;
S_w=0.575;S_fus=0.071*0.075;S_tl=0.0812+0.0235;S_bld=0.0094;S_lg=0.0029;S_wh=0.0013*3;
Cd_w=0.0091;Cd_fus=0.29;Cd_tl=0.021;Cd_bld=0.7;Cd_lg=1.10;Cd_wh=1.460;
Cd_w=1.1*(Cd_w+Cd_fus*S_fus/S_w+Cd_tl*S_tl/S_w+Cd_bld*S_bld/S_w+Cd_lg*S_lg/S_w+Cd_wh*S_wh/S_w);
for i=1:151;
Cd_t(i)=Cd_w+Cl(i)^2/(pi*AR_w);
error=(Cd_w-0.5*Cd_t(i))/Cd_w;
if abs(error)<=0.005;
    fprintf('有利升力系數為 %.3f \n',Cl(i));
    Kmax=max(Cl(i)/Cd_t(i));
    fprintf('全機最大升阻比為 %.3f \n',Kmax);
end
end
figure(3)
plot(Cd_t,Cl);
title('全機極曲線');
xlabel('Cd');
ylabel('Cl');