plotData.m

function figs = plotData(data, opts)
% Make all figures for IMR experimental data (dimensional, nondim, envelopes, CI, strain-rate, FFT, spectrogram, cloud).
%   figs = imr.plotData(data, Name=Value, ...)
%
% INPUT: data struct with any of these (the more you provide, the more plots work):
%   % REQUIRED for most plots
%   data.t_nd        % [N x T] nondimensional time (aka tmatrix)
%   data.R_nd        % [N x T] nondimensional radius (aka Rmatrix)
%
%   % OPTIONAL (recommended)
%   data.t_d         % [N x T] dimensional time
%   data.R_d         % [N x T] dimensional radius
%   data.Rdot_nd     % [N x T] nondim Rdot (aka Rdotmatrix)
%   data.strainrate  % [N x T] = (-2*Rdot_nd./R_nd)
%   data.Rmax_each   % [1 x T] per-trial Rmax (aka Rmax_range)
%   data.tc_each     % [1 x T] per-trial tc (aka tc)
%   data.Req_each    % [1 x T] per-trial equilibrium radius (dimensional)
%   data.noisyR      % noisyR(i,end).tvR = [t, R] for cloud plot
%
% NAME-VALUE OPTS (all false except the common ones):
%   DataDimPlot      (default: true)
%   DataNondimPlot   (default: true)
%   EnvelopePlot     (default: false)
%   ConfidencePlot   (default: false)
%   StrainRatePlot   (default: true)
%   FFTPlot          (default: false)
%   SpectroPlot      (default: false)
%   CloudPlot        (default: false)  % requires data.noisyR
%   LineWidth        (default: 2)
%   XMaxND           (default: 10)
%
% OUTPUT:
%   figs : struct of figure handles (only the ones created)
%
% EXAMPLE (from your prep script):
%   data = struct('t_nd', tmatrix, 'R_nd', Rmatrix, ...
%                 't_d', t_d, 'R_d', R_d, ...
%                 'Rdot_nd', Rdotmatrix, 'strainrate', (-2*Rdotmatrix./Rmatrix), ...
%                 'Rmax_each', Rmax_range, 'tc_each', tc);
%   figs = imr.plotData(data, DataDimPlot=true, DataNondimPlot=true, StrainRatePlot=true);

arguments
    data struct
    opts.DataDimPlot   (1,1) logical = true
    opts.DataNondimPlot(1,1) logical = true
    opts.VelocityPlot  (1,1) logical = false
    opts.EnvelopePlot  (1,1) logical = false
    opts.ConfidencePlot(1,1) logical = false
    opts.StrainRatePlot(1,1) logical = true
    opts.FFTPlot       (1,1) logical = false
    opts.SpectroPlot   (1,1) logical = false
    opts.CloudPlot     (1,1) logical = false
    opts.LineWidth     (1,1) double  {mustBePositive} = 2
    opts.XMaxND        (1,1) double  {mustBePositive} = 10
    opts.StrainPlot    (1,1) logical = true
    opts.StrainPhasePlot(1,1) logical = true
end

figs = struct();

has = @(f) isfield(data, f) && ~isempty(data.(f));
need = @(f) assert(has(f), "plotData:MissingField", "Required field '%s' not found in data.", f);

% Required baseline
need('t_nd'); need('R_nd');
t_nd = data.t_nd; R_nd = data.R_nd;

% Optional fields (use if present)
t_d = []; R_d = []; Rdot_nd = []; strain_rate = [];
Rmax_each = []; tc_each = []; Req_each = []; noisyR = [];
if has('t_d'),          t_d = data.t_d;                 end
if has('R_d'),          R_d = data.R_d;                 end
if has('Rdot_nd'),      Rdot_nd = data.Rdot_nd;         end
if has('strain_rate'),  strain_rate = data.strain_rate; end
if has('Rmax_each'),    Rmax_each = data.Rmax_each;     end
if has('tc_each'),      tc_each   = data.tc_each;       end
if has('Req_each'),     Req_each  = data.Req_each;      end
if has('noisyR'),       noisyR    = data.noisyR;        end

strain = []; strain_H = []; strain_proxy = [];
if has('strain'),       strain   = data.strain;   end
if has('strain_H'),        strain_H        = data.strain_H;        end
if has('strain_proxy'), strain_proxy = data.strain_proxy; end


% Build a shared time base for envelope/CI/cloud
t_end    = min(max(t_nd));
t_common = linspace(0, t_end, 1000)';
R_interp = nan(numel(t_common), size(R_nd,2));
for i = 1:size(R_nd,2)
    R_interp(:,i) = interp1(t_nd(:,i), R_nd(:,i), t_common, 'linear', NaN);
end

lw = opts.LineWidth; xmax = opts.XMaxND;

% ---- Dimensional radius ----
if opts.DataDimPlot && ~isempty(t_d) && ~isempty(R_d)
    figs.dataDim = figure('Name','Dimensional R(t)'); hold on
    for i = 1:size(R_d, 2)
        plot(t_d(:,i), R_d(:,i), 'k', 'LineWidth', lw)
    end
    xlabel('$t$','Interpreter','latex'); ylabel('$R$','Interpreter','latex');
    xlim([0, max(t_d(end,:))]); formatPlot();
end

% ---- Nondimensional radius ----
if opts.DataNondimPlot
    figs.dataND = figure('Name','Nondimensional R/Rmax vs t/tc'); hold on
    for i = 1:size(R_nd,2)
        plot(t_nd(:,i), R_nd(:,i), 'k', 'LineWidth', lw)
    end
    xlabel('$t^*$','Interpreter','latex'); ylabel('$R^*$','Interpreter','latex');
    xlim([0, xmax]); ylim([0, 1.25]); formatPlot();
end

% ---- Nondimensional wall velocity: Rdot_nd vs t/tc ----
if opts.VelocityPlot
    % Ensure Rdot_nd exists; if not, compute it from R_nd and t_nd
    if isempty(Rdot_nd)
        Rdot_nd = zeros(size(R_nd));
        for i = 1:size(R_nd,2)
            Rdot_nd(:,i) = gradient(R_nd(:,i), t_nd(:,i)); % d(R/Rmax)/d(t/tc)
        end
    end

    figs.velocityND = figure('Name','Nondimensional Rdot'); hold on
    for i = 1:size(Rdot_nd,2)
        plot(t_nd(:,i), Rdot_nd(:,i), 'k', 'LineWidth', lw);
    end
    xlabel('$t^*$','Interpreter','latex');
    ylabel('$\dot{R}^*$','Interpreter','latex');
    xlim([0, xmax]);  % use nondimensional window
    formatPlot();
end


% ---- Envelope (min/max/mean) ----
if opts.EnvelopePlot
    figs.envelope = figure('Name','Envelope'); hold on
    Rmin  = nanmin(R_interp, [], 2);
    Rmax  = nanmax(R_interp, [], 2);
    Rmean = nanmean(R_interp, 2);
    fill([t_common; flipud(t_common)], [Rmin; flipud(Rmax)], [0.678, 0.847, 0.902], ...
         'FaceAlpha', 1.0, 'EdgeColor', 'none');
    plot(t_common, Rmin, 'k-', 'LineWidth', lw)
    plot(t_common, Rmax, 'k-', 'LineWidth', lw)
    plot(t_common, Rmean,'r-', 'LineWidth', lw)
    legend({'Experimental Window','Min/Max','', 'Mean'}, 'Interpreter','latex', 'Location','northeast');
    xlabel('$t/t_\mathrm{c}$','Interpreter','latex'); ylabel('$R/R_{\mathrm{max}}$','Interpreter','latex');
    xlim([0, xmax]); ylim([0, 1.25]); formatPlot();
end

% ---- 95% CI ----
if opts.ConfidencePlot
    figs.ci = figure('Name','95% CI'); hold on
    Rmean = nanmean(R_interp, 2);
    Rstd  = nanstd(R_interp, 0, 2);
    ci_upper = Rmean + 1.96*Rstd; ci_lower = Rmean - 1.96*Rstd;
    Rmin = nanmin(R_interp, [], 2); Rmax = nanmax(R_interp, [], 2);
    fill([t_common; flipud(t_common)], [ci_lower; flipud(ci_upper)], [1, 0.6, 0.6], ...
         'FaceAlpha', 1.0, 'EdgeColor', 'none');
    plot(t_common, Rmin, 'k-', 'LineWidth', lw)
    plot(t_common, Rmax, 'k-', 'LineWidth', lw)
    plot(t_common, Rmean,'r-', 'LineWidth', lw)
    legend({'95\% CI','Min/Max','', 'Mean'}, 'Interpreter','latex', 'Location','northeast');
    xlabel('$t/t_c$','Interpreter','latex'); ylabel('$R/R_{\mathrm{max}}$','Interpreter','latex');
    xlim([0, xmax]); ylim([0, 1.25]); formatPlot();
end

% ---- Cloud of noisy realizations (requires data.noisyR) ----
if opts.CloudPlot
    figure; hold on
    set(gcf,'Renderer','opengl');  % ensure FaceAlpha works

    % Interpolate noisy realizations onto t_common
    nRealizations = size(noisyR, 1);
    Rvals_interp  = nan(length(t_common), nRealizations);
    for i = 1:nRealizations
        t_i = noisyR(i,end).tvR(:,1);
        R_i = noisyR(i,end).tvR(:,2);
        [t_i, uniq] = unique(t_i(:),'stable');  % ensure monotonic
        R_i = R_i(uniq);
        Rvals_interp(:, i) = interp1(t_i, R_i, t_common, 'linear', NaN);
    end

    % Compute bounds and mean
    Rmin  = nanmin(Rvals_interp, [], 2);
    Rmax  = nanmax(Rvals_interp, [], 2);
    Rmean = nanmean(Rvals_interp, 2);

    % --- FILL in cyan between min and max ---
    good = isfinite(Rmin) & isfinite(Rmax);
    d = diff([false; good; false]);
    s = find(d==1); e = find(d==-1)-1;
    for k = 1:numel(s)
        g = s(k):e(k);
        fill([t_common(g); flipud(t_common(g))], ...
             [Rmin(g);   flipud(Rmax(g))], ...
             'c', 'FaceAlpha', 0.75, 'EdgeColor', 'none'); % cyan
    end

    % Overlay lines
    plot(t_common, Rmin,  'k-', 'LineWidth', 2);
    plot(t_common, Rmax,  'k-', 'LineWidth', 2);
    plot(t_common, Rmean, 'r-', 'LineWidth', 3);

    legend({'Noisy Envelope','Min/Max','', 'Noisy Mean'}, ...
           'Interpreter','latex', 'Location','northeast');
    xlabel('$t/t_\mathrm{c}$','Interpreter','latex');
    ylabel('$R/R_\mathrm{max}$','Interpreter','latex');
    xlim([0, 10]); ylim([0, 1.25]);
    formatPlot();
end

% ---- Strain-rate (dimensional, get rid of tc_each scaling to go ND) ----
if opts.StrainRatePlot
    % compute from Rdot if not provided
    if isempty(strain_rate) && ~isempty(Rdot_nd)
        strain_rate = (-2 .* Rdot_nd) ./ R_nd;
        strain_rate(~isfinite(strain_rate)) = NaN;
    end
    needSR = ~isempty(strain_rate);
    assert(needSR, 'StrainRatePlot requires data.strainrate or data.Rdot_nd.');
    figs.strain = figure('Name','Strain Rate'); hold on
    for i = 1:size(R_nd,2)
        plot(t_nd(:,i), strain_rate(:,i)./tc_each(i), 'k', 'LineWidth', lw);
    end
    xlabel('$t/t_c$','Interpreter','latex'); ylabel('$\dot{\varepsilon} = -2 \dot{R} / R$','Interpreter','latex');
    xlim([0, xmax]); formatPlot();
end

% ---- Strain (magnitude) plot: |log strain| vs t/tc ----
if opts.StrainPlot
    % Prefer true Hencky if available; fall back to equilibrium log strain.
    if ~isempty(strain_H)
        strain_to_plot = abs(strain_H);
        ylab = '$|\varepsilon_H| = |\ln(R/R_0)|$';
    elseif ~isempty(strain)
        strain_to_plot = abs(strain);
        ylab = '$|\ln(R/R_{\mathrm{eq}})|$';
    else
        % Final fallback: proxy (ND)
        assert(~isempty(strain_proxy), 'StrainPlot needs strain_H, strain, or strain_proxy.');
        strain_to_plot = strain_proxy;   % already magnitude-like
        ylab = '$|R/R_{\mathrm{max}} - (R/R_{\mathrm{max}})_{\mathrm{tail}}|$';
    end

    figs.strainMag = figure('Name','Strain Magnitude'); hold on
    for i = 1:size(R_nd,2)
        plot(t_nd(:,i), strain_to_plot(:,i), 'k', 'LineWidth', lw);
    end
    xlabel('$t/t_c$','Interpreter','latex'); ylabel(ylab,'Interpreter','latex');
    xlim([0, xmax]); formatPlot();
end

% ---- Strain–Strain-Rate Phase Space: ε (x) vs ė (y) ----
if opts.StrainPhasePlot
    % Ensure we have strain rate (dimensionless) available or computable
    if isempty(strain_rate)
        if ~isempty(Rdot_nd)
            strain_rate = (-2 .* Rdot_nd) ./ R_nd;
            strain_rate(~isfinite(strain_rate)) = NaN;
        else
            error('StrainPhasePlot requires data.strainRate (or strain_rate) or Rdot_nd to compute it.');
        end
    end

    % Choose which strain to plot on x-axis (preference: Hencky, then equilibrium-log-strain, then proxy)
    if ~isempty(strain_H)
        strain_to_plot = strain_H;                  % true Hencky strain ln(R/R0)
        xlab = '$\varepsilon_H = \ln(R/R_0)$';
    elseif ~isempty(strain)
        strain_to_plot = strain;                    % equilibrium-referenced ln(R/Req)
        xlab = '$\ln\!\big(R/R_{\mathrm{eq}}\big)$';
    else
        assert(~isempty(strain_proxy), 'StrainPhasePlot needs strain_H, strain, or strain_proxy.');
        strain_to_plot = strain_proxy;              % fallback (dimensionless strain-like proxy)
        xlab = '$\text{strain proxy}$';
    end

    figs.strainPhase = figure('Name','Phase: Strain vs Strain-Rate'); hold on
    for i = 1:size(R_nd,2)
        x = strain_to_plot(:,i);
        y = strain_rate(:,i);
        good = isfinite(x) & isfinite(y);
        if nnz(good) < 2, continue; end

        % Plot trajectory
        plot(x(good), y(good), 'k-', 'LineWidth', lw);

        % Optional: show start and end markers for direction
        % ig = find(good);
        % plot(x(ig(1)),  y(ig(1)), 'k','LineWidth', 1.5);
        % plot(x(ig(end)),y(ig(end)), 'k', 'LineWidth', 1.5);
    end
    xlabel('$\varepsilon^*$', 'Interpreter','latex','FontSize',36);
    ylabel('$\dot{\varepsilon}^*$','Interpreter','latex');
    % Try to keep a sensible box around typical ranges; feel free to adjust:
    xlim auto; ylim auto;
    formatPlot();
end


% ---- FFT of strain rate per trial ----
if opts.FFTPlot
    assert(~isempty(strain_rate) || ~isempty(Rdot_nd), 'FFTPlot needs strainRate or Rdot.');
    if isempty(strain_rate)
        strain_rate = (-2 .* Rdot_nd) ./ R_nd;
        strain_rate(~isfinite(strain_rate)) = NaN;
    end
    figs.fft = figure('Name','FFT of Strain Rate'); hold on
    for i = 1:size(R_nd,2)
        s = strain_rate(:,i); s(~isfinite(s)) = 0;
        L = numel(s); NFFT = 2^nextpow2(L);
        Y = fft(s, NFFT);
        f = (0:NFFT-1)'/NFFT; % normalized frequency
        plot(f, abs(Y), 'LineWidth', 1.5);
    end
    xlabel('Normalized $f$','Interpreter','latex'); ylabel('$|\mathrm{FFT}|$','Interpreter','latex');
    formatPlot();
end

% ---- “Spectrogram” across stretch ratios (FFT magnitude vs freq & stretch) ----
if opts.SpectroPlot
    % stretch = Rmax/Req (dimensional). Estimate Req if missing.
    if isempty(Req_each)
        assert(~isempty(R_d), 'SpectroPlot needs data.Req_each or data.R_d.');
        last10 = ceil(0.10 * size(R_d,1));
        Req_each = mean(R_d(end-last10+1:end, :), 1);
    end
    assert(~isempty(Rmax_each), 'SpectroPlot needs data.Rmax_each.');
    if isempty(strain_rate) && ~isempty(Rdot_nd)
        strain_rate = (-2 .* Rdot_nd) ./ R_nd;
        strain_rate(~isfinite(strain_rate)) = NaN;
    end
    assert(~isempty(strain_rate), 'SpectroPlot needs strainRate.');

    stretch = Rmax_each(:) ./ Req_each(:);
    L = size(strain_rate,1);
    NFFT = 2^nextpow2(L);
    amp = zeros(size(R_nd,2), NFFT/2);
    for i = 1:size(R_nd,2)
        s = strain_rate(:,i); s(~isfinite(s)) = 0;
        Y = fft(s, NFFT);
        amp(i,:) = abs(Y(1:NFFT/2));
    end
    [stretchSorted, si] = sort(stretch);
    figs.spectro = figure('Name','Stretch vs FFT amplitude');
    imagesc( (0:NFFT/2-1), stretchSorted, log10(amp(si,:)+1e-6) );
    axis xy; colormap winter; cb = colorbar;
    cb.Label.String = 'log$_{10}$(Amplitude)'; cb.Label.Interpreter = 'latex'; cb.TickLabelInterpreter = 'latex';
    xlabel('FFT bin','Interpreter','latex'); ylabel('$R_{\max}/R_{\mathrm{eq}}$','Interpreter','latex');
    set(gca,'TickLabelInterpreter','latex'); formatPlot();
end
end

% ---- Local style helper ----
function formatPlot()
    ax  = gca;
    fig = gcf;

    % ---------- Make window as large as possible with fixed aspect ratio ----------
    ar = 1.5;                                 % desired width/height (match pbaspect)
    scr = get(0,'ScreenSize');                % [left bottom width height] (px)
    pad = 60;                                 % margin for OS chrome

    scrW = scr(3) - 2*pad;
    scrH = scr(4) - 2*pad;

    W = scrW; H = W/ar;                       % try width-limited
    if H > scrH                                % if too tall, limit by height instead
        H = scrH; W = H*ar;
    end

    left   = pad + (scr(3) - 2*pad - W)/2;    % center on screen
    bottom = pad + (scr(4) - 2*pad - H)/2;

    set(fig, 'Units','pixels', 'Position',[left bottom W H]);

    % ---------- Axes styling ----------
    box(ax,'on');
    set(ax,'LineWidth',2,'FontSize',36,'TickLabelInterpreter','latex');
    pbaspect(ax,[ar 1 1]);                    % keep same visual ratio as window

    % ---------- Tighten horizontally only (keep vertical unchanged) ----------
    drawnow;                                  % update layout metrics
    set(ax,'Units','normalized');

    % Avoid extra padding when exporting
    set(ax,'LooseInset', get(ax,'TightInset'));

    ti  = get(ax,'TightInset');               % [left bottom right top] (norm)
    pos = get(ax,'Position');                 % [x y w h] (norm)

    % keep current y, h; adjust x and w using left/right insets
    pos(1) = ti(1);                           % left edge clears y labels
    pos(3) = 1 - ti(1) - ti(3);               % fill remaining width
    set(ax,'Position',pos);

    % ---------- Clean vector export defaults ----------
    set(fig,'Renderer','painters','Color','w');
    set(ax,'Color','w');                   % transparent axes bg
end