feat(Polar): add interpolated polars

This allows the user to create polars for interpolated airfoils.
When the shape of the local airfoil is unknown but is between sections
of known airfoils, this allows to create interpolated airfoils polars
based on the known polars of the bounding airfoils.

+af_tools/@Polar/Polar.m

@@ -180,6 +180,8 @@ classdef Polar < handle
         self = loadpolar(self, polarFile) % Load polar from file
         self = polypolar(polyCl, polyCd) % Create polar using polynomial coefficients
 
+        % Create a polar by interpolating between two different polar objects
+        self = interppolar(Polar1, Polar2, position)
     end
 
 end

+af_tools/@Polar/interppolar.m

+function self = interppolar(Polar1, Polar2, weightPolar2)
+    % INTERPPOLAR  Create a new polar object by interpolation of two other polars.
+    %   This static method can be used in place of a normal constructor.
+    %   This method can be used when the user wants to interpolate the polars bewteen two known
+    %   airfoils. Consider that the airfoil at x=0 is a NACA 0012, and the one at x=1 is a CLARY-Y.
+    %   Rather than assuming that all segments in [0,1] are NACA 0012 and all segments after 1 are
+    %   CLARK-Y, this function allows to create polars for the [0,1] interval that are a mix between
+    %   the two bounds. If we are at 0.1, the newly created Polar object, will then be 90% of NACA
+    %   0012 and 10% of CLARK-Y.
+    % Note:
+    %   The two input Polars should have the same range of AOA, Reynolds, Mach, and nCrit.
+    % -----
+    %
+    % Usage:
+    %   MiedPolar = Polar.INTERPPOLAR(Polar1, Polar2, weightPolar2) create an new Polar object that
+    %   is a linear interpolation between Polar1 and Polar2, with a given weight for Polar2.
+    %
+    % Inputs:
+    %   Polar1 : Polar object for one part of the interpolation.
+    %   Polar2 : Polar object for the second part of the interpolation. Muse have the same Re, Mach
+    %            and nCrit values as Polar1.
+    %   weightPolar2: Scalar in [0-1] to skew the interpolation  properly between Polar1 and Polar2.
+    %
+    % Output:
+    %   self : Polar object that results from the interpolation.
+    %
+    % Example:
+    %   MixedPolar = af_tools.interppolar(PolarNACA0012, PolarClarkY, 0.6); Creates an interpolated
+    %                Polar object MixedPolar that results from the interpolation of PolarNACA0012
+    %                and PolarClarkY, asusming that the MixedPolar is made of 60% ClarkY and 40%
+    %                NACA0012.
+    %
+    % See also: af_tools.Polar.
+    %
+    % -----
+    % <a href="https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox">Documentation (online)</a>
+
+    % ----------------------------------------------------------------------------------------------
+    % (c) Copyright 2022-2024 University of Liege
+    % Author: Thomas Lambert <t.lambert@uliege.be>
+    % ULiege - Aeroelasticity and Experimental Aerodynamics
+    % MIT License
+    % Repo: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox
+    % Issues: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/issues
+    % ----------------------------------------------------------------------------------------------
+
+    % --- Defaults and constants
+    DEBUG = true;
+
+    % --- Input checks
+    % Check if both Polars are Polars objects and weightPolar2 is in [0,1]
+    validateattributes(Polar1, {'af_tools.Polar'}, {'scalar'}, mfilename(), 'Polar1', 1);
+    validateattributes(Polar2, {'af_tools.Polar'}, {'scalar'}, mfilename(), 'Polar2', 2);
+    validateattributes(weightPolar2, {'double'}, {'scalar', 'nonnegative', '<=', 1}, ...
+                       mfilename(), 'weightPolar2', 3);
+
+    % Check that both Polars have the same AOA, Reynolds, Mach and nCrit
+    equalfieldscheck(Polar1, Polar2, 'aoa');
+    equalfieldscheck(Polar1, Polar2, 'reynolds');
+    equalfieldscheck(Polar1, Polar2, 'mach');
+    equalfieldscheck(Polar1, Polar2, 'nCrit');
+
+    % --- Interpolate the polars
+
+    self = af_tools.Polar;
+
+    wPol2 = weightPolar2;
+    wPol1 = 1 - wPol2;
+    nPolars = numel(Polar1.reynolds);
+
+    % Set base polar information
+    afName = {sprintf('MIX-%d%%%s-%d%%%s', round(wPol1 * 100), Polar1.airfoil{1}, ...
+                      round(wPol2 * 100), Polar2.airfoil{2})};
+    self.airfoil = repmat(afName, 1, nPolars);
+
+    self.reynolds = Polar1.reynolds;
+    self.mach = Polar1.mach;
+    self.nCrit = Polar1.nCrit;
+
+    % Interpolate the angle of attack, cl, cd and cm using the base polars
+    if Polar1.aoa ~= Polar2.aoa
+        % Reinterpolate the original polars so they share the same range of AOA.
+        aoaMin = min(min(Polar1.aoa(:, 1)), min(Polar2.aoa(:, 1)));
+        aoaMax = max(max(Polar1.aoa(:, 1)), max(Polar2.aoa(:, 1)));
+        aoaStep = min(diff(Polar1.aoa(1:2, 1)), diff(Polar2.aoa(1:2, 1)));
+        aoaVect = aoaMin:aoaStep:aoaMax;
+
+        aoa1 = repmapt(aoaVect', 1, nPolars);
+        aoa2 = repmapt(aoaVect', 1, nPolars);
+        cl1 = interp1(Polar1.aoa(:, 1), Polar1.cl, aoaVect');
+        cl2 = interp1(Polar2.aoa(:, 1), Polar2.cl, aoaVect');
+        cd1 = interp1(Polar1.aoa(:, 1), Polar1.cd, aoaVect');
+        cd2 = interp1(Polar2.aoa(:, 1), Polar2.cd, aoaVect');
+        cm1 = interp1(Polar1.aoa(:, 1), Polar1.cm, aoaVect');
+        cm2 = interp1(Polar2.aoa(:, 1), Polar2.cm, aoaVect');
+    else
+        % Use direclty the polar data
+        aoa1 = Polar1.aoa;
+        aoa2 = Polar2.aoa;
+        cl1 = Polar1.cl;
+        cl2 = Polar2.cl;
+        cd1 = Polar1.cd;
+        cd2 = Polar2.cd;
+        cm1 = Polar1.cm;
+        cm2 = Polar2.cm;
+    end
+    self.aoa = lininterparray(wPol1, aoa1, wPol2, aoa2);
+    self.cl = lininterparray(wPol1, cl1, wPol2, cl2);
+    self.cd = lininterparray(wPol1, cd1, wPol2, cd2);
+    self.cm = lininterparray(wPol1, cm1, wPol2, cm2);
+
+    % --- DEBUGGING
+    if DEBUG
+        % Plot the original polars for each Reynolds as well as the interpolated one
+        for j = 1:nPolars
+            figure('Name', 'Combined Polars');
+            hold on;
+            plot(Polar1.aoa(:, j), Polar1.cl(:, j));
+            plot(Polar2.aoa(:, j), Polar2.cl(:, j));
+            plot(self.aoa(:, j), self.cl(:, j));
+            legend(Polar1.airfoil{1}, Polar2.airfoil{1}, self.airfoil{1}, 'Location', 'NorthWest');
+            grid on;
+        end
+    end
+
+end
+
+function interpArray = lininterparray(weight1, array1, weight2, array2)
+    % LININTERPARRAY Linear interpolation between two arrays
+
+    interpArray = weight1 * array1 + weight2 * array2;
+
+end
+
+function equalfieldscheck(Polar1, Polar2, field)
+    % FIELDCHECK Check if some fields are equal in both Polar objects
+
+    if Polar1.(field) ~= Polar2.(field)
+        error('Polar:interppolar:different%s', ...
+              ['Impossible to interpolate between the two Polar objects as they have '...
+               'different %s.\n'...
+               'Please use two Polar objects with identical AOA, Reynolds, Mach and nCrit'], ...
+              field, field);
+    end
+end