feat: add figures for single flapping

src/fig_setup.m

@@ -40,6 +40,7 @@ FORCE_RELOAD = 'false';  % If 'reload', force reload the results
 FIG_TYPE = 'dissert'; % Can either be 'dissertation' or 'presentation'
 LINE_WIDTH = 1.2;
 SAVE_FIGURES = true;
+KEEP_FIGURES_OPEN = false;
 
 MATLAB_FIG_DIR = '../figures/methodo/';
 THESIS_FIG_DIR = '../../../../../01-Thesis/figures/mecharaptor';

src/fig_singleFlap.m

+% FIG_SINGLEFLAP Create figures for single (front) wing flapping.
+%   This script load the appropriate data and generate the various figures used in my thesis.
+%   The figures are automatically saved using matlab2tikz.
+%
+%   The title and some annotations are then added afterwards for easier analysis when generated with
+%   Matlab, they may not all be present in the exported figure for clarity in the dissertation.
+% -----
+%
+% Syntax:
+%   FIG_SINGLEFLAP;
+%
+% Inputs:
+%   None
+%
+% Output:
+%   The following figures are created and saved
+%       - NACA 0012 - CL, CD vs time
+%       - NACA 0012 - CL, CD vs flapping angle
+%
+% See also: TestCase, TestResult.
+
+% ----------------------------------------------------------------------------------------------
+% TODO:
+% - Use same colors for same frequencies throughougt the graphs!
+% - Same for velocities
+% ----------------------------------------------------------------------------------------------
+% matlab2tikz: https://github.com/matlab2tikz/matlab2tikz
+% ----------------------------------------------------------------------------------------------
+% (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/mecharaptor/data-processing
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% --- Initial setup
+fig_setup;
+
+% --- Load data required for this file
+printsection('0. Data loading');
+
+RL.loadsingle(FORCE_RELOAD);
+logmsg('Data fully loaded', 2, VERBOSITY);
+
+pause(0.1); % So everything displays nicely
+
+% --- Defaults and constants
+ALL_FREQS = [1.25, 2, 2.5, 3, 3.5];
+ALL_AIRSPEED = [0, 2.5, 4.6, 7.7, 9.9];
+
+FIG_DIR = '../figures/results/singleFlap/';
+
+TEXT_POS = 0.03; % Text position as a % of the figure height, from the bottom line
+
+% [FIXME] Should come from the windtunnel data
+WING_CHORD = 0.05;
+WING_SPAN = 2 * 0.20;
+
+%% NACA 0012 - CL, CD vs time
+% - One figure to compare same frequency, different airspeed
+% - No point in comparing different frequencies as no normalisation
+
+printsection('1. CL, CD vs time');
+
+for iFreq = 1:numel(ALL_FREQS)
+    plotBkg = true;
+
+    thisFreq = ALL_FREQS(iFreq);
+    logmsg(sprintf('Plotting for f = %0.2f', thisFreq), 2, VERBOSITY);
+
+    % Figure - ONLY ONE WING
+    initiatefig('N0012 - CL vs Time', KEEP_FIGURES_OPEN);
+
+    hold on;
+    for i = 1:numel(RL.Single0012)
+        ThisRes = RL.Single0012(i);
+
+        if ThisRes.TestCond.airspeed == 0
+            continue
+        end
+
+        thisnondim = 1 / 2 * ThisRes.density * ThisRes.TestCond.airspeed^2 * ...
+            WING_CHORD * WING_SPAN;
+        iSpeed = find(ALL_AIRSPEED == ThisRes.TestCond.airspeed);
+
+        if ~ThisRes.rejected && ThisRes.TestCond.frequency == thisFreq
+            ThisResReorgF = ThisRes.PeriodsF.reorganize('upstroke');
+            fCond = sprintf('f = %0.2f Hz', ThisRes.TestCond.frequency);
+            airCond = sprintf('V = %0.2f m/s', ThisRes.TestCond.airspeed);
+
+            if plotBkg
+                [~, bgIdx] = max(ThisResReorgF.MedianData.angle);
+                timeUp = ThisResReorgF.MedianData.time(bgIdx) / ThisResReorgF.MedianData.time(end);
+                hPatch = bkgpatch(timeUp, 0.1);
+                plotBkg = false;
+            end
+
+            plot(ThisResReorgF.MedianData.time / ThisResReorgF.MedianData.time(end), ...
+                 ThisResReorgF.MedianData.Fz ./ thisnondim, ...
+                 'color', cmap(iSpeed, :), 'DisplayName', airCond);
+        end
+    end
+    hold off;
+
+    % Plot style
+    setgca();
+    figTitle = sprintf('F = %0.2f Hz', thisFreq);
+    ylabel('$\cfZ$ [-]');
+    xlabel('Time, t/T [-]');
+    legend;
+    yLimits = ylim;
+    hPatch.Vertices(:, 2) = [yLimits(1); yLimits(1); yLimits(2); yLimits(2)];
+    txtY = yLimits(1) + TEXT_POS * diff(yLimits);
+    bkgtxt({timeUp / 2, txtY, 'Upstroke'; 1 - timeUp / 2, txtY, 'Downstroke'});
+
+    % Save figure
+    savethisfig(SAVE_FIGURES, figTitle, FIG_DIR, sprintf('singleFlap-F%0.2f_vTime.tex', ...
+                                                         thisFreq));
+end
+
+%% NACA 0012 - CL, CD vs alpha
+% - One figure to compare same airspeed, different frequency
+% - One figure to compare same frequency, different airspeed
+
+% Only plot the following speeds and frequencies on the same graphs
+SELECTED_SPEEDS = [];
+SELECTED_FREQS = [];
+
+printsection('2. CL, CD vs flapping angle');
+
+for iFreq = 1:numel(ALL_FREQS)
+    thisFreq = ALL_FREQS(iFreq);
+    logmsg(sprintf('Plotting for f = %0.2f', thisFreq), 2, VERBOSITY);
+
+    for force = {'Fx', 'Fz'}
+        % Create figure
+        initiatefig(sprintf('ALPHA - %s', force{:}), KEEP_FIGURES_OPEN);
+        plotallspeeds(RL.Single0012, [], force{:}, thisFreq, ...
+                      SELECTED_SPEEDS, WING_CHORD, WING_SPAN, ALL_AIRSPEED, cmap);
+
+        % Save figure
+        savethisfig(SAVE_FIGURES, ...
+                    sprintf('%s - f = %0.2fHz', force{:}, thisFreq), ...
+                    FIG_DIR, ...
+                    sprintf('singleFlap-0012-%s-f%0.2fHz_vAngle.tex', force{:}, thisFreq));
+    end
+
+end
+
+for iAirSpeed = 1:numel(ALL_AIRSPEED)
+    thisSpeed = ALL_AIRSPEED(iAirSpeed);
+    logmsg(sprintf('Plotting for U = %0.2f', thisSpeed), 2, VERBOSITY);
+
+    for force = {'Fx', 'Fz'}
+        % Create figure
+        initiatefig(sprintf('ALPHA - %s', force{:}), KEEP_FIGURES_OPEN);
+        plotallfreqs(RL.Single0012, [], force{:}, thisSpeed, ...
+                     SELECTED_FREQS, WING_CHORD, WING_SPAN, ALL_FREQS, cmap);
+
+        % Save figure
+        savethisfig(SAVE_FIGURES, ...
+                    sprintf('%s - U = %0.2fm/s', force{:}, thisSpeed), ...
+                    FIG_DIR, ...
+                    sprintf('singleFlap-0012-%s-U%0.2fmps_vAngle.tex', force{:}, thisSpeed));
+    end
+
+end
+
+%% NACA 0012 vs 6409 - CL, CD vs alpha
+% - Comparison of NACA plantforms
+
+printsection('3. 0012 vs 6409');
+
+for iFreq = 1:numel(ALL_FREQS)
+    thisFreq = ALL_FREQS(iFreq);
+    logmsg(sprintf('Plotting for f = %0.2f', thisFreq), 2, VERBOSITY);
+
+    for force = {'Fx', 'Fz'}
+        % Create figure
+        initiatefig(sprintf('AIRFOILS - %s', force{:}), KEEP_FIGURES_OPEN);
+        plotallspeeds(RL.Single0012, RL.Single6409, force{:}, thisFreq, ...
+                      SELECTED_SPEEDS, WING_CHORD, WING_SPAN, ALL_AIRSPEED, cmap);
+
+        % Save figure
+        savethisfig(SAVE_FIGURES, ...
+                    sprintf('%s - f = %0.2fHz', force{:}, thisFreq), ...
+                    FIG_DIR, ...
+                    sprintf('singleFlap-multi-%s-f%0.2fHz_vAngle.tex', force{:}, thisFreq));
+    end
+
+end
+
+for iAirSpeed = 1:numel(ALL_AIRSPEED)
+    thisSpeed = ALL_AIRSPEED(iAirSpeed);
+    logmsg(sprintf('Plotting for U = %0.2f', thisSpeed), 2, VERBOSITY);
+
+    for force = {'Fx', 'Fz'}
+
+        % Create figure
+        initiatefig(sprintf('AIRFOILS - %s', force{:}), KEEP_FIGURES_OPEN);
+        plotallfreqs(RL.Single0012, RL.Single6409, force{:}, thisSpeed, ...
+                     SELECTED_FREQS, WING_CHORD, WING_SPAN, ALL_FREQS, cmap);
+
+        % Save figure
+        savethisfig(SAVE_FIGURES, ...
+                    sprintf('%s - U = %0.2fm/s', force{:}, thisSpeed), ...
+                    FIG_DIR, ...
+                    sprintf('singleFlap-multi-%s-U%0.2fmps_vAngle.tex', force{:}, thisSpeed));
+
+    end
+end
+
+%% MAINTENANCE / CHORES
+
+% Copy all images created here to the thesis directory for inclusion in latex document
+if SAVE_FIGURES
+    copydir2dir('../figures/results/', THESIS_FIG_DIR);
+end
+
+%% ================================== HELPER FUNCTIONS =============================================
+
+function initiatefig(figname, keepFigOpen)
+    % INITIATEFIG Initiate a new figure
+
+    if ~keepFigOpen
+        close(findobj('type', 'figure', 'name', figname));
+        pause(0.1);
+    elseif ~isempty(findobj('type', 'figure', 'name', figname))
+        i = 1;
+        newFigName = figname;
+        while ~isempty(findobj('type', 'figure', 'name', newFigName))
+            newFigName = sprintf('%s -%d', figname, i);
+            i = i + 1;
+        end
+        figname = newFigName;
+    end
+    figure('Name', figname);
+
+end
+
+function savethisfig(saveFig, figTitle, figDir, figName)
+    % SAVETHISFIG Save the figure
+
+    condsave2tikz(saveFig, figTitle, [figDir, figName], '\normalsize');
+    pause(0.1); % Do not show them all at once
+
+end
+
+function plotallspeeds(LoadedRes1, LoadedRes2, fType, frequency, ...
+                       selectedSpeeds, chord, span, airspeeds, cmap)
+    % PLOTALLSPEEDS Plot all airspeeds on the same figure
+
+    % Plot everything properly
+    hold on;
+    if isempty(LoadedRes2)
+        plotspeed(LoadedRes1);
+    else
+        plotspeed(LoadedRes1, '-', 'NACA0012');
+        plotspeed(LoadedRes2, '--', 'NACA0064');
+    end
+    hold off;
+
+    % Plot style
+    setgca();
+    xlabel('Flapping angle [°]');
+    if strcmp(fType, 'Fx')
+        ylabel('$\cfX$ [-]');
+    elseif strcmp(fType, 'Fz')
+        ylabel('$\cfZ$ [-]');
+    else
+        ylabel('$\cF$ [-]');
+    end
+    legend;
+
+    function plotspeed(LoadedRes, linestyle, extraLgd)
+        % PLOTFREQ Plot forces vs flapping angle for a single frequency
+
+        if nargin < 2
+            linestyle = '-';
+        end
+
+        hold on;
+        for i = 1:numel(LoadedRes)
+            ThisRes = LoadedRes(i);
+
+            % Name to be displayed on the legend
+            if nargin == 3
+                aCond = sprintf('U = %0.2f Hz, %s', ThisRes.TestCond.airspeed, extraLgd);
+            else
+                aCond = sprintf('U = %0.2f Hz', ThisRes.TestCond.airspeed);
+            end
+
+            if ~ThisRes.rejected && ThisRes.TestCond.frequency == frequency && ...
+                    ThisRes.TestCond.airspeed ~= 0 && ...
+                    (isempty(selectedSpeeds) || any(ThisRes.TestCond.airspeed == selectedSpeeds))
+                ThisResReorgF = ThisRes.PeriodsF.reorganize('upstroke');
+                thisnondim = 1 / 2 * ThisRes.density * ThisRes.TestCond.airspeed^2 * ...
+                    chord * span;
+                plot(ThisResReorgF.MedianData.angle, ...
+                     ThisResReorgF.MedianData.(fType) ./ thisnondim, ...
+                     linestyle, 'color', cmap(airspeeds == ThisRes.TestCond.airspeed, :), ...
+                     'DisplayName', aCond);
+            end
+        end
+        hold off;
+
+    end
+
+end
+
+function plotallfreqs(LoadedRes1, LoadedRes2, fType, airspeed, ...
+                      selectedFreqs, chord, span, freqs, cmap)
+    % PLOTALLFREQ Plot all frequencies on the same figure
+
+    % Plot everything properly
+    hold on;
+    if isempty(LoadedRes2)
+        plotfreq(LoadedRes1);
+    else
+        plotfreq(LoadedRes1, '-', 'NACA0012');
+        plotfreq(LoadedRes2, '--', 'NACA0064');
+    end
+
+    hold off;
+
+    % Plot style
+    setgca();
+    xlabel('Flapping angle [°]');
+    if strcmp(fType, 'Fx')
+        ylabel('$\cfX$ [-]');
+    elseif strcmp(fType, 'Fz')
+        ylabel('$\cfZ$ [-]');
+    else
+        ylabel('$\cF$ [-]');
+    end
+    legend;
+
+    function plotfreq(LoadedRes, linestyle, extraLgd)
+        % PLOTFREQ Plot forces vs flapping angle for a single frequency
+
+        if nargin < 2
+            linestyle = '-';
+        end
+
+        hold on;
+        for i = 1:numel(LoadedRes)
+            ThisRes = LoadedRes(i);
+
+            % Name to be displayed on the legend
+            if nargin == 3
+                fCond = sprintf('f = %0.2f Hz, %s', ThisRes.TestCond.frequency, extraLgd);
+            else
+                fCond = sprintf('f = %0.2f Hz', ThisRes.TestCond.frequency);
+            end
+
+            if ~ThisRes.rejected && ThisRes.TestCond.airspeed == airspeed && ...
+                    (isempty(selectedFreqs) || any(ThisRes.TestCond.frequency == selectedFreqs))
+                ThisResReorgF = ThisRes.PeriodsF.reorganize('upstroke');
+                thisnondim = 1 / 2 * ThisRes.density * ThisRes.TestCond.airspeed^2 * ...
+                    chord * span;
+                plot(ThisResReorgF.MedianData.angle, ...
+                     ThisResReorgF.MedianData.(fType) ./ thisnondim, ...
+                     linestyle, 'color', cmap(freqs == ThisRes.TestCond.frequency, :), ...
+                     'DisplayName', fCond);
+            end
+        end
+        hold off;
+
+    end
+
+end