diff --git a/+af_tools/@Polar/Polar.m b/+af_tools/@Polar/Polar.m
index 196137adcfac9ea8e661ebeb470116429c92dbd8..69fbc62174ccb2da4668d25048bc98d86d3b849d 100644
--- a/+af_tools/@Polar/Polar.m
+++ b/+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
diff --git a/+af_tools/@Polar/interppolar.m b/+af_tools/@Polar/interppolar.m
new file mode 100644
index 0000000000000000000000000000000000000000..054ff051bdb04b7375390b3ed01511f95281be13
--- /dev/null
+++ b/+af_tools/@Polar/interppolar.m
@@ -0,0 +1,145 @@
+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