plotwindonoff.m

function plotwindonoff(ResData)
    % 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
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    POSITION = 'Front';
    INTERESTING_POINTS = [
                          0.10, 2.5, 4.6
                          0.10, 2.5, 7.7
                          0.10, 3.5, 2.5
                          0.10, 3.5, 4.6
                          0.10, 3.5, 7.7
                         ];

    PLOT_TIME = false;
    SAVE_TIKZ_PARAM = INTERESTING_POINTS(1:2, :);

    for i = 1:size(ResData, 1)
        for j = 1:size(ResData, 2)

            if isinteresting(ResData(i, j), INTERESTING_POINTS)

                expParam = sprintf('d = %0.2f, f = %0.2f, V = %0.1f', ...
                                   ResData(i, j).Testcase.dx, ...
                                   ResData(i, j).(POSITION).trueFreq, ...
                                   ResData(i, j).Testcase.airspeed);

                % Index of corresponding wind-off data
                woIdx = (j) - mod(j - 1, 4);

                windOff = ResData(i, woIdx).(POSITION);
                timeOff = ResData(i, woIdx).time;
                sampling = ResData(i, woIdx).sampling;
                windOn = ResData(i, j).(POSITION);
                timeOn = ResData(i, j).time;

                meanPeakAngle = averagepeaks(windOff.angles, sampling, windOff.trueFreq);
                meanPeakFoff = averagepeaks(windOff.forces(:, 3), sampling, windOff.trueFreq);
                meanPeakFon = averagepeaks(windOn.forces(:, 3), sampling, windOn.trueFreq);

                % Trim to same length if needed
                minLen = min(length(meanPeakFoff), length(meanPeakFon));
                meanPeakAngle = meanPeakAngle(1:minLen);
                meanPeakFoff = meanPeakFoff(1:minLen);
                meanPeakFon = meanPeakFon(1:minLen);
                meanDiff = meanPeakFon - meanPeakFoff;

                % True period
                period = 1 / windOff.trueFreq;
                nPeriodIdx = period * sampling;

                if PLOT_TIME
                    time = (0:length(meanDiff) - 1) / nPeriodIdx;

                    figure;
                    hold on;
                    plot(time, meanPeakFoff);
                    plot(time, meanPeakFon);
                    plot(time, meanDiff);
                    yyaxis right;
                    plot(time, meanPeakAngle);
                    title(expParam);
                    legend('Wind off', 'Wind on', 'Diff');
                    grid on;
                    xlim([0, 1]);
                end

                figure;
                ax1 = subplot(3, 1, 1);
                setcolormap();
                hold on;

                plot(meanPeakAngle, meanPeakFon);
                % plot(meanPeakAngle(round(end / 5)), meanPeakFon(round(end / 5)), 'or');
                setgca('XTickLabel', [], ...
                       'YLabel', [], ...
                       'XLabel', []);

                ax2 =  subplot(3, 1, 2);
                setcolormap();
                hold on;
                plot(meanPeakAngle, meanPeakFoff);
                % plot(meanPeakAngle(round(end / 5)), meanPeakFoff(round(end / 5)), 'or');
                setgca('XTickLabel', [], ...
                       'XLabel', []);

                ax3 =  subplot(3, 1, 3);
                setcolormap();
                hold on;
                plot(meanPeakAngle, meanDiff);
                % plot(meanPeakAngle(round(end / 5)), meanDiff(round(end / 5)), 'or');
                linkaxes([ax1 ax2 ax3], 'xy');
                ax1.YLim = 1.25 * ax1.YLim;
                setgca('YLabel', []);

                if isinteresting(ResData(i, j), SAVE_TIKZ_PARAM)
                    figdir = 'figures/results/';
                    filename = sprintf('inertia-Dx%0.2f_F%0.2f_V%0.1f.tex', ...
                                       ResData(i, j).Testcase.dx, ...
                                       ResData(i, j).Testcase.freq, ...
                                       ResData(i, j).Testcase.airspeed);
                    save2tikz([figdir, filename], '\small');
                end
                % Add titles after having saved
                title(ax1, ['Wind on - ', expParam]);
                title(ax2, ['Wind off - ', expParam]);
                title(ax3, ['Diff - ', expParam]);
            end

        end
    end

end

function meanPeak = averagepeaks(signal, sampling, freq)
    % AVERAGEPEAKS Make the average of all peaks in the signal

    period = 1 / freq;
    nPeriodIdx = period * sampling;
    peakLocs = getallpeaks(signal, sampling);

    % Remove first and last peaks
    peakLocs = peakLocs(2:end - 1);

    % Create mean peak
    for i = 1:min(length(peakLocs), 100)
        peakVals(i, :) = signal(round(peakLocs(i) - nPeriodIdx / 2): ...
                                round(peakLocs(i) + nPeriodIdx / 2));
    end
    meanPeak = mean(peakVals);

end

function setgca(varargin)

    xlabel('Flapping angle [deg]');
    ylabel('\cfz [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