% liftdragpolar.m
% 
% LIFTDRAGPOLAR Calculate the lift drag polar for the Aerosonde UAV 
%               using Chebfuns
%
% References:
% J. F. Whidborne, ``Solving optimal control problems using the Chebfun
%    system'', UKACC Control 2016, Belfact, 2016.
% G.~J. Holland, T.~McGeer, and H.~Youngren, ``Autonomous Aerosondes for
%   economical atmospheric soundings anywhere on the globe,'' Bulletin of
%   the American Meteorological Society, vol.~73, no.~12, pp. 1987-1998, 1992. 
% Unmanned Dynamics, AeroSim Blockset User's Guide, Unmanned Dynamics,
%   LLC, 8 Fourth St., Hood River, OR., 2004, version 1.2.
%

% (c) James Whidborne, Cranfield University, 15 March 2016
% MATLAB   Version 8.6.0.267246 (R2015b)
% Chebfun  Version 5.3.0  


%% Basic data 
% Flight conditions
m = 13.5;     % mass (kg) % 13.5 kg for Aerosonde;
g = 9.81;     % gravitational constant (m/s^2)
rho = 1.225;  % mass density of air

% constants for Aerosonde
CD0 = 0.0434; % Minimum drag coefficient
CLmind = 0.23;% Lift at minimum drag
osw = 0.75;   % Oswalds coefficient (wing efficiency factor)
b = 2.8956;   % Wind span (m)
S = 0.55;     % Wing area (m^2)
AR = b^2/S;   % wing aspect ratio
 
%% Calculate lift drag polar
% airspeed over interval (m/s^2)
V=chebfun('V',[15,60]);
% lift coefficient 
CL = 2*m*g./(rho*S*V.^2);  %C_L=\frac{2W}{rho S V^2}, W=mg 
% drag coefficient
CD = CD0 + (CL-CLmind).^2./(pi*osw*AR); % C_D = C_{D0} + \frac{(C_L - C_{L\min_D})^2}{\pi e A\!R}

%% Calculate drag as function of V 
D = (rho*S/2)*CD.*V.^2;
% calculate minimum drag speed
[minD,VminD]=min(D);
disp(['Minimum drag D=',num2str(minD),' N at airspeed V=',num2str(VminD),' m/s'])

%% optimal lift drag ratio calculation
LDR=CL./CD; 
[maxLDR,VmaxLDR]=max(LDR);
disp(['Maximum lift drag ratio L/D=',num2str(maxLDR),' at airspeed V=',num2str(VmaxLDR),' m/s'])

%% plots

% plotting constants
fs=12; % fontsize
lw=1.5; % linewidth
ms=14; % markersize

figure(1)% lift drag polar
plot(CD,CL,'linewidth',lw);
axis auto
xlabel('Lift coefficient, $C_D$','Interpreter' ,'Latex','fontsize',fs); 
ylabel('Drag coefficient, $C_L$','Interpreter' ,'Latex','fontsize',fs);
title('Lift drag polar')
grid on 

figure(2) % drag as function of V 
plot(V,D,'b',VminD,minD,'xr','linewidth',lw,'markersize',ms)
xlabel('Air speed, $V$, (m/s)','Interpreter' ,'Latex','fontsize',fs); 
ylabel('Drag, $D$, (N)','Interpreter' ,'Latex','fontsize',fs); 
title('Drag versus airspeed')  
grid on 

figure(3) % lift drag ratio calculation 
plot(V,LDR,'b',VmaxLDR,maxLDR,'xr','linewidth',lw,'markersize',ms)
xlabel('Air speed, $V$, (m/s)','Interpreter' ,'Latex','fontsize',fs); 
ylabel('Lift drag ratio $L/D$','Interpreter' ,'Latex','fontsize',fs); 
title('Lift drag ratio versus airspeed')  
grid on