% No.04
% 飛行性能計(jì)算－起飛性能
% 確定飛機(jī)的起飛性能，包括滑跑距離、起飛距離、起飛時(shí)間、離地速度
% 氣動(dòng)修正考慮了收放起落架、襟翼的影響，但未考慮地面效應(yīng)
% 安全高度按15m計(jì)算，起飛安全速度取1.2Vld

function qfxn
clc;
clear;
close;

global m;
global Pky;
global Cx;
global Cy;
global D0;
global S;
global f;
global Theta;
global Alpha;


% 起飛質(zhì)量
m = 8150.0;
G = m*9.81;
S = 23.0;
% 干水泥跑道摩擦系數(shù)
f = 0.03;
% 查起降時(shí)氣動(dòng)特性曲線(襟翼25°，起落架放下)
TCxCy = [
    0.0       0.0551
    0.1014    0.0516
    0.2029    0.0553
    0.3007    0.0691
    0.3986    0.0859
    0.5000    0.1084
    0.6014    0.1390
    0.6993    0.1848
    0.7971    0.2400
    0.8986    0.3097]';
CyAlpha = (0.8-0.2)/14.0;

% 在0高度上查可用推力曲線，由于推力變化在小M下比較平緩，可不考慮速度的影響，取平均值
Pky_full = 50000.0; % 全加力狀態(tài)
Pky_min = 40000.0;  % 小加力狀態(tài)
Pky_max = 30000.0;  % 最大狀態(tài)
TPky = [Pky_max, Pky_min, Pky_full];

% 0高度大氣密度、音速
D0 = 1.225;
A0 = 340.294;

% 取前輪離地速度為70m/s
Vql = 70.0;

key = menu('請(qǐng)選擇起飛時(shí)發(fā)動(dòng)機(jī)工作狀態(tài)', '最大', '小加力' , '全加力');
Pky = TPky(key);

% 求滑跑時(shí)間、滑跑距離

% 前輪離地前取滑跑迎角0.33°
Cy = 0.2+CyAlpha*0.33;
Cx = interp1(TCxCy(1,:), TCxCy(2,:), Cy, 'cube');

t11 = quad8(@Fun_T1, 0.0, Vql);
L11 = quad8(@Fun_L1, 0.0, Vql);

% 前輪離地后取滑跑迎角10°
Cy = 0.2+CyAlpha*10.0;
Cx = interp1(TCxCy(1,:), TCxCy(2,:), Cy, 'cube');
% 求離地速度、起飛安全速度
Vld = sqrt(2.0*G/(Cy*D0*S))
Vaq = 1.2*Vld;

t12 = quad8(@Fun_T1, Vql, Vld);
L12 = quad8(@Fun_L1, Vql, Vld);

t1 = t11+t12;
L1 = L11+L12;
disp('滑跑時(shí)間(s)：');
disp(t1);
disp('滑跑距離(m)：');
disp(L1);

% 求上升時(shí)間、滑跑距離

H2 = 0.0;
Theta = 0.0;
Haq = 15.0;
while(H2 < Haq)
    Theta = Theta + 0.001;
    Cy = 2.0*G*cos(Theta)/(D0*Vld.^2*S);
    Alpha = (Cy-0.2)/CyAlpha*pi/180.0;
    Cx = interp1(TCxCy(1,:), TCxCy(2,:), Cy);
    L2 = quad(@Fun_L2, Vld, Vaq);
    H2 = L2*sin(Theta);
end

t2 = quad(@Fun_T2, Vld, Vaq);

disp('上升時(shí)間(s)：');
disp(t2);
disp('上升距離(m)：');
disp(L2);

disp('起飛時(shí)間(s)：');
disp(t1+t2);
disp('起飛距離(m)：');
disp(L1+L2);


% 被積函數(shù)

function y = Fun_T1(V)
	global m;
	global Pky;
	global Cx;
	global Cy;
	global D0;
	global S;
	global f;
	
	y = m./(Pky-Cx*0.5*D0*V.^2*S-(m*9.81-Cy*0.5*D0*V.^2*S)*f);
return

function y = Fun_L1(V)
	global m;
	global Pky;
	global Cx;
	global Cy;
	global D0;
	global S;
	global f;
	
	y = m*V./(Pky-Cx*0.5*D0*V.^2*S-(m*9.81-Cy*0.5*D0*V.^2*S)*f);
return

function y = Fun_T2(V)
	global m;
	global Pky;
	global Cx;
	global Cy;
	global D0;
	global S;
    global Theta;
    global Alpha;
	
	y = m./(Pky*cos(Alpha)-Cx*0.5*D0*V.^2*S-m*9.81*sin(Theta));
return

function y = Fun_L2(V)
	global m;
	global Pky;
	global Cx;
	global Cy;
	global D0;
	global S;
    global Theta;
    global Alpha;
	
	y = m*V./(Pky*cos(Alpha)-Cx*0.5*D0*V.^2*S-m*9.81*sin(Theta));
return