diff --git a/utils/analysis/plotstatictests.m b/utils/analysis/plotstatictests.m
new file mode 100644
index 0000000000000000000000000000000000000000..12c51425b6be6953e4264d1253d6bdb076f45c72
--- /dev/null
+++ b/utils/analysis/plotstatictests.m
@@ -0,0 +1,128 @@
+function ParsedDat = plotstatictests(StatData)
+    % PLOTSTATICTESTS Plot CL and CD for each set of wing w.r.t front and aft dihedral
+    %   Todo
+
+    % ----------------------------------------------------------------------------------------------
+    % (c) Copyright 2022 University of Liege
+    % Author: Thomas Lambert <t.lambert@uliege.be>
+    % ULiege - Aeroelasticity and Experimental Aerodynamics
+    % MIT License
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    close all;
+
+    % Constants and defaults
+    CHORD = 0.05;
+    SPAN = 2 * 0.2;
+    C_TO_K = 273.15;
+    AIRSPEEDS = [5, 8, 10, 14];
+
+    % Create useful vectors for each sets of front/aft position
+    % Vectors dimensions = number of speeds
+
+    % First we get all unique angles for front and aft wings that were tested
+    for i = numel(StatData):-1:1
+        angleF(i) = StatData(i).Front.angles(1);
+        angleA(i) = StatData(i).Aft.angles(1);
+    end
+    angleF = unique(angleF);
+    angleA = unique(angleA);
+    nVeloc = numel(StatData) / (length(angleF) * length(angleA)); % Number of velocities per config
+
+    % Create base mesh grid for all tests cases
+    [gridF, gridA] = meshgrid(angleF, angleA);
+
+    % Initialize force coefficients in X and Y for Front and Aft
+    cfzF = zeros([size(gridF), nVeloc]);
+    cfxF = zeros([size(gridF), nVeloc]);
+    cfzA = zeros([size(gridF), nVeloc]);
+    cfxA = zeros([size(gridF), nVeloc]);
+    stdfz = zeros([size(gridF), nVeloc]);
+
+    for i = 1:numel(StatData)
+
+        idx1 = gridF == StatData(i).Front.angles(1) & gridA == StatData(i).Aft.angles(1);
+        [~, idx2] = min(abs(AIRSPEEDS - StatData(i).Testcase.airspeed));
+
+        for iPos = {'Front', 'Aft'}
+            pos = string(iPos);
+
+            fz = mean(StatData(i).(pos).forces(:, 3));
+            fx = mean(StatData(i).(pos).forces(:, 1));
+            %             stdfz = std(StatData(i).(pos).forces(:, 3));
+            %             stdfx = std(StatData(i).(pos).forces(:, 1));
+
+            dens = airdensity(StatData(i).Testcase.temp + C_TO_K, ...
+                              StatData(i).Testcase.press);
+            speed = StatData(i).Testcase.airspeed;
+
+            cfz = forcetocoeff(fz, dens, speed, CHORD, SPAN, 'force');
+            cfx = forcetocoeff(fx, dens, speed, CHORD, SPAN, 'force');
+
+            % Nondim with projected area
+            cfz = cfz / sind(90 + StatData(i).(pos).angles(1));
+            cfx = cfx / sind(90 + StatData(i).(pos).angles(1));
+
+            if strcmp(pos, 'Front')
+                cfzF(:, :, idx2) = cfzF(:, :, idx2) + idx1 * fz;
+                cfxF(:, :, idx2) = cfxF(:, :, idx2) + idx1 * fx;
+                cfzF(:, :, idx2) = cfzF(:, :, idx2) + idx1 * fz;
+                cfxF(:, :, idx2) = cfxF(:, :, idx2) + idx1 * fx;
+            else
+                cfzA(:, :, idx2) = cfzA(:, :, idx2) + idx1 * fz;
+                cfxA(:, :, idx2) = cfxA(:, :, idx2) + idx1 * fx;
+            end
+
+        end
+
+    end
+
+    ParsedDat.cfzF = cfzF;
+    ParsedDat.cfxF = cfxF;
+    ParsedDat.cfzA = cfzA;
+    ParsedDat.cfxA = cfxA;
+
+    for i = 1:nVeloc
+        figure('Name', ['V=', i]);
+        ax1 = subplot(221); % Front CL
+        surf(gridF, gridA, cfzF(:, :, i));
+        setgca();
+
+        ax2 = subplot(222); % Front CD
+        surf(gridF, gridA, cfxF(:, :, i));
+        setgca();
+
+        ax3 = subplot(223); % Aft CL
+        surf(gridF, gridA, cfzA(:, :, i));
+        setgca();
+
+        ax4 = subplot(224); % Aft CD
+        surf(gridF, gridA, cfxA(:, :, i));
+        setgca();
+
+        %         linkprop([ax1 ax3], 'Zlim'); % CL plots
+        %         linkprop([ax2 ax4], 'Zlim'); % CD plots
+    end
+
+    for i = 1:length(angleF)
+        figure;
+        size(angleA);
+        size(cfzF(i, :, 4));
+        hold on;
+        title([angleF(i)]);
+        errobar(angleA, cfzF(:, i, 3), stdfzF(:, i, 3), '-o');
+        plot(angleA, cfzF(:, i, 3), '-o');
+        plot(angleA, cfzA(:, i, 3), '-o');
+        plot(angleA, cfzF(:, i, 3) + cfzA(:, i, 3), '--x');
+        hold on;
+    end
+
+end
+
+function setgca()
+    % SETGCA Sets the axes parameters
+
+    xlabel('Front');
+    ylabel('Aft');
+    grid on;
+end
diff --git a/utils/plotsemiflap.m b/utils/plotsemiflap.m
new file mode 100644
index 0000000000000000000000000000000000000000..499862df915e2506cf20975ddd3acdcb87ff0499
--- /dev/null
+++ b/utils/plotsemiflap.m
@@ -0,0 +1,418 @@
+function ResData = plotsemiflap(type)
+    % PLOTWINDONON Plots the forces in wind on and wind off conditions
+    %   Todo
+
+    % ----------------------------------------------------------------------------------------------
+    % (c) Copyright 2022 University of Liege
+    % Author: Thomas Lambert <t.lambert@uliege.be>
+    % ULiege - Aeroelasticity and Experimental Aerodynamics
+    % MIT License
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    set(0, 'defaulttextinterpreter', 'latex');
+
+    POSITION = 'Front';
+    INTERESTING_POINTS = [
+                         ];
+    PLOT_TIME = true;
+    SAVE_TIKZ_PARAM = []; % INTERESTING_POINTS(1:2, :);
+    MYGREY = [186, 216, 221] / 255;
+
+    ALLPLOTS = false;
+    SAMPLING = 1000;
+    ERR_CAPSIZE = 3;
+
+    CHORD = 0.05;
+    SPAN = 2 * 0.2;
+    C_TO_K = 273.15;
+
+    switch type
+        case 'Front' % Front wing is flapping
+            dat3cm = '2022-3cm-flapFront.mat';
+            dat5cm = '2022-5cm-flapFront.mat';
+            notType = 'Aft';
+        case 'Aft' % Aft wing is flapping
+            dat3cm = '2022-3cm-flapAft.mat';
+            dat5cm = '2022-5cm-flapAft.mat';
+            notType = 'Front';
+    end
+
+    Tmp = load(['data/cleaned_MAT/2022/', dat3cm]);
+    ResData(1, :) = Tmp.ExpData;
+    Tmp = load(['data/cleaned_MAT/2022/', dat5cm]);
+    ResData(2, :) = Tmp.ExpData;
+
+    % Remove useless data
+    ResData = rmfield(ResData, {'arduStart', 'tunnelStart'}); % Useless
+
+    % FIXME: remove deprecated
+    ResData = rmfield(ResData, {'arduFile', 'tunnelFile', ...
+                                'dx', 'airspeed', 'freq'});
+
+    for i = 1:size(ResData, 1)
+        %     for i = 1
+        for j = 1:size(ResData, 2)
+
+            dens = airdensity(ResData(i, j).Testcase.temp + C_TO_K, ...
+                              ResData(i, j).Testcase.press);
+            speed = ResData(i, j).Testcase.airspeed;
+
+            % Calculate main frequency
+            ResData(i, j).(type).trueFreq = getmainfrequency(ResData(i, j).(type).angles, ...
+                                                             1 / diff(ResData(i, j).time(1:2)));
+            ResData(i, j).(notType).trueFreq = 0;
+
+            % Calculate inertia
+            woIdx = (j) - mod(j - 1, 4); % Corresponding wind-off index
+
+            [ResData(i, j).(type).inertiaF, ...
+                ResData(i, j).(type).inertiaM] = calcinertia(ResData(i, woIdx), type);
+            ResData(i, j).(notType).inertiaF = zeros(size(ResData(i, j).(type).inertiaF));
+            ResData(i, j).(notType).inertiaM = zeros(size(ResData(i, j).(type).inertiaM));
+
+            % Average everything over a mean period for the active wing
+            [avSignal, peakIdx] = averageperiods(ResData(i, j).(type).angles, ...
+                                                 SAMPLING, ...
+                                                 ResData(i, j).(type).trueFreq);
+            ResData(i, j).(type).avAngles = avSignal;
+            [ResData(i, j).(type).avFz, ResData(i, j).(type).stdFz] = getavsignal(ResData(i, j).(type).forces(:, 3), peakIdx); % mh:ignore_style (WONTFIX)
+            [ResData(i, j).(type).avFx, ResData(i, j).(type).stdFx] = getavsignal(ResData(i, j).(type).forces(:, 1), peakIdx); % mh:ignore_style (WONTFIX)
+            [ResData(i, j).(notType).avFz, ResData(i, j).(notType).stdFz] = getavsignal(ResData(i, j).(notType).forces(:, 3), peakIdx); % mh:ignore_style (WONTFIX)
+            [ResData(i, j).(notType).avFx, ResData(i, j).(notType).stdFx] = getavsignal(ResData(i, j).(notType).forces(:, 1), peakIdx); % mh:ignore_style (WONTFIX)
+
+            % Remove inertia and invert FX so it is oriented alongside drag diraction
+            ResData(i, j).(type).avFz = ResData(i, j).(type).avFz - ...
+                ResData(i, j).(type).inertiaF(3);
+            ResData(i, j).(type).avFx = -(ResData(i, j).(type).avFx - ...
+                                          ResData(i, j).(type).inertiaF(1));
+
+            % Calculate coefficient
+            for iPos = {'Aft', 'Front'}
+                pos = string(iPos);
+                ResData(i, j).(pos).cFz = forcetocoeff(ResData(i, j).(pos).avFz, ...
+                                                       dens, speed, CHORD, SPAN, 'force');
+                ResData(i, j).(pos).cFx = forcetocoeff(ResData(i, j).(pos).avFx, ...
+                                                       dens, speed, CHORD, SPAN, 'force');
+                ResData(i, j).(pos).stdFz = forcetocoeff(ResData(i, j).(pos).stdFz, ...
+                                                         dens, speed, CHORD, SPAN, 'force');
+                ResData(i, j).(pos).stdFx = forcetocoeff(ResData(i, j).(pos).stdFx, ...
+                                                         dens, speed, CHORD, SPAN, 'force');
+
+                % Calculate mean CL and mean CD over the average period
+                ResData(i, j).(pos).cFz_mean = mean(ResData(i, j).(pos).cFz);
+                ResData(i, j).(pos).cFx_mean = mean(ResData(i, j).(pos).cFx);
+                % Get mean STD over the period
+                ResData(i, j).(pos).cFz_std = mean(ResData(i, j).(pos).stdFz);
+                ResData(i, j).(pos).cFx_std = mean(ResData(i, j).(pos).stdFx);
+            end
+
+            % -----------------------------------------------
+            % Plot each test case
+            if ALLPLOTS
+                redFreq = redfreq(ResData(i, j).(type).trueFreq, ...
+                                  CHORD, ResData(i, j).Testcase.airspeed);
+
+                expParam = sprintf(['Aft angle = %0.1f - d = %0.2f, f = %0.2f, '...
+                                    'V = %0.1f (k = %0.02f)'], ...
+                                   ResData(i, j).(notType).angles(1), ...
+                                   ResData(i, j).Testcase.dx, ...
+                                   ResData(i, j).(type).trueFreq, ...
+                                   ResData(i, j).Testcase.airspeed, ...
+                                   redFreq);
+
+                [~, donwStrokeEndIdx] = min(ResData(i, j).(type).avAngles);
+                nPeriodIdx = 1 / ResData(i, j).(type).trueFreq * SAMPLING;
+                time = (0:length(ResData(i, j).(type).avAngles) - 1) / nPeriodIdx;
+
+                patchX = [0, donwStrokeEndIdx / nPeriodIdx, donwStrokeEndIdx / nPeriodIdx, 0];
+                patchY = [-5, -5, 10, 10];
+
+                figure('Name', 'HalfFlap');
+                ax1 = subplot(211);
+                setcolormap();
+                hold on;
+                patch(patchX, patchY, [186, 216, 221] / 255, 'EdgeColor', 'none', 'FaceAlpha', 0.5);
+                h1 = plot(time, ResData(i, j).Front.avFz);
+                h2 = plot(time, ResData(i, j).Front.avFx);
+                hold off;
+                setgca('XTickLabel', [], ...
+                       'XLabel', []);
+                legend([h1, h2], {'\fz', '\fx'}, 'Location', 'NorthEast');
+                legend('boxoff');
+                % ax1.Position(3) = 1.25 * ax1.Position(3);
+
+                ax2 = subplot(212);
+                setcolormap();
+                hold on;
+                patch(patchX, patchY, [186, 216, 221] / 255, 'EdgeColor', 'none', 'FaceAlpha', 0.5);
+                h1 = plot(time, ResData(i, j).Aft.avFz);
+                h2 = plot(time, ResData(i, j).Aft.avFx);
+
+                hold off;
+
+                setgca();
+                linkaxes([ax1 ax2], 'xy');
+                ax2.Position(3) = ax1.Position(3);
+
+                % Increase axes Y lim by 10% of the data range before writing text
+                FACTOR = 0.1;
+
+                yDiff = diff(ax1.YLim);
+                ax1.YLim = [ax1.YLim(1), ax1.YLim(2) + FACTOR * yDiff];
+                addtxt(ax1, '\bf Front', ax1.YLim(2));
+                addtxt(ax2, '\bf Aft', ax2.YLim(2));
+
+                text(ax2, (donwStrokeEndIdx / nPeriodIdx) / 2, 0.5, ...
+                     'Downstroke', ...
+                     'HorizontalAlignment', 'center');
+                text(ax2, donwStrokeEndIdx / nPeriodIdx + ...
+                     (1 - (donwStrokeEndIdx / nPeriodIdx)) / 2, 0.5, ...
+                     'Upstroke', ...
+                     'HorizontalAlignment', 'center');
+
+                %
+                %             if isinteresting(ResData(i, j), SAVE_TIKZ_PARAM)
+                %                 pause(1);
+                %                 figdir = 'figures/results/';
+                %                 filename = sprintf('flap%s-Dx%0.2f_F%0.2f_V%0.1f.tex', ...
+                %                     type, ...
+                %                     ResData(i, j).Testcase.dx, ...
+                %                     ResData(i, j).Testcase.freq, ...
+                %                     ResData(i, j).Testcase.airspeed);
+                %                 save2tikz([figdir, filename], '\small');
+                %             end
+                % Update titles after having saved
+                title(ax1, ['Front - ', expParam]);
+                title(ax2, ['Aft - ', expParam]);
+            end
+        end
+    end
+
+    % Loop on all cases and create useful vectors
+    nSpeeds = 4; % Number of airspeeds tested
+    for iDx = 1:size(ResData, 1)
+        for j = 2:nSpeeds
+
+            speedIdx = j:4:size(ResData, 2);
+            for k = 1:length(speedIdx)
+                speeds(j, k) = ResData(iDx, speedIdx(k)).Testcase.airspeed;
+                freq(j, k) = ResData(iDx, speedIdx(k)).Testcase.freq;
+                dx(j, k) = ResData(iDx, speedIdx(k)).Testcase.dx;
+                ang2(j, k) = mean(ResData(iDx, speedIdx(k)).Aft.angles);
+                cfzF(j, k) = ResData(iDx, speedIdx(k)).Front.cFz_mean;
+                cfxF(j, k) = ResData(iDx, speedIdx(k)).Front.cFx_mean;
+                cfzA(j, k) = ResData(iDx, speedIdx(k)).Aft.cFz_mean;
+                cfxA(j, k) = ResData(iDx, speedIdx(k)).Aft.cFx_mean;
+                stdcfzF(j, k) = ResData(iDx, speedIdx(k)).Front.cFz_std;
+                stdcfxF(j, k) = ResData(iDx, speedIdx(k)).Front.cFx_std;
+                stdcfzA(j, k) = ResData(iDx, speedIdx(k)).Aft.cFz_std;
+                stdcfxA(j, k) = ResData(iDx, speedIdx(k)).Aft.cFx_std;
+            end
+
+            figure('Name', 'Dx');
+            ax1 = subplot(221);
+            setcolormap();
+            hold on;
+            e11 = errorbar(ang2(j, 1:3), cfzF(j, 1:3), ...
+                           stdcfzF(j, 1:3), '-s', 'Capsize', ERR_CAPSIZE);
+            e12 = errorbar(ang2(j, 1:3), cfxF(j, 1:3), ...
+                           stdcfxF(j, 1:3), '-o', 'Capsize', ERR_CAPSIZE);
+            e11.MarkerFaceColor = e11.Color;
+            e12.MarkerFaceColor = e12.Color;
+            hold off;
+            setgca2('XTickLabel', [], ...
+                    'XLabel', []);
+
+            legend({'\cl', '\cd'}, 'Location', 'SouthWest');
+            legend('boxoff');
+
+            ax2 = subplot(222);
+            setcolormap();
+            hold on;
+            e21 = errorbar(ang2(j, 4:end), cfzF(j, 4:end), ...
+                           stdcfzF(j, 4:end), '-s', 'Capsize', ERR_CAPSIZE);
+            e22 = errorbar(ang2(j, 4:end), cfxF(j, 4:end), ...
+                           stdcfxF(j, 4:end), '-o', 'Capsize', ERR_CAPSIZE);
+            e21.MarkerFaceColor = e21.Color;
+            e22.MarkerFaceColor = e22.Color;
+            hold off;
+            setgca2('XTickLabel', [], ...
+                    'XLabel', [], ...
+                    'YTickLabel', [], ...
+                    'YLabel', []);
+            ax3 = subplot(223);
+            setcolormap();
+            hold on;
+            cosd(ang2(j, 4:end));
+            e31 = errorbar(ang2(j, 1:3), cfzA(j, 1:3) ./ cosd(ang2(j, 4:end)), ...
+                           stdcfzA(j, 1:3), '-s', 'Capsize', ERR_CAPSIZE);
+            e32 = errorbar(ang2(j, 1:3), cfxA(j, 1:3) ./ cosd(ang2(j, 4:end)), ...
+                           stdcfxA(j, 1:3), '-o', 'Capsize', ERR_CAPSIZE);
+            e31.MarkerFaceColor = e31.Color;
+            e32.MarkerFaceColor = e32.Color;
+            hold off;
+            setgca2();
+            ylabel('$\cl^*$ and $\cd^*$ [-]');
+
+            ax4 = subplot(224);
+            setcolormap();
+            hold on;
+            e41 = errorbar(ang2(j, 4:end), cfzA(j, 4:end) ./ cosd(ang2(j, 4:end)), ...
+                           stdcfzA(j, 4:end), '-s', 'Capsize', ERR_CAPSIZE);
+            e42 = errorbar(ang2(j, 4:end), cfxA(j, 4:end) ./ cosd(ang2(j, 4:end)), ...
+                           stdcfxA(j, 4:end), '-o', 'Capsize', ERR_CAPSIZE);
+            e41.MarkerFaceColor = e41.Color;
+            e42.MarkerFaceColor = e42.Color;
+            hold off;
+            setgca2('YTickLabel', [], ...
+                    'YLabel', []);
+
+            linkaxes([ax1 ax2 ax3 ax4], 'xy');
+
+            % Add annotations
+            % Increase axes Y lim by 10% of the data range before writing text
+            FACTOR = 0.1;
+
+            yDiff = diff(ax1.YLim);
+            ax1.YLim = [ax1.YLim(1), ax1.YLim(2) + FACTOR * yDiff];
+            addtxt(ax1, '$\mathbf{\freq}$ = 1.25 Hz (Front)', ax1.YLim(2));
+            addtxt(ax2, '$\mathbf{\freq}$ = 2.00 Hz (Front)', ax2.YLim(2));
+            addtxt(ax3, '$\mathbf{\freq}$ = 1.25Hz (Aft)', ax3.YLim(2));
+            addtxt(ax4, '$\mathbf{\freq}$ = 2.00 Hz (Aft)', ax4.YLim(2));
+            pause(1);
+            figdir = 'figures/results/';
+            filename = sprintf('half%s-AllF-Dx%0.2f_V%0.1f.tex', ...
+                               type, ...
+                               dx(j, 1), ...
+                               speeds(j, 1));
+            save2tikz([figdir, filename], '\small');
+
+            expParam = sprintf('Dx = %0.02f, U = %0.2f', ...
+                               dx(j, 1), ...
+                               speeds(j, 1));
+            title(ax1, expParam);
+            title(ax3, expParam);
+
+        end
+
+    end
+
+end
+
+function addtxt(ax, txt, ypos)
+
+    text(ax, -10, ypos, ...
+         txt, ...
+         'HorizontalAlignment', 'left', 'VerticalAlignment', 'top', 'FontWeight', 'bold');
+end
+
+function [inertiaF, inertiaM] = calcinertia(ResDataWO, type)
+    % CALCINERTIA Calcualte inertia forces
+
+    SAMPLING = 1000;
+
+    % Get first and last peaks for wind-off conditions
+    woBounds = getboundingpeaks(ResDataWO.(type).angles, SAMPLING);
+    woIdx = woBounds(1):woBounds(2);
+
+    % Calculate mean force and moments due to inertia
+    inertiaF = mean(ResDataWO.(type).forces(woIdx, :));
+    inertiaM = mean(ResDataWO.(type).moments(woIdx, :));
+end
+
+function [avSignal, peakIdx] = averageperiods(signal, sampling, freq)
+    % AVERAGEPERIOD Average all periods of the signal and return the average signal
+
+    % Find all peaks location
+    period = 1 / freq;
+    nPeriodIdx = period * sampling;
+    peakLocs = getallpeaks(signal, sampling);
+
+    % Remove first and last peaks
+    peakLocs = peakLocs(2:end - 1);
+
+    % Split signal in a matrix whith rows for the peak number and cols to the INDEXES of values in
+    for iPeak = length(peakLocs):-1:1
+        peakIdx(iPeak, :) = peakLocs(iPeak):peakLocs(iPeak) + round(nPeriodIdx);
+        peakSignal(iPeak, :) = signal(peakIdx(iPeak, :));
+    end
+
+    % Take the mean value of the signal formed with all these peaks
+    avSignal = mean(peakSignal);
+end
+
+function [avSignal, stdSignal] = getavsignal(signal, peakIdx)
+
+    for iPeak = size(peakIdx, 1):-1:1
+        peakSignal(iPeak, :) = signal(peakIdx(iPeak, :));
+    end
+
+    % Take the mean value of the signal formed with all these peaks
+    avSignal = mean(peakSignal);
+    stdSignal = std(peakSignal);
+
+end
+
+function setgca(varargin)
+
+    ylim([-1.5, 1.5]);
+    xlim([0, 1]);
+
+    xlabel('t/T [-]');
+    ylabel('\fz and \fx [N]');
+
+    set(gca, ...
+        'Box', 'off', ...
+        'TickDir', 'out', ...
+        'TickLength', [.02 .02], ...
+        'XMinorTick', 'off', ...
+        'YMinorTick', 'off', ...
+        'YGrid', 'off', ...
+        'XGrid', 'off', ...
+        'XColor', 'k', ...
+        'YColor', 'k', ...
+        'GridLineStyle', ':', ...
+        'GridColor', 'k', ...
+        'GridAlpha', 0.25, ...
+        'LineWidth', 1, ...
+        'FontName', 'Helvetica', ...
+        'Fontsize', 14);
+
+    if nargin > 0
+        set(gca, varargin{:});
+    end
+
+end
+
+function setgca2(varargin)
+
+    %     ylim([-1.5, 1.5])
+    xlim([-13, 40]);
+    xticks([-10:10:40]);
+
+    %
+    xlabel('Aft wing dihedral [deg]');
+    ylabel('\cl and \cd [-]');
+
+    set(gca, ...
+        'Box', 'off', ...
+        'TickDir', 'out', ...
+        'TickLength', [.02 .02], ...
+        'XMinorTick', 'off', ...
+        'YMinorTick', 'off', ...
+        'YGrid', 'off', ...
+        'XGrid', 'off', ...
+        'XColor', 'k', ...
+        'YColor', 'k', ...
+        'GridLineStyle', ':', ...
+        'GridColor', 'k', ...
+        'GridAlpha', 0.25, ...
+        'LineWidth', 1, ...
+        'FontName', 'Helvetica', ...
+        'Fontsize', 14);
+
+    if nargin > 0
+        set(gca, varargin{:});
+    end
+
+end