.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/*

.gitlab-ci.yml

@@ -27,7 +27,7 @@ style_check:
   <<: *miss_hit
   stage: static analysis
   script:
-    - mh_style . --fix
+    - mh_style .
 
 metric_check:
   <<: *miss_hit

.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

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

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
 

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)

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

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"

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:
-        print("%d/%d - Calculating polar for Re=%ge6" % (i, len(reynolds), re / 1e6))
+    for re in REYNOLDS:
+        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")
 
 

src/classes/@Blade/Blade.m

@@ -19,31 +19,35 @@ classdef Blade < handle
     % -----
     %
     % Blade properties:
-    %   nElem - Number of elements, [-]
-    %   dy    - Element span, [m]
-    %   y     - Element absolute radial position, [m]
-    %   r     - Element relative radial position (relative to the total radius), [-]
-    %   area  - Element area, [m^2]
-    %   chord - Element chord, [m]
-    %   twist - Element twist (also called stagger angle), [rad]
-    %   iAf   - Index of airfoil to use for each element, [-]
-    %   sol   - Element solidity, [-]
+    %   nElem    - Number of elements, [-]
+    %   dy       - Element span, [m]
+    %   y        - Element absolute radial position, [m]
+    %   r        - Element relative radial position (relative to the total radius), [-]
+    %   area     - Element area, [m^2]
+    %   chord    - Element chord, [m]
+    %   twist    - Element twist (also called stagger angle), [rad]
+    %   iAf      - Index of airfoil to use for each element, [-]
+    %   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
-    %   provided (see list below).
+    %   Bl = Blade(nBlades, rad, chord, twist, iAf, nElem, afInterp) creates a Blade object based on
+    %   the input provided (see list below).
     %
     % Constructor inputs:
-    %   nBlades : Number of blades, [-]
-    %   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)
-    %   twist   : Twist of the guide stations, [deg] (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, [-]
+    %   nBlades  : Number of blades, [-]
+    %   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)
+    %   twist    : Twist of the guide stations, [deg] (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, [-]
+    %   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
-                    self.iAf(self.y >= rad(i) & self.y <= rad(i + 1)) = iAf(i);
+                if self.afInterp
+                    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

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)

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);
 

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

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)

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)
-        curOp = DEF.allowed_oper{i};
+    for i = 1:length(DEF.base_oper)
+        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)

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
-                self.Af = Af;
-                self.Bl = Blade(nBlades, rad, chord, twist, iAf, nElem);
+                % Create blade discretization
+                self.Bl = Blade(nBlades, rad, chord, twist, iAf, nElem, afInterp);
+
+                % 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
 

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];

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, [-]

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, [-]

matlab_airfoil_toolbox @ 0995fa2a

-Subproject commit ec8f9f9b31c077d1e1082c1e0a85d72dfa3a2c06
+Subproject commit 0995fa2ae733a5f7df81c6d9eb0cae7b801b9e2e