% MATLAB Script to summarize the 2-engine climb performance for Concorde
% using best guess drag polar.
% AAE 190 Introduction to Aerospace Engineering
%

clear all
clf
disp(' '); disp('Start here for 2-engine analysis'); disp(' ')
ifig=0;

% Given quantities and constants
Vindicated=140:10:500; %Indicated airspeed, knots
H=0 % Altitude, feet
W=408000  %lbf
S=3856 %ft^2
B=83.8 % Wing span ft
A=B*B/S % Aspect ratio
rho=.002378 %slugs/ft^3 for standard sea level conditions.
V2=220 % V2 engine out climb speed knots
VzrcOEFGU=193 % knots, Zero rate of climb speed with one engine failed gear up
VzrcOEFGD=205 % knots, Zero rate of climb speed with one engine failed gear down

% Computed quanities
Vtrueknots=Vindicated; %Assumes that the airspeed indicator is perfect 
% and flight at sea level in standard atmosphere.
Vtruefps=Vtrueknots*1.687808; % Convert knots to ft/sec

% Get 2 engine performance with gear down <<<<<<<<<<<<<<<<<<
disp(' '); disp('Start 2 engine gear down performance'); disp(' ')
nengines=2
% We assume that thrust is constant with airspeed. OK for jets.
thrust=nengines*38050*ones(size(Vtruefps));
qbar=.5*rho*Vtruefps.*Vtruefps; % Dynamic pressure, lbf/ft^2
CL=W./(qbar*S); % Lift coefficient in level flight
% Drag Polar
k= 0.31155 
%CD0=0.011876*2
CD0=.036
CD=CD0+CL.*CL*k;
%k4  = 0.31155, CD04 = 0.011876

%Plot Drag polar and CLmax, gear down
ifig=ifig+1; figure(ifig); 
plot(CD,CL)
grid on
xlabel('Drag Coefficient, CD')
ylabel('Lift Coefficient, CL')
Title('Gear DOWN Drag Polar for Concorde')

grid on
VRindicatedknots=195 % knots, Rotation speed indicated at sea level
VRtrueknots=VRindicatedknots  %Assumes that the airspeed indicator is perfect 
% and flight at sea level in standard atmosphere.
VRtruefps=VRtrueknots*1.687808 % Convert knots to ft/sec
Vstallfps=VRtruefps/1.15 % Rough guess of stall speed
CLmax=W/(.5*.002378*Vstallfps*Vstallfps*S)
CDstall=CD0+CLmax*CLmax*k
hold on; plot(CDstall,CLmax,'x')
tstring4=['CLmax = ',num2str(CLmax),' at Vstall= ',num2str(Vstallfps/1.68780),' knots']
text(.02,1.3,tstring4);
tstring6=['k = ',num2str(k),', CD0 = ',num2str(CD0)]
text(.02,1.1,tstring6);
hold off

% Power required
tstring=[num2str(nengines),' engines, gear down']
PR=S*qbar.*CD.*Vtruefps; % ft-lbf/sec

% Power Available
% We assume that thrust is constant with airspeed. OK for jets.
PA=thrust.*Vtruefps; % ft-lbf/sec
%Plot results
ifig=ifig+1; figure(ifig); 
subplot(2,2,1)
plot(Vindicated,PR,'x',Vindicated,PA,'o')
grid on
xlabel('Indicated airspeed, knots')
ylabel('Power, ft-lbf/sec')
legend('Power Required','Power Available')
title(['Concorde Power for ',tstring])
text(102,12000000,'Conclusion:')
text(102, 4000000,'Concorde cannot climb on two eng.!')
text(300,40000000,'Frontside of PR curve')


% Vertical velocity
hdot=(PA-PR)/W; % ft/sec
RC4engGD=hdot*60
%Plot results
subplot(2,2,2)
plot(Vindicated,hdot*60)
grid on
xlabel('Indicated airspeed, knots')
ylabel('dh/dt, ft/min')
%YI = INTERP1(X,Y,XI) interpolates to find YI, the values of
%    the underlying function Y at the points in the vector XI.
V2climboutspeed=220 %knots
RCV2=interp1(Vindicated,hdot*60,V2climboutspeed);
hold on; plot(V2climboutspeed,RCV2,'x'); hold off
hold on; plot(V2climboutspeed,5000,'rx'); hold off
tstring2=['Rate of climb = ',num2str(RCV2),' at V2 = ',num2str(V2climboutspeed),' knots'];
text(150,5500,tstring2)
text(102,-6000,'Conclusion:')
text(102,-8000,'Concorde cannot climb on two eng.!')
text(150,2500,'k= 0.31155, CD0=.036')



% Flight path Angle
gamma=asin(hdot./Vtruefps);
%Plot results
subplot(2,2,3) 
plot(Vindicated,gamma*57.3)
grid on
xlabel('Indicated airspeed, knots')
ylabel('Flight path angle, deg')
text(102,-11,'Conclusion:')
text(102,-13,'Concorde cannot climb on two eng.!')


% Drag
Drag=S*qbar.*CD;
%Plot results
subplot(2,2,4)
plot(Vindicated,Drag)
grid on
xlabel('Indicated airspeed, knots')
ylabel('Drag, lbf')
minDrag=min(Drag)
ii=find(Drag==minDrag);
VminDrag1=Vindicated(ii)
hold on; plot(VminDrag1,minDrag,'x'); hold off
hold on; plot(400,minDrag,'rx'); hold off
tstring3a=['Minimum Drag = ',num2str(minDrag),' lbf.'];
tstring3b=['at VminDrag = ',num2str(VminDrag1),' knots'];
text(150,125000,tstring3a)
text(150,105000,tstring3b)

