Compare revisions

Thomas Lambert · Thomas Lambert · Thomas Lambert · Thomas Lambert · Thomas Lambert · Thomas Lambert
--- a/.gitignore
+++ b/.gitignore
@@ -7,6 +7,7 @@
 # Folders to ignore
 results/
+validation/private
 # Ignore all configs by default and only add some to be kept
 src/configs/*

--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -27,7 +27,7 @@ style_check:
  <<: *miss_hit
  stage: static analysis
  script:
-    - mh_style . --fix
+    - mh_style .
 metric_check:
  <<: *miss_hit

--- a/.gitmodules
+++ b/.gitmodules
@@ -4,3 +4,6 @@
 [submodule "src/libs/matlab_airfoil_toolbox"]
 	path = src/libs/matlab_airfoil_toolbox
 	url = ../../am-dept/matlab_airfoil_toolbox
+[submodule "src/libs/octave-atmosisa"]
+	path = src/libs/octave-atmosisa
+	url = https://github.com/lentzi90/octave-atmosisa
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,10 +9,6 @@ to [Semantic Versioning][sem_ver].
 ### Added
- Basic result summary table
- Basic result plotting function
- Result filtering function
 ### Changed
 ### Deprecated
@@ -21,6 +17,59 @@ to [Semantic Versioning][sem_ver].
 ### Fixed
+## [0.2.2] - 2024-03-13
+### Changed
+- Allow user to decide if airfoils should be interpolated between reference
+  sections (slow) or if same airfoil should be used until next section (fast)
+## [0.2.1] - 2024-03-13
+### Fixed
+- Issue with Polar interpolation between different airfoils
+## [0.2.0] - 2024-03-13
+### Added
+- Proper interpolation of Airfoils and Polars between reference sections
+## [0.1.6] - 2024-02-12
+### Fixed
+- Issue with toolbox check
+## [0.1.5] - 2023-12-14
+### Changed
+- Add `octave-atmosisa` in libs so we can use it if _Aerospace Toolbox_ is not
+  installed
+- Update `matlab_airfoil_toolbox` lib to remove dependency to the _Signal
+  Processing Toolbox_ in `findstall`.
+## [0.1.4] - 2023-12-14
+### Added
+- Result summary table
+- Result plotting function
+- Result filtering function
+- Preliminary support for coaxial rotors
+### Changed
+- Implement rotation matrices to remove useless dependency
+- Defaults values in polar generator script
+### Fixed
+- Typo in Caradonna config
+- Remove `--fix` flag for miss_hit in the pipeline
 ## [0.1.3] - 2023-05-25
 ### Added
@@ -103,7 +152,12 @@ _Initial commit_: single rotor in hover/axial flows
 [sem_ver]:<https://semver.org/spec/v2.0.0.html>
 [keep_chglog]: <https://keepachangelog.com/en/1.0.0/>
-[Unreleased]: https://gitlab.uliege.be/rotare/rotare/compare/0.1.3...main
+[Unreleased]: https://gitlab.uliege.be/rotare/rotare/compare/0.2.2...main
+[0.2.2]: https://gitlab.uliege.be/rotare/rotare/compare/0.2.1...0.2.2
+[0.2.1]: https://gitlab.uliege.be/rotare/rotare/compare/0.2.0...0.2.1
+[0.1.6]: https://gitlab.uliege.be/rotare/rotare/compare/0.1.5...0.1.6
+[0.1.5]: https://gitlab.uliege.be/rotare/rotare/compare/0.1.4...0.1.5
+[0.1.4]: https://gitlab.uliege.be/rotare/rotare/compare/0.1.3...0.1.4
 [0.1.3]: https://gitlab.uliege.be/rotare/rotare/compare/0.1.2...0.1.3
 [0.1.2]: https://gitlab.uliege.be/rotare/rotare/compare/0.1.1...0.1.2
 [0.1.1]: https://gitlab.uliege.be/rotare/rotare/compare/0.1.0...0.1.1

--- a/README.md
+++ b/README.md
@@ -35,6 +35,7 @@ A more exhaustive list of features can be found in the [complete
 documentation][rotare-doc] that will be made available with version 1.0.0.
 ### Types of rotors
 - [x] Helicopters main or tail rotors
 - [x] Aircraft propellers
 - [ ] Wind/tidal turbines
@@ -91,6 +92,8 @@ name `COMPLETE CODE`. This archive contains all the dependencies needed for the
 proper use of Rotare. Always download this one and not the automatically
 generated _source code_ archives.
+In theory, no additional toolboxes should be required to work with Rotare.
 ## Documentation
 Rotare comes with a [detailed documentation][rotare-doc]. This documentation is
@@ -120,7 +123,7 @@ troubleshooting section of the documentation.
 If you encounter any other issue with this code, please check [the issue
 tracker][rotare-issues] and fill a new issue report if applicable. You can also
-contact me directly at tlambert@uliege.be.
+contact me directly at <tlambert@uliege.be>.
 ## Other

--- a/RELEASE.md
+++ b/RELEASE.md
@@ -10,4 +10,5 @@ automatically generated and do not contain all dependencies.
 This release features:
- Result output and analysis methods
+- Allow user to decide if airfoils should be interpolated between reference
+  sections (slow) or if same airfoil should be used until next section (fast)
--- a/miss_hit.cfg
+++ b/miss_hit.cfg
@@ -29,6 +29,6 @@ regex_method_name: "[a-z](_?[a-zA-Z0-9]+)*"
 metric "cnest": limit 5
 metric "file_length": limit 800
 metric "cyc": limit 15
-metric "parameters": limit 10
+metric "parameters": limit 11
 metric "globals": limit 0
 metric "persistent": limit 3
--- a/src/.miss_hit
+++ b/src/.miss_hit
+# miss_hi style checker
+# (https://florianschanda.github.io/miss_hit/style_checker.html)
+#
+# Exclude libraries directory from style checking
+exclude_dir: "libs"
--- a/src/airfoil_data/polargenerator.py
+++ b/src/airfoil_data/polargenerator.py
+#!/usr/bin/python3
 """Generate airfoil polars automatically using XFOIL
 Rotare requires the input of airfoil polars to interpolate properly the lift
@@ -27,25 +28,30 @@ import os
 import numpy as np
 import aeropy.xfoil_module as xf
+# Defaults for polar generation
+# AOA
+REYNOLDS = [re * 1000000 for re in [0.02, 0.1, 0.2, 0.5, 1, 2, 5, 10]]
+AOA = list(np.arange(-16, 22, 0.5))
+MACH = 0.1
 def main():
    """Creates polar files for each reynolds"""
    airfoil = parse_inputs()
-    reynolds = [re * 1000000 for re in [0.2, 0.5, 1, 1.5, 2, 5, 10]]
-    aoa = list(np.arange(-18, 25, 0.5))
    i = 1
-    for re in reynolds:
+    for re in REYNOLDS:
-        print("%d/%d - Calculating polar for Re=%ge6" % (i, len(reynolds), re / 1e6))
+        print("%d/%d - Calculating polar for Re=%ge6" % (i, len(REYNOLDS), re / 1e6))
        i = i + 1
        if airfoil.startswith("naca"):
            xf.call(
                airfoil,
-                alfas=aoa,
+                alfas=AOA,
                output="Polar",
                Reynolds=re,
+                Mach=MACH,
                plots=False,
                NACA=True,
                iteration=500,
@@ -54,9 +60,10 @@ def main():
            if os.path.exists(airfoil) or os.path.exists(airfoil + ".dat"):
                xf.call(
                    airfoil,
-                    alfas=aoa,
+                    alfas=AOA,
                    output="Polar",
                    Reynolds=re,
+                    Mach=MACH,
                    plots=False,
                    NACA=False,
                    iteration=500,
@@ -86,7 +93,7 @@ def cleanup(airfoil):
    for file in os.listdir(cwd):
        if file.startswith("Coordinates_"):
            safe_rename(file, airfoil + ".dat")
-        elif file.startswith("Polar_"):
+        elif file.startswith("Polar_") and not (file.endswith(".txt")):
            safe_rename(file, file + ".txt")

--- a/src/classes/@Blade/Blade.m
+++ b/src/classes/@Blade/Blade.m
@@ -19,31 +19,35 @@ classdef Blade < handle
    % -----
    %
    % Blade properties:
-    %   nElem - Number of elements, [-]
+    %   nElem    - Number of elements, [-]
-    %   dy    - Element span, [m]
+    %   dy       - Element span, [m]
-    %   y     - Element absolute radial position, [m]
+    %   y        - Element absolute radial position, [m]
-    %   r     - Element relative radial position (relative to the total radius), [-]
+    %   r        - Element relative radial position (relative to the total radius), [-]
-    %   area  - Element area, [m^2]
+    %   area     - Element area, [m^2]
-    %   chord - Element chord, [m]
+    %   chord    - Element chord, [m]
-    %   twist - Element twist (also called stagger angle), [rad]
+    %   twist    - Element twist (also called stagger angle), [rad]
-    %   iAf   - Index of airfoil to use for each element, [-]
+    %   iAf      - Index of airfoil to use for each element, [-]
-    %   sol   - Element solidity, [-]
+    %   sol      - Element solidity, [-]
+    %   afInterp - Interpolate airfoil between base sections, [true/false]
    %
    % Blade methods:
+    %   optimizeairfoils - Optimize airfoils to avoid unecessary calls for interpolated airfoils
    %
    % Blade constructor:
    %   Bl = Blade() creates an empty object.
    %
-    %   Bl = Blade(nBlades, rad, chord, twist, iAf, nElem) creates a Blade object based on the input
+    %   Bl = Blade(nBlades, rad, chord, twist, iAf, nElem, afInterp) creates a Blade object based on
-    %   provided (see list below).
+    %   the input provided (see list below).
    %
    % Constructor inputs:
-    %   nBlades : Number of blades, [-]
+    %   nBlades  : Number of blades, [-]
-    %   rad     : Radius of the guide stations, [m] (at least 2 elements (root and tip))
+    %   rad      : Radius of the guide stations, [m] (at least 2 elements (root and tip))
-    %   chord   : Chord of the guide stations, [m] (same size as rad)
+    %   chord    : Chord of the guide stations, [m] (same size as rad)
-    %   twist   : Twist of the guide stations, [deg] (same size as rad)
+    %   twist    : Twist of the guide stations, [deg] (same size as rad)
-    %   iAf     : Index of the airfoil for the guide stations, [-] (same size as rad)
+    %   iAf      : Index of the airfoil for the guide stations, [-] (same size as rad)
-    %   nElem   : Number of elements to mesh the whole blade, [-]
+    %   nElem    : Number of elements to mesh the whole blade, [-]
+    %   afInterp : Interpolate airfoil between base sections, [true/false]
+    %
    %
    % See also: Rotor, rotare, template, af_tools.Airfoil.
    %
@@ -69,6 +73,7 @@ classdef Blade < handle
        chord (1, :) double {mustBePositive} % Element chord, [m]
        twist (1, :) double {mustBeFinite}   % Element twist (also called stagger angle), [rad]
        iAf (1, :) double {mustBePositive}   % Index of airfoil to use for each element, [-]
+        afInterp (1, 1) logical = false      % Interpolate airfoil between base sections
    end
    properties (SetAccess = private, Dependent)
@@ -82,7 +87,7 @@ classdef Blade < handle
    methods
-        function self = Blade(nBlades, rad, chord, twist, iAf, nElem)
+        function self = Blade(nBlades, rad, chord, twist, iAf, nElem, afInterp)
            % BLADE Constructor.
            %   Discretizes the blade in elements, using a spline interpolation between the base
            %   stations passed as input. See main class help for details.
@@ -94,6 +99,7 @@ classdef Blade < handle
                self.nElem = nElem;
                self.nBlades = nBlades;
+                self.afInterp = afInterp;
                % Determine element radial position
                self.dy = (rad(end) - rad(1)) / nElem;
@@ -105,8 +111,12 @@ classdef Blade < handle
                self.chord = interp1(rad, chord, self.y, INTERP_METHOD);
                % Assign airfoil to each element
-                for i = 1:length(rad) - 1
+                if self.afInterp
-                    self.iAf(self.y >= rad(i) & self.y <= rad(i + 1)) = iAf(i);
+                    self.iAf = 1:1:nElem;
+                else
+                    for i = 1:length(rad) - 1
+                        self.iAf(self.y >= rad(i) & self.y <= rad(i + 1)) = iAf(i);
+                    end
                end
                % Element area (generic elements considered trapezoids)
@@ -122,5 +132,13 @@ classdef Blade < handle
            sol = self.nBlades .* self.chord ./ (2 * pi * self.y);
        end
+        % ---------------------------------------------
+        % Other methods
+        function self = optimizeairfoils(self, newiAf)
+            % OPTIMIZEAIRFOILS Optimize airfoils to avoid unecessary calls for interpolated airfoils
+            self.iAf = newiAf;
+        end
    end
 end
--- a/src/classes/@ElemPerf/ElemPerf.m
+++ b/src/classes/@ElemPerf/ElemPerf.m
@@ -56,7 +56,8 @@ classdef ElemPerf < handle
        Rot (1, 1) Rotor % Handle of the rotor linked to this specific ElemPerf instance
        Op (1, 1)  Oper  % Handle of the Operating point linked to this specific ElemPerf instance
-        reynolds (1, :) double % Reynolds number, [-]
+        upstreamVelAx (1, :) double % Upstream axial velocity above rotor disk, [m/s]
+        upstreamVelTg (1, :) double % Upstream tangential velocity above rotor disk, [m/s]
        tgSpeed   (1, :) double % Tangential speed, [m/s]
        pitch     (1, :) double % Pitch (twist + collective), [rad]
@@ -83,6 +84,7 @@ classdef ElemPerf < handle
    % from a single set method otherwise.
    properties (Dependent)
        alpha    (1, :) double % Angle of attack, [rad]
+        reynolds (1, :) double % Reynolds number, [-]
        indVelAx (1, :) double % Induced axial velocity at rotor disk, [m/s]
        indVelTg (1, :) double % Induced tangential velocity at rotor disk, [m/s]
        inflowRat (1, :) double % Inflow ratio, [-]
@@ -94,6 +96,7 @@ classdef ElemPerf < handle
    % Cache for dependent properties to avoid excessive recalculation
    properties (Access = private, Hidden)
        alpha_    (1, :) double % Angle of attack, [rad]
+        reynolds_ (1, :) double % Reynolds number, [-]
        indVelAx_ (1, :) double % Induced axial velocity at rotor disk, [m/s]
        indVelTg_ (1, :) double % Induced tangential velocity at rotor disk, [m/s]
        inflowRat_ (1, :) double % Inflow ratio, [-]
@@ -117,13 +120,15 @@ classdef ElemPerf < handle
                % Blade pitch
                self.pitch = Rot.Bl.twist + Op.coll;
+                % Upstream axial velocity
+                self.upstreamVelAx = ones(size(self.pitch)) * Op.speed;
+                % Upstream tangential velocity
+                self.upstreamVelTg = zeros(size(self.pitch));
                % In-plane velocities
                self.tgSpeed = Op.omega .* Rot.Bl.y;
-                % Reynolds
-                relVel = sqrt(self.tgSpeed.^2 + Op.speed.^2);
-                self.reynolds = Flow.reynolds(relVel, Rot.Bl.chord, Op.Flow.mu, Op.Flow.rho);
                % Get alpha_0 from reynolds
                self.alpha0 = zeros(size(self.pitch));
                for i = 1:length(Rot.Af)
@@ -154,6 +159,23 @@ classdef ElemPerf < handle
            alpha = self.alpha_;
        end
+        function reynolds = get.reynolds(self)
+            % Always recalculate Reynolds based on the most accurate data
+            velAx = self.upstreamVelAx;
+            velTg = self.tgSpeed;
+            if ~isempty(self.indVelAx)
+                velAx = velAx + self.indVelAx;
+            end
+            if ~isempty(self.indVelTg)
+                velTg =  velTg - self.indVelTg;
+            end
+            relVel = sqrt(velAx.^2 + velTg.^2);
+            reynolds = Flow.reynolds(relVel, self.Rot.Bl.chord, ...
+                                     self.Op.Flow.mu, self.Op.Flow.rho);
+        end
        function indVelAx = get.indVelAx(self)
            indVelAx = self.indVelAx_;
        end
@@ -213,7 +235,7 @@ classdef ElemPerf < handle
            % Calculate induced velocity and ratio
            self.indVelAx_ = val .* self.Op.omega * self.Rot.radius - self.Op.speed;
-            self.indInflowRat_ = self.indVelAx_ ./ (self.Op.omega * self.Rot.radius);
+            self.indInflowRat_ = self.indVelAx ./ (self.Op.omega * self.Rot.radius);
        end
        function set.swirlRat(self, val)
@@ -222,7 +244,7 @@ classdef ElemPerf < handle
            % Calculate induced swirl velocity and ratio
            self.indVelTg_ = self.tgSpeed - val .* (self.Op.omega * self.Rot.radius);
-            self.indSwirlRat_ = self.tgSpeed ./ (self.Op.omega * self.Rot.radius) - val;
+            self.indSwirlRat_ = self.indVelTg ./ (self.Op.omega * self.Rot.radius);
        end
        function set.indInflowRat(self, val)

--- a/src/classes/@ElemPerf/calcforces.m
+++ b/src/classes/@ElemPerf/calcforces.m
@@ -38,7 +38,7 @@ function calcforces(self)
    if ~isempty(self.cl)
-        axVel = self.Op.speed + self.indVelAx;
+        axVel = self.upstreamVelAx + self.indVelAx;
        tgVel = self.tgSpeed - self.indVelTg;
        relVel = sqrt(axVel.^2 + tgVel.^2);

--- a/src/classes/@ElemPerf/plotveltriangles.m
+++ b/src/classes/@ElemPerf/plotveltriangles.m
@@ -67,13 +67,13 @@ function plotveltriangles(self, nTriangles, varargin)
    % Velocity triangles
    % FIXME: Needs adaptation for coaxial
-    vAx_up = ones(size(self.tgSpeed)) * self.Op.speed;
+    vAx_up = ones(size(self.tgSpeed)) .* self.upstreamVelAx;
    vTg_up = self.tgSpeed;
-    vAx_rot = self.Op.speed + self.indVelAx;
+    vAx_rot = self.upstreamVelAx + self.indVelAx;
    vTg_rot = self.tgSpeed - self.indVelTg;
-    vAx_down = self.Op.speed + 2 * self.indVelAx;
+    vAx_down = self.upstreamVelAx + 2 * self.indVelAx;
    vTg_down = self.tgSpeed - 2 * self.indVelTg;
    % Determine sections to be plotted

--- a/src/classes/@OperRotor/OperRotor.m
+++ b/src/classes/@OperRotor/OperRotor.m
@@ -100,7 +100,7 @@ classdef OperRotor < handle
        end
        function relSpeed = get.relTipSpeed(self)
-            relSpeed = sqrt(self.tgTipSpeed^2 + self.Op.speed^2);
+            relSpeed = sqrt(self.tgTipSpeed.^2 + self.upstreamVelAx.^2);
        end
        function advRat = get.advanceRatio(self)

--- a/src/classes/@Result/plotperf.m
+++ b/src/classes/@Result/plotperf.m
@@ -45,11 +45,16 @@ function plotperf(self, varargin)
    % Defaults and constants
    DEF.ALLOWED_PROPERTIES = {'thrust', 'torque', 'power', 'eff', 'cT', 'cQ', 'cP'};
-    DEF.allowed_oper = self.operPts.Properties.VariableNames;
+    DEF.base_oper = self.operPts.Properties.VariableNames;
+    DEF.allowed_oper = [DEF.base_oper, 'advRat']; % Also allows advance ratio
    DEF.ALLOWED_SOLVERS = {'leishman', 'indfact', 'indvel', 'stahlhut'};
+    DEF.LEISHMAN_LINE = {'LineStyle', '-', 'color', 'red'};
+    DEF.INDFACT_LINE = {'LineStyle', ':', 'color', 'magenta'};
+    DEF.INDVEL_LINE = {'LineStyle', '--', 'color', 'blue'};
+    DEF.STAHLHUT_LINE = {'LineStyle', '-.', 'color', 'green'};
    % Parse inputs
-    [perf, oper, solvers, operPt, newFig] = parseinputs(DEF, varargin{:});
+    [perf, oper, solvers, operPt, newFig, lineSpec] = parseinputs(DEF, varargin{:});
    % Plot the figure
    if newFig
@@ -60,18 +65,22 @@ function plotperf(self, varargin)
    for iSolv = 1:length(solvers)
        thissolv = solvers{iSolv};
        if ~isempty(self.(thissolv))
+            if strcmpi(oper, 'advRat')
+                operVect = vectorize(filtered.(thissolv), 'advanceRatio');
+            else
+                operVect = self.operPts(idx, :).(oper);
+            end
            dataVect = vectorize(filtered.(thissolv), perf);
-            plot(self.operPts(idx, :).(oper), dataVect, 'DisplayName', solvers{iSolv});
+            plot(operVect, dataVect, lineSpec.(thissolv){:}, 'DisplayName', thissolv);
            hold on;
        end
    end
 end
-function [perf, oper, solvers, operPt, newFig] = parseinputs(DEF, varargin)
+function [perf, oper, solvers, operPt, newFig, lineSpec] = parseinputs(DEF, varargin)
    % PARSEINPUTS Parse the varargin inputs
    %   - Only the parameter corresponding to oper can be nan, all others need to be defined
@@ -87,18 +96,33 @@ function [perf, oper, solvers, operPt, newFig] = parseinputs(DEF, varargin)
    addRequired(p, 'collective', valData);
    addParameter(p, 'solvers', DEF.ALLOWED_SOLVERS);
    addParameter(p, 'newFig', false, valLogi);
+    addParameter(p, 'leishmanLine', DEF.LEISHMAN_LINE);
+    addParameter(p, 'indfactLine', DEF.INDFACT_LINE);
+    addParameter(p, 'indvelLine', DEF.INDVEL_LINE);
+    addParameter(p, 'stahlhutLine', DEF.STAHLHUT_LINE);
    parse(p, varargin{:});
    perf = p.Results.perf;
    oper = p.Results.oper;
    solvers = p.Results.solvers;
    operPt = struct('altitude', [], 'speed', [], 'rpm', [], 'collective', []);
-    for i = 1:length(DEF.allowed_oper)
+    for i = 1:length(DEF.base_oper)
-        curOp = DEF.allowed_oper{i};
+        curOp = DEF.base_oper{i};
        operPt.(curOp) = p.Results.(curOp);
    end
-    operPt.(p.Results.oper) = NaN;
+    if strcmpi(p.Results.oper, 'advRat')
+        operPt.rpm = NaN;
+        operPt.speed = NaN;
+    else
+        operPt.(p.Results.oper) = NaN;
+    end
    newFig = p.Results.newFig;
+    lineSpec.leishman = p.Results.leishmanLine;
+    lineSpec.indfact = p.Results.indfactLine;
+    lineSpec.indvel = p.Results.indvelLine;
+    lineSpec.stahlhut = p.Results.stahlhutLine;
 end
 function vect = vectorize(struct, prop)

--- a/src/classes/@Rotor/Rotor.m
+++ b/src/classes/@Rotor/Rotor.m
@@ -92,18 +92,22 @@ classdef Rotor < handle
    methods
-        function self = Rotor(nBlades, Af, rad, chord, twist, iAf, nElem, position, name)
+        function self = Rotor(nBlades, BaseAf, afInterp, rad, chord, twist, iAf, nElem, ...
+                              position, name)
            % ROTOR Constructor.
            %   Constructs the object, stores the handle of the airfoil used and creates the blade
            %   object to represent the elements. See main class help for details.
            % Note: Angles should be passed IN DEGREE.
            if nargin > 0
+                import af_tools.Airfoil
                self.nBlades = nBlades;
-                if nargin >= 8
+                if nargin >= 9
                    self.position = position;
                end
-                if nargin == 9
+                if nargin == 10
                    self.name = name;
                end
@@ -112,9 +116,59 @@ classdef Rotor < handle
                self.cutout = rad(1);
                self.r0 = rad(1) / rad(end);
-                % Saves airfoil and create discretization
+                % Create blade discretization
-                self.Af = Af;
+                self.Bl = Blade(nBlades, rad, chord, twist, iAf, nElem, afInterp);
-                self.Bl = Blade(nBlades, rad, chord, twist, iAf, nElem);
+                % Interpolate airfoils (polars) between reference sections or use only reference
+                if afInterp
+                    self.Af = Airfoil;
+                    self.Af(1) = BaseAf(iAf(1));
+                    newiAf = ones(1, nElem);
+                    for i = 2:nElem - 1
+                        afCount = numel(self.Af);
+                        % Get previous/next refSection
+                        prevRef = find(self.Bl.y(i) > rad, 1, 'last');
+                        nextRef = find(self.Bl.y(i) < rad, 1, 'first');
+                        % Interpolate airfoils only when necessary
+                        if iAf(prevRef) ~= iAf(nextRef) % Interp between two reference sections
+                            % Calculate weight of second refsection
+                            weightNextAf = (self.Bl.y(i) - rad(prevRef)) / ...
+                                (rad(nextRef) - rad(prevRef));
+                            % Interpolate airfoils between reference sections
+                            self.Af(afCount + 1) = Airfoil.interpairfoil(BaseAf(iAf(prevRef)), ...
+                                                                         BaseAf(iAf(nextRef)), ...
+                                                                         weightNextAf);
+                            self.Af(afCount + 1).Polar.analyze();
+                            self.Af(afCount + 1).Polar.extrapMethod = BaseAf(1).Polar.extrapMethod;
+                            newiAf(i) = newiAf(i - 1) + 1;
+                        elseif self.Af(afCount) ~= BaseAf(iAf(prevRef)) % New airfoil Section
+                            self.Af(afCount + 1) = BaseAf(iAf(prevRef));
+                            newiAf(i) = newiAf(i - 1) + 1;
+                        else  % No need to interpolate
+                            newiAf(i) = newiAf(i - 1);
+                        end
+                    end
+                    % Last airfoil
+                    if  iAf(end) ~= iAf(end - 1)
+                        self.Af(end + 1) = BaseAf(iAf(end));
+                        newiAf(end) = newiAf(end - 1) + 1;
+                    else
+                        newiAf(end) = newiAf(end - 1);
+                    end
+                    % Optimize Blade structure to limit calls to interpolated airfoils
+                    self.Bl.optimizeairfoils(newiAf);
+                else
+                    self.Af = BaseAf;
+                end
            end
        end

--- a/src/configs/caradonna1981.m
+++ b/src/configs/caradonna1981.m
@@ -49,7 +49,6 @@ Sim.Misc.appli  = 'heli'; % Type of application ('helicopter', 'propeller', 'win
 % Solvers
 Mod.solvers = {'all'};   % BEMT Solver ('leishman', 'indfact', 'indvel', 'stahlhut', 'all')
-% Mod.solvers = {'leishman'};   % BEMT Solver ('leishman', 'indfact', 'indvel', 'stahlhut', 'all')
 % Extensions/corrections
 Mod.Ext.losses = 'all';  % Include losses using Prandtl formula ('none', 'hub', 'tip', 'both')
@@ -91,7 +90,7 @@ Airfoil.polarType = 'file'; % Type of polar to use ('file', 'polynomial')
 Airfoil.coordFile = 'airfoil_data/naca0012.dat';
 % If Airfoil.polarType == 'file'
-Airfoil.polarFile = 'airfoil_data/NACA_0012-Re_1e5-1e7.mat';
+Airfoil.polarFile = 'airfoil_data/NACA_0012-Re_2e5-1e7.mat';
 Airfoil.extrapMethod = 'viterna'; % Polar extrapol. ('none', 'spline', 'Viterna')
 % ==================================================================================================
@@ -110,9 +109,18 @@ Blade.radius   = [0.1905, 1.143];  % Spanwise position of the blade base station
 Blade.chord    = [0.1905, 0.1905]; % Chord at the blade base stations; [m]
 Blade.twist    = [0, 0];           % Twist at the blade base stations, [deg]
 Blade.iAirfoil = [1, 1];           % Index of the airfoil to use for the base stations, [-]
+% Interpolate the airfoil between sections (slow) or apply same airfoil until next section (fast)
+Blade.afInterp = false;
 % Discretization
 Blade.nElem = 100;  % Number of blade elements, [-]
 % Rotor base position
 Blade.hubPos = [0, 0, 0]; % Rotor center position (used for coaxial rotors), [m]
+% ==================================================================================================
+% ======================================= Overwrites ===============================================
+% ==================================================================================================
+% Minor overwrites to get more insightful filenames
+Sim.Save.filename = [Sim.Save.filename, '-loss_', Mod.Ext.losses];
--- a/src/configs/knight1937.m
+++ b/src/configs/knight1937.m
@@ -104,6 +104,8 @@ Blade.radius   = [0.127, 0.7620];  % Spanwise position of the blade base station
 Blade.chord    = [0.0508, 0.0508]; % Chord at the blade base stations; [m]
 Blade.twist    = [0, 0];           % Twist at the blade base stations, [deg]
 Blade.iAirfoil = [1, 1];           % Index of the airfoil to use for the base stations, [-]
+% Interpolate the airfoil between sections (slow) or apply same airfoil until next section (fast)
+Blade.afInterp = false;
 % Discretization
 Blade.nElem = 100;  % Number of blade elements, [-]

--- a/src/configs/template.m
+++ b/src/configs/template.m
@@ -59,7 +59,7 @@
 Sim.Save.autosave    = true;       % Auto-save the simulation results in a mat file
 Sim.Save.overwrite   = false;       % Overwrite previous result if filename is the same
 Sim.Save.dir         = '../results/';  % Directory where the results are saved
-Sim.Save.filename    = 'tempalte';  % File name of the saved result
+Sim.Save.filename    = 'template';  % File name of the saved result
 % Outputs
 Sim.Out.showPlots = true;  % Show all plots (forces, angles, speed, ...)
@@ -156,6 +156,8 @@ Blade.radius   = [0.1, 0.4, 2];     % Spanwise position of the blade base statio
 Blade.chord    = [0.3, 0.27, 0.10]; % Chord at the blade base stations; [m]
 Blade.twist    = [40, 30, 3];       % Twist at the blade base stations, [deg]
 Blade.iAirfoil = [1, 2, 2];         % Index of the airfoil to use for the base stations, [-]
+% Interpolate the airfoil between sections (slow) or apply same airfoil until next section (fast)
+Blade.afInterp = false;
 % Discretization
 Blade.nElem = 100;  % Number of blade elements, [-]

--- a/matlab_airfoil_toolbox @ 0995fa2a
+++ b/matlab_airfoil_toolbox @ 0995fa2a
-Subproject commit ec8f9f9b31c077d1e1082c1e0a85d72dfa3a2c06
+Subproject commit 0995fa2ae733a5f7df81c6d9eb0cae7b801b9e2e
No results found