diff --git a/utils/plotwindonoff.m b/utils/plotwindonoff.m index a69e88dbd0698911fa9f202fbd5f53a70772012b..7e957ad2130431369bf44a23a466601a1721675f 100644 --- a/utils/plotwindonoff.m +++ b/utils/plotwindonoff.m @@ -9,17 +9,15 @@ function plotwindonoff(ResData) % MIT License %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + set(0, 'defaulttextinterpreter', 'latex'); + 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 + 0.10, 2, 7.7 ]; - - PLOT_TIME = false; - SAVE_TIKZ_PARAM = INTERESTING_POINTS(1:2, :); + PLOT_TIME = true; + SAVE_TIKZ_PARAM = []; % INTERESTING_POINTS(1:2, :); + MYGREY = [186, 216, 221] / 255; for i = 1:size(ResData, 1) for j = 1:size(ResData, 2) @@ -40,44 +38,134 @@ function plotwindonoff(ResData) 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); + %% PROBABLY BETTER + windOffPeriodIdx = getmeanperiodidx(windOff.angles, sampling, windOff.trueFreq); + windOnPeriodIdx = getmeanperiodidx(windOn.angles, sampling, windOn.trueFreq); + + windOffIdx = windOffPeriodIdx(1):windOffPeriodIdx(2); + windOnIdx = windOnPeriodIdx(1):windOnPeriodIdx(2); % 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; + minLen = min(length(windOffIdx), length(windOnIdx)); + windOffIdx = windOffIdx(1:minLen); + windOnIdx = windOnIdx(1:minLen); % 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]); + fzOff = windOff.forces(windOffIdx, 3); + anglesOff = windOff.angles(windOffIdx); + anglesOn = windOn.angles(windOnIdx); + fzOn = windOn.forces(windOnIdx, 3); + meanDiff = fzOn - fzOff; + + time = (0:length(meanDiff) - 1) / nPeriodIdx; + [~, donwStrokeEndIdx] = min(anglesOff); + + figure; + setcolormap(); + hold on; + patchX = [0, donwStrokeEndIdx / nPeriodIdx, donwStrokeEndIdx / nPeriodIdx, 0]; + patchY = [-5, -5, 10, 10]; + + patch(patchX, patchY, [186, 216, 221] / 255, 'EdgeColor', 'none', 'FaceAlpha', 0.5); + h1 = plot(time, fzOff, 'LineWidth', 1.2); + h2 = plot(time, fzOn, 'LineWidth', 1.2); + h3 = plot(time, meanDiff, 'LineWidth', 1.2); + text((donwStrokeEndIdx / nPeriodIdx) / 2, -1.25, ... + 'Downstroke', ... + 'HorizontalAlignment', 'center'); + text(donwStrokeEndIdx / nPeriodIdx + ... + (1 - (donwStrokeEndIdx / nPeriodIdx)) / 2, -1.25, ... + 'Upstroke', ... + 'HorizontalAlignment', 'center'); + setgca(); + legend([h1, h2, h3], {'\airspeed = 0 m/s', '\airspeed = 7.7 m/s', 'Difference'}); + legend('boxoff'); + xlim([0, 1]); + ylim([-1.5 1.5]); + + if isinteresting(ResData(i, j), SAVE_TIKZ_PARAM) + figdir = 'figures/results/'; + filename = sprintf('inertia-period-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 + title(expParam); + + figure; + hold on; + plot(time, fzOff); + plot(time, fzOn); + plot(time, meanDiff); + plot(time(round(end / 5)), fzOn(round(end / 5)), 'or'); + yyaxis right; + plot(time, windOff.angles(windOffIdx)); + title(expParam); + legend('Wind off', 'Wind on', 'Diff'); + grid on; + xlim([0, 1]); + + figure; + ax1 = subplot(3, 1, 1); + setcolormap(); + hold on; + + plot(anglesOn, fzOn); + plot(anglesOn(round(end / 5)), fzOn(round(end / 5)), 'or'); + setgca('XTickLabel', [], ... + 'YLabel', [], ... + 'XLabel', []); + + ax2 = subplot(3, 1, 2); + setcolormap(); + hold on; + plot(anglesOff, fzOff); + % plot(meanPeakAngle(round(end / 5)), meanPeakFoff(round(end / 5)), 'or'); + setgca('XTickLabel', [], ... + 'XLabel', []); + + ax3 = subplot(3, 1, 3); + setcolormap(); + hold on; + plot(anglesOff, 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 % + + figure; + hold on; + plot(time, fzOff); + plot(time, fzOn); + plot(time, meanDiff); + plot(time(round(end / 5)), fzOn(round(end / 5)), 'or'); + yyaxis right; + plot(time, windOff.angles(windOffIdx)); + title(expParam); + legend('Wind off', 'Wind on', 'Diff'); + grid on; + xlim([0, 1]); figure; ax1 = subplot(3, 1, 1); setcolormap(); hold on; - plot(meanPeakAngle, meanPeakFon); - % plot(meanPeakAngle(round(end / 5)), meanPeakFon(round(end / 5)), 'or'); + plot(anglesOn, fzOn); + plot(anglesOn(round(end / 5)), fzOn(round(end / 5)), 'or'); setgca('XTickLabel', [], ... 'YLabel', [], ... 'XLabel', []); @@ -85,7 +173,7 @@ function plotwindonoff(ResData) ax2 = subplot(3, 1, 2); setcolormap(); hold on; - plot(meanPeakAngle, meanPeakFoff); + plot(anglesOff, fzOff); % plot(meanPeakAngle(round(end / 5)), meanPeakFoff(round(end / 5)), 'or'); setgca('XTickLabel', [], ... 'XLabel', []); @@ -93,7 +181,7 @@ function plotwindonoff(ResData) ax3 = subplot(3, 1, 3); setcolormap(); hold on; - plot(meanPeakAngle, meanDiff); + plot(anglesOff, meanDiff); % plot(meanPeakAngle(round(end / 5)), meanDiff(round(end / 5)), 'or'); linkaxes([ax1 ax2 ax3], 'xy'); ax1.YLim = 1.25 * ax1.YLim; @@ -111,6 +199,12 @@ function plotwindonoff(ResData) title(ax1, ['Wind on - ', expParam]); title(ax2, ['Wind off - ', expParam]); title(ax3, ['Diff - ', expParam]); + + % Add titles after having saved + title(ax1, ['Wind on - ', expParam]); + title(ax2, ['Wind off - ', expParam]); + title(ax3, ['Diff - ', expParam]); + end end @@ -118,8 +212,8 @@ function plotwindonoff(ResData) end -function meanPeak = averagepeaks(signal, sampling, freq) - % AVERAGEPEAKS Make the average of all peaks in the signal +function meanPeriod = getmeanperiodidx(signal, sampling, freq) + % GETMEANPEAK Return bounds of a period somewhere in the signal period = 1 / freq; nPeriodIdx = period * sampling; @@ -128,19 +222,14 @@ function meanPeak = averagepeaks(signal, sampling, freq) % 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); + meanPeriod = [peakLocs(round(end / 3)), peakLocs(round(end / 3)) + nPeriodIdx]; end function setgca(varargin) - xlabel('Flapping angle [deg]'); - ylabel('\cfz [N]'); + xlabel('t/T [-]'); + ylabel('\fz [N]'); set(gca, ... 'Box', 'off', ...