Merge branch 'master' of gitlab.uliege.be:am-dept/matlab_airfoil_toolbox

+af_tools/+utils/parsefileinputs.m

@@ -36,6 +36,7 @@ function [filenames, filepaths, idxOpts] = parsefileinputs(optList, filetype, va
         ext = ['*', filetype];
         [filenames, filepaths] = uigetfile(ext, 'Select all dat-files to aggregate', ...
                                            'MultiSelect', 'on');
+        filenames = string(filenames);
         filepaths = repmat(string(filepaths), 1, length(filenames));
 
         idxOpts = 1;

+af_tools/+utils/spacedvector.m

 function vect = spacedvector(spacing, nPoints)
     % SPACEDVECTOR Returns a vector based on chosen spacing and number of points
+    %
+    % Spacing:
+    %   - linear     : points spaced linearily
+    %   - halfcosine : more dense at beginning and end of vector
+    %   - cosine     : more dense at the beginning and more spaced out at the end
+    %   - sine       : more dense at the end and more spaced out at the beginning
+    %
     % -----
-    % (c) Copyright 2022 University of Liege
+    % (c) Copyright 2022-2023 University of Liege
     % Author: Thomas Lambert <t.lambert@uliege.be>
     % ULiege - Aeroelasticity and Experimental Aerodynamics
     % Apache 2.0 License
@@ -18,6 +25,9 @@ function vect = spacedvector(spacing, nPoints)
         case 'cosine'
             beta = linspace(0, pi / 2, nPoints);
             vect = 1 - cos(beta);
+        case 'sine'
+            beta = linspace(0, pi / 2, nPoints);
+            vect = sin(beta);
     end
 
 end

+af_tools/@Polar/Polar.m

@@ -81,6 +81,7 @@ classdef Polar < handle
         Stall = struct('aoa', [], 'cl', [], 'cd', []) % Stall point data
         Zero  = struct('aoa', [], 'cd', []) % Zero-lift point data
         Lin   = struct('slope', [], 'aoaRange', [],  'clRange', []) % Lift linear range data
+        extrapMethod  (1, :) {mustBeMember(extrapMethod, {'none', 'spline', 'viterna'})} = 'spline'
     end
 
     % ----------------------------------------------------------------------------------------------
@@ -144,7 +145,7 @@ classdef Polar < handle
             %
             % See also: af_tools.findcllinearrange.
 
-            [self.Lin.clRange, self.Lin.aoaRange, self.Lin.slope] = ...
+            [self.Lin.aoaRange, self.Lin.clRange, self.Lin.slope] = ...
                 af_tools.findcllinearrange(self);
         end
 

+af_tools/@Polar/getcoeffs.m

-function [cl, cd, cm] = getcoeffs(self, aoa, re, polarType)
+function [cl, cd, cm] = getcoeffs(self, aoa, re, extrapMethod)
     % getcoeffs Interpolate polars and returns coefficients for given AoA and Reynolds.
     %   This function returns does mesh interpolation of the polars and returns the three
     %   coefficients (lift, drag, moment) for the input angles of attacks and reynolds.
+    %   It is possible to specify an extrapolation method as well if results outside the original
+    %   polars need to be retrieved. In that case, it is strongly suggested to use the 'Viterna'
+    %   method as it is the only one that tries to emulate the proper aerodynamic behavior of the
+    %   airfoil in the extreme ranges. If no extrapMethod is specified, self.extrapMethod will be
+    %   used (default is 'spline').
     %
     % Note
     %   The inputs can be supplied as vectors of same size. In that case, the output will be the
@@ -16,9 +21,9 @@ function [cl, cd, cm] = getcoeffs(self, aoa, re, polarType)
     %   [cl, cd, cm] = Airfoil.getcoeffs(aoa, re)
     %
     % Input
-    %   aoa        : Angles of attacks, [rad]
-    %   re         : Reynolds numbers, [-]
-    %   polarType  : Type of Polar to use ('original', 'extended')
+    %   aoa          : Angles of attacks, [rad]
+    %   re           : Reynolds numbers, [-]
+    %   extrapMethod : Extrapolation method ('none', 'spline', 'viterna')
     %
     % Output
     %   cl : Lift coefficients, [-]
@@ -45,28 +50,36 @@ function [cl, cd, cm] = getcoeffs(self, aoa, re, polarType)
     narginchk(3, 4);
 
     % Defaults and constants
-    ALLOWED_TYPES = {'original', 'extended'};
-    DEF.POLAR_TYPE = 'extended';
+    ALLOWED_EXTRAP = {'none', 'spline', 'viterna'};
+    EXTENDED_POLARS = {'viterna'};
 
     % Input validation
     if nargin == 3
-        polarType = DEF.POLAR_TYPE;
+        extrapMethod = self.extrapMethod;
     end
 
     validateattributes(aoa, {'double'}, {'row', 'nonempty'}, mfilename(), 'aoa');
     validateattributes(re, {'double'}, {'row', 'nonempty', 'size', size(aoa)}, mfilename(), 're');
-    polarType = validatestring(polarType, ALLOWED_TYPES, mfilename(), 'polarType');
+    extrapMethod = validatestring(extrapMethod, ALLOWED_EXTRAP, mfilename(), 'extrapMethod');
+
+    if any(strcmpi(extrapMethod, EXTENDED_POLARS))
+        polarType = 'extended';
+        if isempty(self.Ext.aoa)
+            self.extendpolar;
+        end
+    else
+        polarType = 'original';
+    end
 
     % Interpolate coefficients
     % -------------------------------------------
 
     % Select correct polar
     Pol = selectpolar(self, polarType);
-
     if numel(self.reynolds) == 1
-        cl = checkandinterp1(Pol, 'cl', aoa);
-        cd = checkandinterp1(Pol, 'cd', aoa);
-        cm = checkandinterp1(Pol, 'cm', aoa);
+        cl = checkandinterp1(Pol, 'cl', aoa, extrapMethod);
+        cd = checkandinterp1(Pol, 'cd', aoa, extrapMethod);
+        cm = checkandinterp1(Pol, 'cm', aoa, extrapMethod);
     else
         % Create meshgrid with the polar data points
         [X, Y] = meshgrid(Pol.aoa(:, 1), Pol.reynolds);
@@ -75,9 +88,9 @@ function [cl, cd, cm] = getcoeffs(self, aoa, re, polarType)
         [Xq, Yq] = meshgrid(aoa', re);
 
         % 2D interpolation
-        cl = checkandinterp2(X, Y, Pol.cl, Xq, Yq);
-        cd = checkandinterp2(X, Y, Pol.cd, Xq, Yq);
-        cm = checkandinterp2(X, Y, Pol.cm, Xq, Yq);
+        cl = checkandinterp2(X, Y, Pol.cl, Xq, Yq, extrapMethod);
+        cd = checkandinterp2(X, Y, Pol.cd, Xq, Yq, extrapMethod);
+        cm = checkandinterp2(X, Y, Pol.cm, Xq, Yq, extrapMethod);
     end
 
 end
@@ -99,28 +112,39 @@ function Polar = selectpolar(self, polarType)
 
 end
 
-function out = checkandinterp1(Pol, field, aoa)
+function out = checkandinterp1(Pol, field, aoa, extrapMethod)
     % checkandinterp1 Checks if field is not empty and interpolate (1D)
     %   If field is empty, return NaN vector of the size of AoA
 
     val = Pol.(field);
 
     if ~isempty(val)
-        out = interp1(Pol.aoa, val, aoa, 'spline');
+        if strcmpi(extrapMethod, 'none')
+            out = interp1(Pol.aoa, val, aoa, 'spline', NaN);
+        else
+            out = interp1(Pol.aoa, val, aoa, 'spline');
+        end
     else
         out = nan(size(aoa));
     end
 
 end
 
-function out = checkandinterp2(X, Y, val, Xq, Yq)
+function out = checkandinterp2(X, Y, val, Xq, Yq, extrapMethod)
     % checkandinterp1 Checks if field is not empty and interpolate (1D)
     %   If field is empty, return NaN vector of the size of AoA
 
+    warning('off', 'MATLAB:interp2:NaNstrip'); % Warning if original polars contain NaN already
+
     if ~isempty(val)
-        out = diag(interp2(X, Y, val', Xq, Yq, 'spline'));
+        if strcmpi(extrapMethod, 'none')
+            out = diag(interp2(X, Y, val', Xq, Yq, 'spline', NaN));
+        else
+            out = diag(interp2(X, Y, val', Xq, Yq, 'spline'));
+        end
     else
         out = nan(length(Xq));
     end
 
+    warning('on', 'MATLAB:interp2:NaNstrip');
 end

+af_tools/formatairfoilcoord.m

@@ -139,7 +139,7 @@ function Dat = formatairfoilcoord(varargin)
         end
 
         % Output structure
-        Dat(i).path = fullpaths;
+        Dat(i).path = fullpaths{i};
         Dat(i).file = allFileNames{i};
         Dat(i).airfoil = char(cellstr(tmpAirfoil{:}));
         Dat(i).format = outputFormat;

+af_tools/nacacamber.m

@@ -6,13 +6,13 @@ function [coord, gradY] = nacacamber(varargin)
     % -----
     %
     % Usage:
-    %   [X, Y] = NACACAMBER prompts the user for the airfoil digits, then calculates the coordinates
-    %   of the camberline using 100 points spaced with the half-cosine rule.
+    %   [coord, gradY] = NACACAMBER prompts the user for the airfoil digits, then calculates the
+    %   coordinates of the camberline using 100 points spaced with the half-cosine rule.
     %
-    %   [X, Y] = NACACAMBER(DIGITS) calculates the coordinates of the camberline for the NACA
-    %   airfoil specified by DIGITS, without further input.
+    %   [coord, gradY] = NACACAMBER(DIGITS) calculates the coordinates of the camberline for the
+    %   NACA airfoil specified by DIGITS, without further input.
     %
-    %   [X, Y] = NACACAMBER(DIGITS, NPOINTS) calculates the coordinates of the NACA airfoil
+    %   [coord, gradY] = NACACAMBER(DIGITS, NPOINTS) calculates the coordinates of the NACA airfoil
     %   specified by DIGITS, using NPOINTS points spaced with the half-cosine rule.
     %
     %   [...] = NACACAMBER(..., 'spacing', SP) uses the defined spacing rule to determine to
@@ -20,14 +20,14 @@ function [coord, gradY] = nacacamber(varargin)
     %     'linear'     - linear
     %     'halfcosine' - (default) finer at leading and trailing edges
     %     'cosine'     - finer at leading edge, coarser at the trailing edge
+    %     'sine'       - finer at trailing edge, coarser at the leading edge
     %
     % Inputs:
     %   digits  : Character vector representing the NACA airfoil (ex: '0012', '24120', '24012')
     %   nPoints : Number of output coordinates
     %
     % Output:
-    %   xc : Camberline X coordinates (0 <= X <= 1)
-    %   yc : Camberline Y coordinates
+    %   coord : Camberline coordinates [nPts x 2]
     %   gradY : Gradient of Y (required to get the whole airfoil coordinates)
     %
     % Example:
@@ -41,7 +41,7 @@ function [coord, gradY] = nacacamber(varargin)
     % <a href="https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox">Documentation (online)</a>
 
     % -------------------------------------------
-    % (c) Copyright 2022 University of Liege
+    % (c) Copyright 2022-2023 University of Liege
     % Author: Thomas Lambert <t.lambert@uliege.be>
     % ULiege - Aeroelasticity and Experimental Aerodynamics
     % MIT License