% Alternative solution to Concorde Drag Polar
clear all
disp(' '); disp('Start here'); disp(' ')
ifig=0;

% Given quantities and constants
H=0 % Altitude, feet
W=408000  %lbf
S=3856 %ft^2
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

Vtrueknots=220                %knots
Vtruefps=Vtrueknots*1.687808 % Convert knots to ft/sec

% Get 4 engine performance with gear down <<<<<<<<<<<<<<<<<<
nengines=4
% We assume that thrust is constant with airspeed. OK for jets.
thrust=nengines*38050
qbar=.5*rho*Vtruefps*Vtruefps % Dynamic pressure, lbf/ft^2
hdot=5000/60 % ft/sec
CD1=2*(thrust-W*hdot/Vtruefps)/(rho*Vtruefps*Vtruefps*S)   %eqn 1
CL1=W/(qbar*S) % Lift coefficient in level flight          %eqn 2


Vmindragknots=400 %knots
Vmindragfps=Vmindragknots*1.687808 % Convert knots to ft/sec

c=(Vmindragfps^4)*rho*rho*S*S/(4*W*W)                      %eqn 5
CD0=CD1/(1+c*CL1*CL1)                                      %eqn 6 
k=CD0*c                                                    %eqn 7


CL2=0:.1:1.2
CD2=CD0+CL2.*CL2*k;

%Plot Drag polar, gear down
ifig=ifig+1; figure(ifig); 
plot(CD2,CL2,'-');                           % plot with solid line
tstring7=[' k  = ',num2str(k),', CD0 = ',num2str(CD0),' solid line, match 5000fpm and Vmindrag, gear down']
z=axis;
z(1)=-.4;
text(z(1),z(4)-.1,tstring7);
grid on
xlabel('Drag Coefficient, CD')
ylabel('Lift Coefficient, CL')
Title('Drag Polar for Concorde')

%CLmindrag2=W/(.5*.002378*Vmindragfps*Vmindragfps*S)
%CDmindrag2=CD0+CLmindrag2*CLmindrag2*k
%hold on; plot(CDmindrag2,CLmindrag2,'rx')
%dstring4=[' CD2 = ',num2str(CDmindrag2),',  CL2 = ',num2str(CLmindrag2),' at Vmindrag= ',num2str(Vmindragfps/1.68780),' knots']
%text(CDmindrag2,CLmindrag2,dstring4);

hold on; plot(CD1,CL1,'rx'); hold off
dstring6=[' CD1 = ',num2str(CD1),',  CL1 = ',num2str(CL1),' at V5000= ',num2str(Vtrueknots),' knots']
text(CD1,CL1,dstring6);


% Analyze  zero rate of climb case (gear down, 3 engines)
CDzrc3=2*3*38050/(rho*(VzrcOEFGD*1.687808)^2*S)               %eqn 8
CLzrc3=2*W/(rho*(VzrcOEFGD*1.687808)^2*S)                     %eqn 9
dstring8=[' CD3 = ',num2str(CDzrc3),',   CL3 = ',num2str(CLzrc3),' at VzrcOEFGD= ',num2str(VzrcOEFGD),' knots']
text(CDzrc3,CLzrc3,dstring8);
hold on; plot(CDzrc3,CLzrc3,'rx'); hold off

% Do least squares fit on two data points (4 engine 5000 ft/min RC and 3 engine zrc)
y=[CD1,CDzrc3]
x=[CL1*CL1,CLzrc3*CLzrc3]
[k3,CD03]=straightline(x,y)
CL3=0:.1:1.2
CD3=CD03+CL3.*CL3*k3;
hold on; plot(CD3,CL3,':'); hold off                    % plot with dotted line
tstring10=[' k3  = ',num2str(k3),', CD03 = ',num2str(CD03),' dotted line, match 5000fpm and VzrcOEF, gear down']
text(z(1),z(4)-.2,tstring10);

% Using min drag constraint, Do least squares fit on two data points (4 engine 5000 ft/min RC and 3 engine zrc)
[m]=straightline2(x,y)
CD04=m/(1+c)                                                   %eqn 10
k4=c*CD04                                                      %eqn 7
CL4=0:.1:1.2
CD4=CD04+CL4.*CL4*k4;
hold on; plot(CD4,CL4,'-.'); hold off                    % plot with dash dot line
tstring11=[' k4  = ',num2str(k4),', CD04 = ',num2str(CD04),' dash dot line,  match 5000fpm and VzrcOEF and Vmindrag, gear down']
text(z(1),z(4)-.3,tstring11);

% Use only the zero rate of climb data point and the Vmindrag=400 constraint
% Using min drag constraint, Do least squares fit on two data points (4 engine 5000 ft/min RC and 3 engine zrc)
[m]=straightline2(x(2),y(2))  % Note only the second data point of the vectors x and y are used.
CD05=m/(1+c)                                                   %eqn 10
k5=c*CD05                                                      %eqn 7
CL5=0:.1:1.2
CD5=CD05+CL5.*CL5*k5;
hold on; plot(CD5,CL5,'--'); hold off                 % plot with  dashed line
tstring12=[' k5  = ',num2str(k5),', CD05 = ',num2str(CD05),' dashed line, match VzrcOEF and Vmindrag, gear down']
text(z(1),z(4)-.4,tstring12);

% Analyze  zero rate of climb case (gear up, 3 engines)
CDzrc4=2*3*38050/(rho*(VzrcOEFGU*1.687808)^2*S)                  %eqn 8
CLzrc4=2*W/(rho*(VzrcOEFGU*1.687808)^2*S)                        %eqn 9
dstring13=[' CD4 = ',num2str(CDzrc4),',  CL4 = ',num2str(CLzrc4),' at VzrcOEFGU= ',num2str(VzrcOEFGU),' knots']
text(CDzrc4,CLzrc4,dstring13);
hold on; plot(CDzrc4,CLzrc4,'ro'); hold off

% Plot drag polar for Boeing 707
CD06=.0131                                                   %eqn 10
k6=1/(pi*7.4*.85)                                                      %eqn 7
CL6=0:.1:1.2
CD6=CD06+CL6.*CL6*k6;
hold on; plot(CD6,CL6,'g-'); hold off                 % plot with  green solid line
tstring13=[' k6  = ',num2str(k6),', CD06 = ',num2str(CD06),' green solid line, Boeing 707-320B']
text(z(1),z(4)-.5,tstring13);
% Ref: Lan and Roskam p.190
LoverDmax=max(CL6./CD6)

% Plot best guess drag polar, gear down
CD07=.036                                                     
k7=.31155                                                  
CL7=0:.1:1.2
CD7=CD07+CL7.*CL7*k7;
hold on; plot(CD7,CL7,'m-'); hold off                 % plot with  magenta solid line
tstring14=[' k7  = ',num2str(k7),', CD07 = ',num2str(CD07),' magenta solid line, best guess, gear down']
text(z(1),z(4)-.6,tstring14);
% Plot best guess drag polar, gear up
CD08=.016                                                     
k8=.31155                                                  
CL8=0:.1:1.2
CD8=CD08+CL8.*CL8*k8;
hold on; plot(CD8,CL8,'c-'); hold off                 % plot with  cyan solid line
tstring15=[' k8  = ',num2str(k8),', CD08 = ',num2str(CD08),' cyan solid line, best guess, gear up']
text(z(1),z(4)-.7,tstring15);
