% The FDC toolbox. Force/moment coefficient plotting routine PLOTCOEFF.
% =====================================================================
% PLOTCOEFF is a Matlab macro that was used to verify the Aerodynamic and
% propulsive force & moments models by comparing the results to ref.[1].
% (Ref.[1] contained similar data from flightsimulator software by Delft
% University of Technology, and comparison data from windtunnel measure-
% ments.)
%
% PLOTCOEFF plots the dymensional force and moment coefficients for a 
% range of angle of attack and sideslip angle values and it visualizes the 
% influence of the propeller slipstream that is taken into account by the 
% 'dpt' value in the propulsion model (the dimensionless pressure increase
% in the plane of the propeller). 
%
% Since the results also provide additional general insight in the typical
% force and moment characteristics of the 'Beaver' aircraft, the macro 
% was included in the FDC 1.4 distribution. See the source code for details.
%
% Usage: PLOTCOEFF will call the Simulink system aeropropmodel.mdl to
% determine the force and moment coefficients for a range of conditions. 
% NOTE: the macro does not clean up the workspace afterwards and may over-
% write variables already existing in the workspace without warning!
%
% REFERENCE:
% [1] Tjee, R.T.H., Mulder, J.A.: Stability and Control Derivatives of 
%     the De Havilland DHC-2 'Beaver' aircraft.
% ------------------------------------------------------------------------

% PLOTCOEFF uses the Simulink model 'aeropropmodel.mdl' to find the output
% values. A dummy time-axis is used to be able to compute the outputs for
% a complete alpha or beta range without resorting to a nested FOR-loop. 
% The remaining FOR loops take into account three different values of 'dpt'.
%
% Variables:
% ==========
% alpha        angle of attack range (degrees)
% alpharad     angle of attack range (radians)
% AM           parameter matrix for aerodynamic model
% beta         sideslip angle range (degrees)
% betarad      sideslip angle range (radians)
% dpt          non-dimensional pressure increase across the plane of the propeller
% EM           parameter matrix for engine propulsion model
% h_fig        height of figure
% i            counter
% screen_size  coordinates of lower left corner, screenwidth and screenheight
% t            time, returned from the Simulink model (not used)
% time         dummy time-axis
% u            input range (matrix with four columns: time (dummy), alpha column,
%                 beta column, and dpt column)
% w_fig        width of figure
% x            state variables, obtained from the Simulink model (n/a, so empty matrix)
% x_fig        X-coordinate of lower left figure corner
% y_fig        Y-coordinate of lower left figure corner
%
% Force  coefficients along X, Y, and Z-axes: CX, CY, CZ
% Moment coefficients along X, Y, and Z-axes: Cl, Cm, Cn


% Check if AM and EM have been defined in the Matlab workspace. If not, run 
% DATLOAD to load them from file.
% -------------------------------------------------------------------------
if exist('AM')==0 | exist('EM')==0
    h=newMsgBox(['First, the model parameters need to be retrieved from file ', ...
                 '(e.g. AIRCRAFT.DAT). Click ''OK'' or press Enter to continue.']);
    uiwait(h);
    datload;
end

% If model parameters are still not present in the workspace,
% e.g. due to an incorrect datafile, force an abort of PLOTCOEFF.
% ---------------------------------------------------------------
if exist('AM')==0 | exist('EM')==0
    error(['ERROR: the model parameters are still not present in the workspace! ', ...
                 'Aborting PLOTCOEFF.']);
end

% Define plot symbols
% -------------------
symbols={'-',':','b--'};

% dpt range
% ---------
dpt = [0 0.5 1];

% Define alpha and beta axes (use rad values for actual computations and
% the other values as dummy time-axis and for plotting). Dummy time-
% axis is used to obtain all results with one function call only.
% ---------------------------------------------------------------------
time     = [0:1:20]';
alpha    = [0:1:20]';
beta     = [-10:1:10]';
alpharad = pi/180*alpha;
betarad  = pi/180*beta;

% Initialize input range
% ----------------------
u = zeros(21,4);

% Determine suitable figure dimensions
% ------------------------------------
screen_size = screensize;
x_fig = screen_size(1)+100;
y_fig = screen_size(2)+100;
w_fig = screen_size(3)-200;
h_fig = screen_size(4)-200;

% Initiate figure
% ---------------
figure('DefaultAxesPosition',[0.08 0.1 0.86 0.8], ...
     'Name','Force and moment coefficients', ...
     'NumberTitle','off', ...
     'Position',[x_fig,y_fig,w_fig,h_fig]);

% Determine CX, CZ, Cm as a function of alpha for three dpt values
% ----------------------------------------------------------------
for i=1:3;

   % Set input range
   % ---------------
   u = [time alpharad zeros(21,1) dpt(i)*ones(21,1)];

   % Determine outputs (only y is relevant)
   % --------------------------------------
   [t,x,y]=sim('aeropropmodel',time,[],u);

   % Plot results
   % ------------
   subplot(2,3,1); plot(alpha,y(:,1),symbols{i}); 
      grid; hold on; xlabel('\alpha [deg]','FontAngle','italic'); ylabel('C_X','FontAngle','italic');
   subplot(2,3,2); plot(alpha,y(:,3),symbols{i});
      grid; hold on; xlabel('\alpha [deg]','FontAngle','italic'); ylabel('C_Z','FontAngle','italic');
   subplot(2,3,3); plot(alpha,y(:,5),symbols{i}); 
      grid; hold on; xlabel('\alpha [deg]','FontAngle','italic'); ylabel('C_m','FontAngle','italic'); 

end

% Determine CY, Cl, Cn as a function of beta for three dpt values
% ---------------------------------------------------------------
for i=1:3;

   % input range
   u = [time zeros(21,1) betarad dpt(i)*ones(21,1)];

   % determine outputs
   [t,x,y]=sim('aeropropmodel',time,[],u);

   % plot results
   subplot(2,3,4); plot(beta,y(:,2),symbols{i});
      grid; hold on; xlabel('\beta [deg]','FontAngle','italic'); ylabel('C_Y','FontAngle','italic');
   subplot(2,3,5); plot(beta,y(:,4),symbols{i});
      grid; hold on; xlabel('\beta [deg]','FontAngle','italic'); ylabel('C_l','FontAngle','italic');
   subplot(2,3,6); plot(beta,y(:,6),symbols{i}); 
      grid; hold on; xlabel('\beta [deg]','FontAngle','italic'); ylabel('C_n','FontAngle','italic');

end

% Include legend in upper left-hand corner of upper left-hand plot
% ----------------------------------------------------------------
subplot(2,3,1);
legend('dpt = 0', 'dpt = 0.5', 'dpt = 1.0', 2);

% Include title
% -------------
subplot(2,3,2);
title('Force and moment coefficients as a function of \alpha, \beta, and dpt', ...
     'FontSize',15, ...
     'Position',[10 0.2], ...
     'VerticalAlignment','bottom');


%-----------------------------------------------------------------------
% The Flight Dynamics and Control Toolbox version 1.4.0. 
% (c) Copyright Marc Rauw, 1994-2005. Licensed under the Open Software 
% License version 2.1; see COPYING.TXT and LICENSE.TXT for more details.
% Last revision of this file: December 31, 2004. 
%-----------------------------------------------------------------------
