feat(Polar): add getcoeffs method

+af_tools/@Polar/Polar.m

@@ -146,6 +146,9 @@ classdef Polar
             self = findcllinrange(self);
         end
 
+        % Interpolate coefficients based on AoAs and Reynolds
+        [cl, cd, cm] = getcoeffs(self, aoa, reynolds)
+
     end
 
     % ----------------------------------------------------------------------------------------------

+af_tools/@Polar/getcoeffs.m

+function [cl, cd, cm] = getcoeffs(self, aoa, re, polarType)
+    % 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.
+    %
+    % Note
+    %   The inputs can be supplied as vectors of same size. In that case, the output will be the
+    %   coefficients for each couple of (aoa,reynolds)
+    %
+    % Note
+    %   This function uses 'spline' interpolation. This is because
+    %     - It is closer to the true curve anyway
+    %     - MATLAB automatically extrapolates to find values out of bounds, including for 2D interp.
+    % -----
+    % Usage
+    %   [cl, cd, cm] = Airfoil.getcoeffs(aoa, re)
+    %
+    % Input
+    %   aoa        : Angles of attacks, [rad]
+    %   re         : Reynolds numbers, [-]
+    %   polarType  : Type of Polar to use ('original', 'extended')
+    %
+    % Output
+    %   cl : Lift coefficients, [-]
+    %   cd : Drag coefficients, [-]
+    %   cm : Moment coefficients, [-]
+    %
+    % Examples
+    %   [cl, cd, cm] = Airfoil.getcoeffs(0.17, 1e6)
+    %   [cl, cd, cm] = Airfoil.getcoeffs([0.17,0.18], [1e5,1e6])
+    %
+    % See also: Airfoil.
+    % -----
+    % (c) Copyright 2022 University of Liege
+    % Author: Thomas Lambert <t.lambert@uliege.be>
+    % ULiege - Aeroelasticity and Experimental Aerodynamics
+    % MIT License
+    % https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox
+    % ----------------------------------------------------------------------------------------------
+
+    narginchk(3, 4);
+
+    % Defaults and constants
+    ALLOWED_TYPES = {'original', 'extended'};
+    DEF.POLAR_TYPE = 'extended';
+
+    % Input validation
+    if nargin == 3
+        polarType = DEF.POLAR_TYPE;
+    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');
+
+    % Interpolate coefficients
+    % -------------------------------------------
+
+    % Select correct polar
+    Pol = selectpolar(self, polarType);
+    plot(diff(Pol.reynolds));
+
+    if numel(self.reynolds) == 1
+        cl = checkandinterp1(Pol, 'cl', aoa);
+        cd = checkandinterp1(Pol, 'cd', aoa);
+        cm = checkandinterp1(Pol, 'cm', aoa);
+    else
+        % Create meshgrid with the polar data points
+        [X, Y] = meshgrid(Pol.aoa(:, 1), Pol.reynolds);
+
+        % Create meshgrid with the actual values for each element
+        [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);
+    end
+
+end
+
+function Polar = selectpolar(self, polarType)
+    % selectpolar Select the correct polar to use for the interpolation.
+
+    if strcmp(polarType, 'original')
+        tmpPol = self;
+    elseif strcmp(polarType, 'extended')
+        tmpPol = self.Ext;
+    end
+
+    Polar = struct('reynolds', self.reynolds, ...
+                   'aoa', tmpPol.aoa, ...
+                   'cl', tmpPol.cl, ...
+                   'cd', tmpPol.cd, ...
+                   'cm', tmpPol.cm);
+
+end
+
+function out = checkandinterp1(Pol, field, aoa)
+    % 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');
+    else
+        out = nan(size(aoa));
+    end
+
+end
+
+function out = checkandinterp2(X, Y, val, Xq, Yq)
+    % checkandinterp1 Checks if field is not empty and interpolate (1D)
+    %   If field is empty, return NaN vector of the size of AoA
+
+    if ~isempty(val)
+        out = diag(interp2(X, Y, val', Xq, Yq, 'spline'));
+    else
+        out = nan(length(Xq));
+    end
+
+end