From 40a17792e850eef2febb096a042069c58309a000 Mon Sep 17 00:00:00 2001
From: Thomas Lambert <t.lambert@uliege.be>
Date: Tue, 2 Aug 2022 21:59:19 +0200
Subject: [PATCH] feat(Polar): add getcoeffs method
---
+af_tools/@Polar/Polar.m | 3 +
+af_tools/@Polar/getcoeffs.m | 123 +++++++++++++++++++++++++++++++++++
2 files changed, 126 insertions(+)
create mode 100644 +af_tools/@Polar/getcoeffs.m
diff --git a/+af_tools/@Polar/Polar.m b/+af_tools/@Polar/Polar.m
index 2dd433e..74f0bf8 100644
--- a/+af_tools/@Polar/Polar.m
+++ b/+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
% ----------------------------------------------------------------------------------------------
diff --git a/+af_tools/@Polar/getcoeffs.m b/+af_tools/@Polar/getcoeffs.m
new file mode 100644
index 0000000..a1d774e
--- /dev/null
+++ b/+af_tools/@Polar/getcoeffs.m
@@ -0,0 +1,123 @@
+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
--
GitLab