%本程序計算全機各部分參數
%符號說明
%A_ht：平尾尾容量；Angle_to：起飛迎角（deg）；A_vt：垂尾尾容量
%AR_ht：平尾展弦比；AR_vt：垂尾展弦比；AR_w：機翼展弦比
%Cl_max：全機最大升力系數（失速時）；Cl_min：全機最小升力系數（最大速度時）
%Cl_cru：全機巡航升力系數（即設計升力系數）
%H_fus：機身中段高度（m）
%K_to：起飛升阻比；K_V_max：最大平飛速度時的升阻比
%L_fus：機身長度（m）；L_ht：平尾尾力臂（m）；L_vt：垂尾尾力臂（m）
%P_to：起飛時螺旋槳輸出功率（w）；P_mot_to：起飛時電機輸出功率（w）
%P_max：螺旋槳最大輸出功率（w）；P_mot_max：電機最大輸出功率（w）
%PWR_to：起飛功重比（W/kg）；PWR_max：最大功率重量比（W/kg）
%Re_max：全機最大雷諾數；Re_cru：全機巡航雷諾數；Re_min：全機最小雷諾數
%S_w：機翼面積(m^2)；S_ail：副翼面積（m^2）；S_ht：平尾面積(m^2)
%S_ele：升降舵面積(m^2)；S_vt：垂尾面積(m^2)；S_rud：方向舵面積(m^2)
%TR_w：機翼根梢比；TR_ht：平尾根梢比；TR_vt：垂尾根梢比
%T_to：起飛時螺旋槳拉力（kg）；T_max：最大速度時螺旋槳拉力（kg）
%V_max：最大飛行速度（m/s）；V_min：最小飛行速度（m/s）；V_cru：巡航飛行速度（m/s）
%WL：翼載荷(kg/m^2)；Wto：起飛重量(g)；W_fus：機身寬度（m）
%W_w：機翼結構重量（g）；W_fus：機身結構重量（g）；W_ht：平尾結構重量（g）
%W_vt：垂尾結構重量（g）；W_lg：起落架結構重量（g）；Wstr：結構重量（g）
%Wfix：固定設備重量（g）；Wele：電子設備（舵機、接收機、連桿等）重量（g）
%Weff：有效載荷（g）；Wthr：推進系統重量（g）；Wcel：電池重量（g）
%Wmot：電機、減速器、調速器、安裝架重量和（g）；Wbld：螺旋槳和槳罩重量（g）
%X0_vt：垂尾前緣后掠角（deg）；X4_vt：垂尾1/4弦線后掠角（deg）
%c_w：機翼平均氣動弦長(m)；c_ail：副翼弦長（m）；c_ht：平尾平均氣動弦長(m)
%c_ele：升降舵弦長(m)；ca_vt：垂尾平均氣動弦長(m)；c_rud：方向舵弦長(m)
%den:大氣密度（kg/m^3）;dih_i:機翼中央翼段上反角（deg）;dih_o:機翼外翼段上反角(deg)
%eta_bld：螺旋槳效率；eta_gear：減速器效率；eta_mot：電機效率
%g：重力加速度（m/s^2）
%l_w_i_eff：機翼中央翼段有效展長（m）；l_w_i：機翼中央翼段實際長度（m）
%l_w_o_eff：機翼外翼段有效展長（m）；；l_w_o：機翼外翼段實際長度（m）
%l_w：機翼展長(m)；l_ail：副翼展長（m）；l_ht：平尾展長(m)
%l_ele：升降舵展長(m)；l_vt：垂尾展長(m)；l_rud：方向舵展長(m)
%mu：空氣粘性系數（kg/(m*s)）
%r1：副翼面積與機翼面積之比；r2：副翼弦長與機翼平均氣動弦長之比
%r3：升降舵面積與平尾面積之比 
%tc：機翼相對厚度(%)
%全機總體參數及機翼
clc;
clear;
disp('*************全機主要參數及機翼參數計算**************');
disp(' ');
disp('************本程序只支持無扭轉的矩形機翼*************');
disp(' ');
Wto=input('請輸入初設全機重量值（g） Wto=');
WL=input('請輸入初設翼載荷值（kg/m2） WL=');
AR_w=input('請輸入機翼初設展弦比 AR_w=');
l_w_i_eff=input('請輸入機翼中央翼段初設有效展長（m） l_w_i_eff=');
dih_i=input('請輸入機翼中央翼段初設上反角（deg） dih_i=');
dih_o=input('請輸入機翼外翼段初設上反角（deg） dih_o=');
V_min=input('請輸入最小平飛速度（m/s） V_min=');
V_cru=input('請輸入巡航速度（m/s） V_cru=');
V_max=input('請輸入最大平飛速度（m/s） V_max=');
r1=input('請輸入副翼面積與機翼面積之比 r1=');
r2=input('請輸入副翼弦長與機翼平均氣動弦長之比 r2=');
disp(' ');
disp('****************************************************');
%Wto=2400;WL=4.17;AR_w=8;dih_i=0;dih_o=3;
den=1.225;mu=1.7894*10^(-5);g=9.81;
%V_min=8;V_cru=15;V_max=20;
S_w=Wto/1000/WL;
l_w=(S_w*AR_w)^0.5;
%l_w_i=1.00;
l_w_i=l_w_i_eff/cos(dih_i*pi/180);
l_w_o_eff=(l_w-l_w_i_eff)/2;
l_w_o=l_w_o_eff/cos(dih_o*pi/180);
c_w=l_w/AR_w;
Cl_max=2*Wto*g/(den*V_min^2*S_w)/1000;
Cl_cru=2*Wto*g/(den*V_cru^2*S_w)/1000;
Cl_min=2*Wto*g/(den*V_max^2*S_w)/1000;
Re_max=den*V_max*c_w/mu;
Re_cru=den*V_cru*c_w/mu;
Re_min=den*V_min*c_w/mu;
fprintf('\n翼載荷為 %.3f kg/m^2\n\n',WL);
fprintf('最大升力系數為 %.3f \n巡航升力系數為 %.3f\n最小升力系數為 %.3f\n\n',Cl_max,Cl_cru,Cl_min);
fprintf('最大雷諾數為 %d \n巡航雷諾數為 %d\n最小雷諾數為 %d\n\n',Re_max,Re_cru,Re_min);
fprintf('機翼面積 %.4f m^2\n機翼展長 %.3f m\n機翼平均氣動弦長 %.3f m\n\n',S_w,l_w,c_w);
fprintf('中央翼段有效展長 %.3f m\n中央翼段上反角 %.3f deg\n中央翼段實際長度 %.3f m\n\n',l_w_i_eff,dih_i,l_w_i);
fprintf('外翼段有效展長 %.3f m\n外翼段上反角 %.3f deg\n外翼段實際長度 %.3f m\n\n',l_w_o_eff,dih_o,l_w_o);
%r1=0.08;r2=0.19;
S_ail=r1*S_w/2;
c_ail=r2*c_w;
l_ail=S_ail/c_ail;
fprintf('副翼面積 %.4f m^2\n副翼展長 %.3f m\n副翼弦長 %.3f m\n\n',S_ail,l_ail,c_ail);
%動力系統
disp('********************動力裝置計算********************');
disp(' ');
eta_bld=input('請輸入螺旋槳效率 eta_bld＝');
eta_gear=input('請輸入減速器效率 eta_gear＝');
eta_mot=input('請輸入電機效率 eta_mot＝');
K_to=input('請輸入起飛離地時升阻比 K_to＝');
Angle_to=input('請輸入起飛離地時迎角（deg） Angle_to＝');
K_V_max=input('請輸入最大速度時升阻比 K_V_max＝');
disp(' ');
disp('****************************************************');
%eta_bld=0.6;eta_gear=0.98;eta_mot=0.7;K_to=9;Angle_to=10;K_V_max=3.8;
Cl_ld=0.85*Cl_max;       %Cl_ld與Cl_max的比值因飛機而異，約在0.8～0.9之間
V_ld=(2*Wto*g/(1000*den*S_w*Cl_ld))^0.5;
T_to=Wto*(cos(Angle_to*pi/180)+K_to*sin(Angle_to*pi/180))/(1000*K_to);
P_to=T_to*g*V_ld;
eta=eta_bld*eta_gear*eta_mot;
P_mot_to=P_to/eta;
PWR_to=P_mot_to/(Wto/1000);
fprintf('\n起飛離地時螺旋槳拉力為 %.3f kg\n起飛離地時螺旋槳輸出功率為 %.3f W\n\n',T_to,P_to);
fprintf('起飛離地時電機輸出功率為 %.3f W\n起飛離地時功率重量比為 %.3f W/kg\n\n',P_mot_to,PWR_to);
T_max=Wto/1000/K_V_max;
P_max=T_max*g*V_max;
P_mot_max=P_max/eta;
PWR_max=P_mot_max/(Wto/1000);
fprintf('最大速度時螺旋槳拉力為 %.3f kg\n螺旋槳最大輸出功率為 %.3f W\n\n',T_max,P_max);
fprintf('電機最大輸出功率為 %.3f W\n最大功率重量比為 %.3f W/kg\n\n',P_mot_max,PWR_max);
%平尾
disp('********************平尾參數計算**********************');
disp(' ');
disp('*****************此程序只支持矩形平尾*****************');
disp(' ');
A_ht=input('請輸入平尾尾容量 A_ht＝');
L_ht=input('請輸入平尾尾力臂（m） L_ht＝');
AR_ht=input('請輸入平尾展弦比 AR_ht＝');
r3=input('請輸入升降舵面積與平尾面積之比 r3=');
disp(' ');
disp('*****************************************************');
%A_ht=0.39;L_ht=0.741;AR_ht=3.4;r3=0.2;
S_ht=S_w*c_w*A_ht/L_ht;
l_ht=(S_ht*AR_ht)^0.5;
c_ht=l_ht/AR_ht;
S_ele=r3*S_ht;
c_ele=r3*c_ht;
l_ele=S_ele/c_ele;
fprintf('\n平尾面積 %.4f m^2\n平尾展長 %.3f m\n平尾平均氣動弦長 %.3f m\n\n',S_ht,l_ht,c_ht);
fprintf('升降舵面積 %.4f m^2\n升降舵展長 %.3f m\n升降舵弦長 %.3f m\n\n',S_ele,l_ele,c_ele);
%垂尾
disp('********************垂尾參數計算*********************');
disp(' ');
A_vt=input('請輸入垂尾尾容量 A_vt＝');
L_vt=input('請輸入垂尾尾力臂（m） L_vt＝');
AR_vt=input('請輸入垂尾展弦比 AR_vt＝');
TR_vt=input('請輸入垂尾根梢比 TR_vt＝');
X4_vt=input('請輸入垂尾1/4弦線后掠角（deg） X4_vt＝');
r4=input('請輸入方向舵面積與垂尾面積之比 r4=');
disp(' ');
disp('*****************************************************');
%A_vt=0.014;L_vt=0.735;AR_vt=0.92;TR_vt=2.55;X4_vt=30;r4=0.2;
S_vt=S_w*l_w*A_vt/L_vt;
l_vt=(S_vt*AR_vt)^0.5;
c0_vt=2*S_vt/(l_vt*(1+1/TR_vt));
X0_vt=atan(tan(X4_vt*pi/180)+(1-1/TR_vt)/(AR_vt*(1+1/TR_vt)))*180/pi;
ca_vt=2*c0_vt*(1+TR_vt+TR_vt^2)/(3*(1+TR_vt));
S_rud=r4*S_vt; 
l_rud=l_vt;
c_rud=S_rud/l_rud;
fprintf('\n垂尾面積 %.4f m^2\n垂尾展長 %.3f m\n',S_vt,l_vt);
fprintf('垂尾根弦長 %.3f m\n垂尾梢弦長 %.3f m\n垂尾平均氣動弦長 %.3f m\n',c0_vt,c0_vt/TR_vt,ca_vt);
fprintf('垂尾前緣后掠角 %.3f deg\n1/4弦線后掠角 %.3f deg\n\n',X0_vt,X4_vt);
fprintf('方向舵面積 %.4f m^2\n方向舵展長 %.3f m\n方向舵弦長 %.3f m\n\n',S_rud,l_rud,c_rud);
%全機及各部分重量
disp('****************全機及各部分重量計算******************');
disp(' ');
W_fus=input('請輸入中段機身寬度（m） W_fus＝');
H_fus=input('請輸入中段機身高度（m） H_fus＝');
L_fus=input('請輸入機身長度（m） L_fus＝');
tc=input('請輸入機翼翼型最大相對厚度(%) tc＝');
Wele=input('請輸入電子設備（舵機、接收機等）重量（g） Wele=');
Weff=input('請輸入有效載荷重量（g） Weff=');
Wmot=input('請輸入電機、調速器、減速器重量和（g） Wmot=');
Wcel=input('請輸入電池組重量（g） Wcel=');
Wbld=input('請輸入螺旋槳重量（g） Wbld=');
disp(' ');
disp('*****************************************************');
%H_fus=0.075;L_fus=1.124;W_fus=0.071;tc=9.154;
%機翼
W_w=453.6*0.14676*(S_w*3.2808^2)^0.4852*AR_w^0.7082*tc^(-0.2210);
%機身
W_fus=453.6*0.07092*(W_fus/H_fus)^0.04832*(L_fus*3.2808)^1.6566;
%平尾
W_ht=453.6*0.1570*(S_ht*3.2808^2)^0.1939;
%垂尾
W_vt=453.6*0.1393*(S_vt*3.2808^2)^0.6729;
%起落架
W_lg=0.07*Wto;
%全機結構
Wstr=W_w+W_fus+W_ht+W_vt+W_lg;
%固定設備
%Wele=220;Weff=500;
Wfix=Wele+Weff;
%推進系統
%Wmot=311;Wcel=332;Wbld=64;
Wthr=Wmot+Wcel+Wbld;
%全機重量
Wto_1=Wstr+Wfix+Wthr;
error=(Wto-Wto_1)/Wto;
if abs(error)<=0.005;
    fprintf('\n機翼重量為 %d g\n',W_w);
    fprintf('機身重量為 %d g\n',W_fus);
    fprintf('平尾重量為 %d g\n',W_ht);
    fprintf('垂尾重量為 %d g\n',W_vt);
    fprintf('起落架重量為 %d g\n',W_lg);
    fprintf('結構重量為 %d g\n',Wstr);
    fprintf('電子設備（操縱系統）重量為 %d g\n',Wele);
    fprintf('有效載荷重量為 %d g\n',Weff);
    fprintf('固定設備重量為 %d g\n',Wfix);
    fprintf('電機、調速器、減速器重量和為 %d g\n',Wmot);
    fprintf('電池組重量為 %d g\n',Wcel);
    fprintf('螺旋槳重量為 %d g\n',Wbld);
    fprintf('推進系統重量為 %d g\n',Wthr);
    fprintf('全機重量為 %d g\n\n',Wto);
else if error<-0.005;
    fprintf('初設重量太輕，請重新確定\n');
end
if error>0.005;
    fprintf('初設重量太重，請重新確定\n');
end
end