# 飛設小工具

###### tags: `飛設`
> [name=F44096128吳宇晟]
## 目錄
> [color=#c4812b]
[TOC]
## 本學期目標
1. 人性化的飛設表格設計
2. 小飛機redesign(製造前)
放隻可愛的小虎鯨

---

## 進度表

| 月份 | 週次(日期) | 內容 |
| ---- | --------------------- | --------------------------------------------------------------------------------------------- |
| Oct | I (Oct.24~Oct.30) | 進度表規劃|
| Nov | II (Oct.31~Nov.6) | 飛設表格設計與公式整理 |
| | III (Nov.7~Nov.13) | 限制分析 |
| | IV (Nov.14~Nov.20) | 期中考 |
| | V (Nov.21~Nov.27) | 主翼設計 |
| Dec | VI (Nov.28~Dec.4) | 安定面與控制面設計&馬達分析 |
| | VII (Dec.5~Dec.11) | 縱向穩定分析 |
| | VIII (Dec.12~Dec.18) | 阻力分析 |
| | IX (Dec.19~Dec.26) | Catia 氣動力外型設計 |
| | X (Dec.26~Jan.1) | Catia 結構設計 |
| Jan | XI (Jan.2~Jan.8) | 期末考 |

---

## 公式整理 (附 matlab code)
### 限制分析
P.S.水平加速可以先不用考慮
#### 小註記 :
* k : 1/πeAR
(e 為奧斯華常數，詳細請看阻力設計)
* n : 附載因子(√2)
* Vv : 爬升速度
* μ(TO) : 摩擦係數(0.1)
#### :warning: 注意 : W單位為(N)，不是(KG)
#### :warning: 注意 : 起飛距離q中的速度是起飛速度不是巡航速度
#### :warning: 注意 : 降落距離q中的速度是失速速度不是巡航速度
#### 小註記 :
* η : 動力整體效率 ( 馬達η * 槳η )

#### 這裡要做的事 :
1. 畫出推重比(T/W)與翼負載(W/S)的關係圖
2. 畫出馬達功率(Power)與翼負載(W/S)的關係圖
3. 找出你的飛機設計點

#### 限制分析matlab code:
等速巡航 :

function [TW]=Constant_Speed_Cruise(WS, q,Cd,k)
TW=q*Cd./WS+k/q.*WS
end
%----------------------------------------------------
% Total 4 inputs :
%  WS : wing loading
%  q : Dynamic Pressure
%  Cd : drag coefficient
%  k : 1/πeAR

---

等速爬升 :

function [TW]=Constant_Speed_Climb(WS, q,Cd,k,Vv,V)
TW=Vv/V+q*Cd./WS+k/q.*WS
end
%----------------------------------------------------
% Total 6 inputs :
%  WS : wing loading
%  q : Dynamic Pressure
%  Cd : drag coefficient
%  k : 1/πeAR
%  Vv : vertical speed
%  V : curise speed

---

等速轉彎 :

function [TW]=Constant_Speed_Level_Turn(WS, q,Cd,k,n)
TW=q*(Cd./WS+k*(n/q)^2.*WS) 
end
%----------------------------------------------------
% Total 4 inputs :
%  WS : wing loading
%  q : Dynamic Pressure
%  Cd : drag coefficient
%  k : 1/πeAR
%  n : 附載因子(通常為 sqrt(2))

---

起飛距離 :

function [TW]=Takeoff_Distance(WS,q,V_to,d,Cd_to,Cl_to,mu)
TW=((V_to^2)/(2*9.81*d))+((q*Cd_to)./WS)+(mu*(1-q*Cl_to./WS));
end

%----------------------------------------------------
% Total 6 inputs :
%  V_to : takeoff speed
%  d : takeoff distance
%  q : Dynamic Pressure
%  Cd_to : takeoff drag coefficient
%  Cl_to : take off lift coeficient
%  WS : wing loading
%  mu : 摩擦係數(通常為0.1)

---

降落速度 :

function WS=Approach_Speed(lo,V_stall,cl_max)
WS=0.5*lo*V_stall^2*cl_max
WS=round(WS,2)
%-----------------------------------------------------
% Total 3 inputs :
%  lo : cruise air density
%  V_stall : stall speed
%  cl_max : max lift coefficient

---

限制分析(運用前方子方程式) :

function [W_S,T_W,P_motor]=Constraint_Analysis(W,AR,lo_air,lo_cruise,V_cruise,V_stall,Vv,cl_max,Cd,Cd_to,Cl_to,d)

%設定數值
WS=24.5:0.1:200;  %------------設定翼負載區間
n=sqrt(2);        %------------附載因子
mu=0.1;           %------------μ摩擦係數
eta=0.6*0.8;      %------------η動力整體效率

e=1.78*(1-0.045*AR^0.68)-0.64;    %------------奧斯華常數
k=1/(pi()*e*AR);                  
q=0.5*lo_cruise*(V_cruise)^2;
V_to=V_stall*1.3;

foldername='D:\NCKU\UAV\matlab\Chart';   %---------圖片儲存路徑

%%
%畫出推重比與翼負載間的關係圖
%--------------------------------------------------Plot Cruise T/W v.s. W/S
figure(1);
Cruise=Constant_Speed_Cruise(WS,q,Cd,k);   
plot(WS,Cruise,'black');
hold on;
%--------------------------------------------------Plot Climb T/W v.s. W/S
Climb=Constant_Speed_Climb(WS,q,Cd,k,Vv,V_cruise);   
plot(WS,Climb,'blue');
hold on;
%--------------------------------------------------Plot Turn T/W v.s. W/S
Level_Turn=Constant_Speed_Level_Turn(WS,q,Cd,k,n);
plot(WS,Level_Turn,'r');
hold on;
%--------------------------------------------------Plot Takeoff T/W v.s. W/S
Takeoff=Takeoff_Distance(WS,q,V_to,d,Cd_to,Cl_to,mu);
plot(WS,Takeoff,'green');
hold on;
%--------------------------------------------------Plot Landing T/W v.s. W/S
Landing=Approach_Speed(lo_air,V_stall,cl_max);
WS1=Landing;
K=Landing*ones(length(WS),1);
plot(K,linspace(0.05,0.35,length(K)));
hold on;
%--------------------------------------------------Find design doint and plot
W_S=floor(WS1);
TW1=Constant_Speed_Cruise(WS1,q,Cd,k);
TW2=Constant_Speed_Climb(WS1,q,Cd,k,Vv,V_cruise);
TW3=Constant_Speed_Level_Turn(WS1,q,Cd,k,n);
TW4=Takeoff_Distance(WS1,q,V_to,d,Cd_to,Cl_to,mu);
A=[TW1 TW2 TW3 TW4];
T_W=ceil(max(A)*1000)/1000;
Design_point=W_S*ones(length(WS),1);
figure(1);
plot(Design_point,T_W,'p','MarkerFaceColor','red','MarkerSize',15);
hold on;
%--------------------------------------------------標題與圖示標記
xlabel('W/S');
ylabel('T/W');
legend('Cruise','Climb','Level Turn','Takeoff','Landing','Design point');
title('Constraint Analysis (T/W v.s.W/S)');
hold off;
saveas(figure(1),fullfile(foldername,'W_S_T_W'), 'png');
%%
%畫出推重比與翼負載間的關係圖
figure(2);
W=W*9.81;     %-----KG to N
%--------------------------------------Plot Cruise P_motor v.s. W/S
P_Cruise=Cruise.*W*V_cruise/eta;   
plot(WS,P_Cruise,'black');
hold on;
%--------------------------------------Plot Climb P_motor v.s. W/S
P_Climb=Climb.*W*V_cruise/eta;    
plot(WS,P_Climb,'blue');
hold on;
%--------------------------------------Plot Turn P_motor v.s. W/S
P_Level_Turn=Level_Turn.*W*V_cruise/eta;   
plot(WS,P_Level_Turn,'r');
hold on;
%--------------------------------------Plot Takeoff P_motor v.s. W/S
P_Takeoff=Takeoff.*W*V_to/eta;
plot(WS,P_Takeoff,'green');
hold on;
%--------------------------------------Plot Landing P_motor v.s. W/S
plot(K,linspace(50,400,length(K)));
hold on;
%--------------------------------------Find design doint and plot
B=[TW1*W*V_cruise/eta TW2*W*V_cruise/eta TW3*W*V_cruise/eta TW4*W*V_to/eta];
P_motor=ceil(max(B));
plot(Design_point,P_motor,'p','MarkerFaceColor','red','MarkerSize',15);
hold on;
%--------------------------------------------------標題與圖示標記
xlabel('W/S');
ylabel('P motor');
legend('Cruise','Climb','Level Turn','Takeoff','Landing','Design point');
title('Constraint Analysis (P motor v.s.W/S)');
saveas(figure(2),fullfile(foldername,'P_motor_W_S'), 'png');
hold off;

---

#### 可求出 :
1. 設計點的T/W、motor P、W/S

2. 推重比(T/W)與翼負載(W/S)的關係圖(下方示意圖)
:star: : 設計點


3. 馬達功率(Power)與翼負載(W/S)的關係圖(下方示意圖)


---
#### :100: 恭喜做完第一步了 ! ! ! ! 接下來請繼續加油
---

### 主翼設計第一階段
#### 這裡要來找機翼的基本需求參數



#### 這裡要做的事 :
1. 找出在巡航時預估的二維升力係數
2. 找出巡航時的雷諾數
3. 找出機翼的wing span 與 chord line
4. 透過1&2的數值去挑選符合的翼型 ( 攻角小於五度時升力係數要大於預估值 )

---

#### 主翼設計第一階段 matlab code：

function [S_ref,b,c,Cl_wing,Re_cruise]=Wing_Design1(AR,W_S,W,lo_cruise,V_cruise,mu)
W=W*9.81;           %kg to N
S_ref=W/W_S;
c=sqrt(S_ref/AR);
MAC=c;
b=AR*c;
CL_cruise=W/(0.5*lo_cruise*V_cruise^2*S_ref);
Cl_wing=CL_cruise/0.9/0.9;
Re_cruise=(lo_cruise*V_cruise*(MAC))/mu;

%--------------------------------------------------
% Total 6 inputs :
%  AR : Aspect ratio
%  W_S : wing loading
%  W : weight(kg)
%  lo_cruise : cruise air density
%  mu : viscosity

---

### 主翼設計第二階段
#### 這裡已經要找好翼型了


#### :warning: 注意 : 這裡有關角度的所有東西都是弧度(斜率也是除弧度)

#### 這裡要做的事 :
1. 找出二維升力斜率與二維力矩斜率(後面N.P.點要用)
2. 算出三維升力係數

---

#### 主翼設計第二階段 matlab code：

function [CLa,CLmax,CMa]=Wing_Design2(AR,Cla,Clmax,Cma)

kappa=Cla/(2*pi());           
CLa=(2*pi()*AR)/(2+sqrt(4+(AR/kappa)^2));    % 三維升力斜率
CLmax=Clmax*0.9;
CMa=CLa*Cma/Cla;
end
%--------------------------------------------------
% Total 4 inputs :
%  AR : Aspect ratio
%  Cla : 二維升力斜率(rad)
%  Clmax : 最大升力係數
%  Cma : 二維力矩斜率(rad)

---

#### 完蛋，控制期中寄了 :cry:

---

### 尾翼設計

#### 體積係數 :
1. 水平尾翼體積係數公式 :


2. 垂直尾翼體積係數公式 :

#### :warning: 注意 : 力臂因為還沒有重心，所以改為主翼氣動力中心

體積係數參考值(:accept:看這邊) :
| Type | l_HT | L_VT |
| --------------- |:----:|:----:|
| 滑翔翼 | 0.5 | 0.2 |
| 手做飛機 | 0.5 | 0.04 |
| 民用單引擎 | 0.7 | 0.04 |
| 民用雙引擎 | 0.8 | 0.07 |
| 農用 | 0.5 | 0.04 |
| 雙渦輪螺槳 | 0.9 | 0.08 |
| 噴射教練機 | 0.7 | 0.06 |
| 噴射戰鬥機 | 0.4 | 0.07 |
| 軍用運輸/轟炸機 | 1 | 0.08 |
| 民用運輸機 | 1 | 0.09 |
| 飛行艇(水用) | 0.7 | 0.06 |

#### 控制面 :
參照下方表抓大概來設計 :8ball:

副翼 :


升降舵(水平尾翼的) :


方向舵(垂直尾翼的) :


#### 這裡要做的事 :
1. 透過體積係數參考表與公式，求出水尾面積與垂尾面積
2. 透過控制面參考表，找出適合的控制面設計

---

#### 求尾翼面積 matlab code :

function [S_HT,S_VT]=Tail_area(V_H,V_V,S_wing,b,c,l_HT,l_VT)
S_HT=V_H*S_wing*c/l_HT;
S_VT=V_V*S_wing*b/l_VT;
end
%--------------------------------------------------
% Total 7 inputs :
%  V_H : Horizontal Tail Volume Coefficient
%  V_V : Vertical Tail Volume Coefficient
%  S_wing : Wing area
%  b : Wing span
%  c : Wing chord line length
%  l_HT : 主翼氣動力中心到水尾氣動力中心的距離
%  l_VT : 主翼氣動力中心到垂尾氣動力中心的距離

---

### 縱向穩定分析


#### 小註記 :
* η : 約為0.9
* 𝜕𝜖/𝜕𝛼 : (2×CL,α) / (AR×π)
#### :warning: 注意 : X(bar)→長度/M.A.C長
#### :warning: 注意 : 全部都是三維(升力、力矩)斜率喔

其中機身三維力矩斜率公式如下 :


#### 小註記 :
* 機身長度為總長 ( 包含碳管 )

:memo: Kf參考表如下 :

| 主翼c/4位置於機身總長之比例 | Kf |
|:------------ | ---- |
| 0.1 | 0.115|
| 0.2 |0.172 |
| 0.3 |0.344 |
| 0.4 |0.487 |
| 0.5 |0.688 |
| 0.6 |0.888 |
| 0.7 |1.146 |
#### 這裡要做的事 :
1. 透過公式，求出整架飛機的中性點
2. 設定穩定裕度，找出飛機重心
3. 開始設計機身參數

---

#### 找N.P.點 matlab code :
:warning: S.M. 預設值為13%，要改請至程式改

function [NP,CG,CMa_f]=Longitudinal_Stability_Analysis(AR,CLa_w,CLa_t,c,l_HT,S_HT,S_wing,CMa_w,CMa_t,Kf,Wf,Lf)

SM=0.13;     % 穩定域度預設值為13%
eta=0.9;
Xwac=c/4;    % 主翼翼前緣到主翼A.C.距離
Xtac=Xwac+l_HT;     % 主翼翼前緣到尾翼A.C.距離

CMa_f=(Kf*Wf^2*Lf)/(S_wing*c);      % 機身三維力矩斜率(rad)
downwash_effect=(2*CLa_w)/(pi()*AR); % 𝜕𝜖/𝜕𝛼

NP=(CLa_w*(Xwac/c)+eta*(S_HT/S_wing)*CLa_t*(1-downwash_effect)*(Xtac/c)-(CMa_w+CMa_f+CMa_t))/(CLa_w+eta*(S_HT/S_wing)*CLa_t*(1-downwash_effect));
CG=NP-SM;
end
%--------------------------------------------------
% Total 13 inputs :
%  AR : Aspect ratio
%  CLa_w : 主翼三維升力斜率
%  CLa_t : 尾翼三維升力斜率
%  c : M.A.C
%  l_HT : 主翼氣動力中心到水尾氣動力中心的距離
%  S_HT : 水尾面積
%  S_wing : Wing area
%  CMa_w : 主翼三維力矩斜率(rad)
%  CMa_t : 尾翼三維力矩斜率(rad)
%  SM : 穩定裕度

% ----- ↓ 以下算機身力矩斜率用 ↓ -----
%  Kf : 機身位置修正值(查表)
%  Wf : 機身最大寬度
%  Lf : 機身長度

### 阻力分析

飛機的阻力以共分成以下幾種 :


---

#### 1. 摩擦/形狀阻力 :

總體公式如下，各部件依照不同計算方式而有不同結果 :


* #### CF : 分成層流邊界層與紊流邊界層兩種計算方式 :


而飛機在飛行時邊界層不會全部都一樣，大概只有前10%~20%是層流，所以我們會分配比例去做加成，公式如下 :


#### 小註記 :
* k 約為0.1
* ρ 為巡航空氣密度
##### :warning: 注意 : 每個部件的雷諾數不一樣，例如機身c是指總長，主翼c是指平均氣動力弦長

* #### FF : 每個部件也會有不同的形狀修正係數，需要個別分開來求 :
#### 機翼 ( 主翼尾翼分開求 ) :


#### 小註記 :
* t/c : 厚弦比，指機翼最厚的地方為弦長幾%
* x/c : 厚度最厚位置佔弦長的比例(%)，(有控制面之尾翼，會因轉軸間隙而增加約10%)
* For 平板翼 : t/c=厚度/弦長，x/c=1

#### 機身 : (根據雷諾數不同代不同公式)


* #### S_wet : 濕面積


S_exp 可由繪圖軟體求得，有暴露在外面的都算
#### :warning: 注意 : 濕面積不是暴露面積，S_exp才是

* #### 干擾阻力修正值 :


#### →CDF(total)=ΣCDF(i)×Q(i)，每一項都要乘對應修正值

---

#### 2. 雜項阻力 :

最後有以下情況的再加上雜項阻力，有再加，隨你高興(起落架、輪胎可加可不加)


---

#### 3. 寄生阻力 :

把剛剛算的東西加一加就變成了我們的CD0了，讚啦 !:+1:


---

#### 4. 誘導阻力 :

誘導阻力為機翼產生升力時的下洗流所伴隨產生的阻力 → 跟你的翼型、攻角有關係

公式如下 :


其中 e 為奧斯華常數 :


---
#### 5. 總阻力
:wc: 全部加一加就對了


---

#### 這裡要做的事 :
1. 找出各項阻力係數，最後求出CD0
2. 機身的詳細數據要出來

---

#### 阻力分析 matlab code :

M=V_cruise/340;     % 預設聲速為340m/s
k_w=0.12;           % 預設機翼層流只發生在前12%
k_f=0.1;            % 預設機身層流只發生在前10%
xlsFile = 'D:\NCKU\UAV\matlab\Aircraft_design.xlsx';   % 參數匯出路徑

%%  fuselage

Re_f=Re_w*Lf/c;    
CF_f_L=1.328/sqrt(Re_f);
CF_f_T=0.455/(log10(Re_f)^2.58*(1+0.144*M^2)^0.65);
CF_f_avg=k_f*CF_f_L+(1-k_f)*CF_f_T;

f=Lf/sqrt(Amax_f*4/pi());
FF_f=(1+60/f^3+f/400);

S_wet_f=pi()*Df*Lf*(0.5+0.135*Ln/Lf)^(2/3)*(1.015+0.3/(Lf/Df)^1.5);

CDF_f=CF_f_avg*FF_f*S_wet_f/S_ref;

%%  wing

CF_w_L=1.328/sqrt(Re_w);
CF_w_T=0.455/(log10(Re_w)^2.58*(1+0.144*M^2)^0.65);
CF_w_avg=k_w*CF_w_L+(1-k_w)*CF_w_T;

FF_w=(1+0.6*(thickness_w)/(thick_position_w)+100*(thickness_w)^4);

S_wet_w=(2+0.5*thickness_w)*S_exp_w;

CDF_w=CF_w_avg*FF_w*S_wet_w/S_ref;
%%  tail

Re_t=Re_w*Lf/c;
CF_t_L=1.328/sqrt(Re_t);
CF_t_T=0.455/(log10(Re_t)^2.58*(1+0.144*M^2)^0.65);
CF_t_avg=k_w*CF_t_L+(1-k_w)*CF_t_T;

FF_t=(1+0.6*(thickness_t)/(thick_position_t+0.1)+100*(thickness_t)^4);

S_wet_t=(2+0.5*thickness_t)*S_exp_t;

CDF_t=CF_t_avg*FF_t*S_wet_t/S_ref;

%%  total parasite drag

CDF_total=CDF_f*Qf+CDF_w*Qw+CDF_t*Qt;
CD_upsweep=Amax_f/S_ref*3.83*deg2rad(u)^2.5;
CD0=CDF_total+CD_upsweep;

%%  export data to excel (選配)

excel_matrix=[CF_w_L,CF_t_L,CF_f_L;CF_w_T,CF_t_T,CF_f_T;k_w,k_w,k_f;CF_w_avg,CF_t_avg,CF_f_avg;FF_w,FF_t,FF_f;S_wet_w,S_wet_t,S_wet_f;Qw,Qt,Qf;CDF_w,CDF_t,CDF_f];
writematrix(round(excel_matrix,3),xlsFile,'Sheet','Results','Range','J4:L11');
writematrix(round([CD_upsweep;CD0],3),xlsFile,'Sheet','Results','Range','J12:J13');

end
%--------------------------------------------------
% Total 18 inputs :
%  V_cruise : Cruise Speed
%  S_exp_w : 主翼暴露在空氣的面積
%  S_exp_t : 尾翼暴露在空氣的面積
%  S_ref : 參考翼面積
%  Re_w : 主翼雷諾數
%  c : 主翼弦長
%  thickness_w : 主翼厚弦比(%)
%  thick_position_w : 主翼最厚的地方位於主翼的位置(%)
%  thickness_t : 尾翼厚弦比(%)
%  thick_position_t : 尾翼最厚的地方位於尾翼的位置(%)
%  Lf : 機身總長
%  Df : 機身等效直徑
%  Ln : 機鼻長度
%  A_max_f : 機身最大截面積
%  u : 機身尾段內縮角度(rad)
%  Qf : 機身干擾阻力修正值
%  Qw : 主翼干擾阻力修正值
%  Qt : 尾翼干擾阻力修正值 

---

### 巡航分析
透過巡航分析我們希望找出最小阻力之速度與最大升阻比之巡航攻角，在有限的電池電量下達到最有效率的飛行。
#### 這裡要做的事 :
1. 找出尾翼的裝置角
2. 找出巡航攻角
3. 找出主翼下洗角
4. 畫出L/D ( 升阻比 ) 的圖
5. 找出最大升阻比
---

#### 巡航(L/D)分析 matlab code :

function [V_D_min,D_min,angle_cruise,tail_angle,downwash_angle,L_D,e1,e2]=LD_Analysis(CD0,CLa_w,CLa_t,a0_w,CG,c,CMa_w,CMa_t,CMa_f,St,Sw,AR,ARt,l_HT,lo_cruise,W);

angle=-1:0.01:8;
a=deg2rad(angle);     % deg to rad
a0_w=deg2rad(a0_w);
eta=0.9;
e1=1.78*(1-0.045*AR^0.68)-0.64;      % 主翼奧斯華常數
e2=1.78*(1-0.045*ARt^0.68)-0.64;     % 尾翼奧斯華常數

% 計算速度與阻力之關係

figure(3);            % 畫速度與阻力的關係圖
    V=5:0.1:30;
    D0=0.5*lo_cruise*V.^2*Sw*CD0;      % 寄生阻力
    k=(1./(pi()*e1*AR)+(St^2./Sw^2)./(pi()*e2*ARt));       % 1/pi*e*AR
    Di=0.5.*lo_cruise.*V.^2.*Sw.*k.*(W*9.81./(0.5*lo_cruise.*V.^2*Sw)).^2;     %誘導阻力
    D_total=D0+Di;
    plot(V,D0,"blue");
    hold on;
    plot(V,Di,"red");
    hold on;
    plot(V,D_total,"black");
    xlabel('Speed(m/s)');
    ylabel('Drag(N)');
    legend('Parasite Drag(N)','Induced Drag(N)','Total Drag(N)');
    title('Speed v.s. Drag');

[value1,index1]=min(D_total);       %找出最低阻力
V_D_min=V(index1);
D_min=value1;

% 計算L/D最大值

for i=1:length(angle) 

        %(下三行)簡化公式
        CLw=CLa_w*(a(i)-a0_w);       %  主翼升力
        epsilon=2*CLw/(pi()*AR);     %  下洗角
        l=(CG*c-(c/4+l_HT))/c;       %  重心到尾翼氣動力中心距離/c

        ai=-(CLw*(CG*c-c/4)/c+CMa_w*a(i)+CMa_t*(a(i)-epsilon)+CMa_f*a(i)+eta*St/Sw*l*CLa_t*(a(i)-epsilon))/(CMa_t+eta*St/Sw*CLa_t*l);
             %  裝置角公式
        CLt=CLa_t*(a(i)-epsilon+ai);

        for k=1:length(angle)        %  找出當下裝置角的升阻比

             CLw1=CLa_w*(a(k)-a0_w);               %  尾翼升力係數
             epsilon1=2*CLw1/(pi()*AR);            %  下洗角
             CLt1=CLa_t*(a(k)-epsilon1+ai);        %  尾翼升力係數
             CD=CD0+CLw1^2/(pi()*e1*AR)+(St^2/Sw^2)*CLt1^2/(pi()*e2*ARt);   %  阻力係數
             LD(k)=(Sw*CLw1+eta*St*CLt1)/(Sw*(CD));             %  升阻比公式

        end

        [value,index]=max(LD);         %  找出最大升阻比

        if index==i         %  判斷最大值是否與巡航攻角相等
            
            figure(4);      %  畫L/D圖
            plot(angle,LD);
            xlabel('Angle of attack (deg)');     
            ylabel('L/D');
            title('L/D v.s. AOA');
            hold off;

            angle_cruise=angle(i);        %  巡航攻角
            tail_angle=rad2deg(ai);       %  裝置角
            downwash_angle=rad2deg(epsilon);   %  主翼下洗角
            L_D=value;
                      
        end
    end

foldername='D:\NCKU\UAV\matlab\Chart';  % 圖片匯出路徑
saveas(figure(3),fullfile(foldername,'Speed_vs_Drag'),'png');
saveas(figure(4),fullfile(foldername,'L_D_vs_AOA'),'png');

end

---
#### 用以上程式可求出 :
1. 阻力與速度之關係圖
2. 巡航速度
3. 最大升阻比
4. 巡航攻角
5. 裝置角
6. 下洗角
7. L/D曲線圖



---
### 電池分析
為了得知完成整個飛行任務所需要的電池重量，須根據各個飛行任務階段來計算出各階段所造成的功率損耗，並再由此計算途中出所需的能源。在這些飛行階段中，最主要耗費能源的任務階段分別為爬升以及巡航，而其他階段如下降、落地由於功率消耗較小，以及起飛由於時間過短因此不列入計算。
1. 巡航所需能量

2. 爬升所需能量

3. 電池重量計算

---
#### 電池分析 matlab code :

Fly_eta=0.8*0.6;               % 修正值
Battery_energy_density=130;    % 電池能量密度
Margin=1.2;                    % 電池裕度
E_cruise=W*g/(L_Dmax)*V_cruise/Fly_eta*Cruise_time/60;               % 巡航所需能量
E_climb=(W*g/(L_Dmax)*(V_cruise-Vv)+W*g*Vv)/Fly_eta*Climb_time/60;   % 爬升所需能量
W_battery=(E_cruise+E_climb)*Margin/Battery_energy_density;          % 電池所需能量

---
#### :warning: 注意 : 在算W時的單位是牛頓(N)而不是公斤(kg)
#### 小註記 :
* η : 約為0.48 (𝜂_𝑝𝑟𝑜𝑝*𝜂_𝑚𝑜𝑡𝑜𝑟)
#### 這裡要做的事 :
1. 找出巡航與爬升所需的最小能量
2. 找出電池重量
---
### 起降分析

---
## 非常簡單的飛設 Excel & Matlab User guide

---