diff --git a/+af_tools/@Polar/Polar.m b/+af_tools/@Polar/Polar.m
index ecedbd173811c1f3c0eb25ad8271c40baa501e51..0ad9dd6cfe488e239c498f544b36ca64b7bce80b 100644
--- a/+af_tools/@Polar/Polar.m
+++ b/+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
 
     % ----------------------------------------------------------------------------------------------
diff --git a/+af_tools/@Polar/getcoeffs.m b/+af_tools/@Polar/getcoeffs.m
index e4da05fb7894d1174e8787ef74ec1115cea19435..88096a2ade63ffedacae6d12bf667558ebb633c4 100644
--- a/+af_tools/@Polar/getcoeffs.m
+++ b/+af_tools/@Polar/getcoeffs.m
@@ -1,7 +1,12 @@
-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