untitled2.m

clear;clc;close all;
expDataPrep
t = tmatrix(:,1);
R = Rmatrix(:,1);
Rdot = Rdotmatrix(:,1);
animate_strainrate_live(R, Rdot, t, ...
    'UseND', true, 'NDHalfWidth', 1.5, ...
    'Nx', 600, 'Speedup', 5, ...
    'LockColorbar', false, ...
    'TimeScale', 1, 'TimeUnit', 'ND', 'PaceScale', 1);

function animate_strainrate_live(R, Rdot, t, varargin)
% animate_strainrate_live (LaTeX everywhere)
% eps_eq(r,t) = 2*sqrt(3)*|R(t)^2 * Rdot(t)| / r^3   for r >= R(t)
%
% New (time handling):
%   TimeScale, TimeUnit, TimeOffset, PaceScale (display vs pacing decoupled)

% ---------- options ----------
p = inputParser;
p.addParameter('UseND', true);
p.addParameter('NDHalfWidth', 1.5);
p.addParameter('SpanInRmax', 8);
p.addParameter('Nx', 500);
p.addParameter('Speedup', 10);
p.addParameter('SaveAs', '');
p.addParameter('UseLog', true);
p.addParameter('LogFloor', 1e-20);
p.addParameter('LockColorbar', true);
p.addParameter('FixedCLim', []);
p.addParameter('GlobalPercentiles', [5 99.5]);

% explicit time handling
p.addParameter('TimeScale', 1);     % DISPLAY: t_display = t*TimeScale + TimeOffset
p.addParameter('TimeUnit', '');     % e.g., 's', '\mu s', 'ND' (LaTeX-ready or plain)
p.addParameter('TimeOffset', 0);
p.addParameter('PaceScale', 1);     % PACING: pause uses Δt * PaceScale / Speedup

p.parse(varargin{:});

UseND    = p.Results.UseND;
NDL      = p.Results.NDHalfWidth;
span     = p.Results.SpanInRmax;
Nx       = p.Results.Nx;
spd      = max(p.Results.Speedup, eps);
base     = p.Results.SaveAs;
UseLog   = p.Results.UseLog;
LogFloor = p.Results.LogFloor;
LockCB   = p.Results.LockColorbar;
FixedCL  = p.Results.FixedCLim;
globPr   = p.Results.GlobalPercentiles;

TimeScale  = p.Results.TimeScale;
TimeUnit   = p.Results.TimeUnit;
TimeOffset = p.Results.TimeOffset;
PaceScale  = p.Results.PaceScale;

% ---------- data prep ----------
R    = R(:).'; 
Rdot = Rdot(:).'; 
t    = t(:).';

% field definition
eps_eq_fun = @(Rp,Rdp,r) 2*sqrt(3)*abs(Rp.^2 .* Rdp) ./ (r.^3 + eps);

% choose domain (ND or physical)
if UseND
    L = NDL;
    x = linspace(-L, L, Nx); 
    y = x;
    toUnits = 1;
    xLabel = '$x$ (ND)'; 
    yLabel = '$y$ (ND)';
else
    Rmax = max(R);
    L = span*Rmax;             % meters
    x = linspace(-L, L, Nx); 
    y = x;
    toUnits = 1e6;             % display in micrometers
    xLabel = '$x~(\mu\mathrm{m})$'; 
    yLabel = '$y~(\mu\mathrm{m})$';
end
[XX,YY] = meshgrid(x, y);
RR = hypot(XX, YY);

% ---------- one-time global CLim if locking and not provided ----------
if LockCB && isempty(FixedCL)
    valsAll = [];
    nR = 800;
    for j = 1:numel(t)
        Rp  = R(j); 
        Rdp = Rdot(j);
        rmin = max(Rp*1.001, L/Nx);  % just outside bubble
        if rmin >= L, continue; end
        r = logspace(log10(max(rmin,1e-12)), log10(L), nR);
        Eline = eps_eq_fun(Rp, Rdp, r);
        v = UseLog * log10(max(Eline, LogFloor)) + (~UseLog) * Eline;
        v(~isfinite(v)) = [];
        valsAll = [valsAll, v]; %#ok<AGROW>
    end
    if isempty(valsAll)
        FixedCL = UseLog * [-10 -9] + (~UseLog) * [0 1];
    else
        pr = prctile(valsAll, globPr);
        if ~(isfinite(pr(1)) && isfinite(pr(2)) && pr(2) > pr(1))
            mu = mean(valsAll, 'omitnan'); 
            sd = std(valsAll, 0, 'omitnan'); 
            if ~isfinite(sd) || sd == 0, sd = 0.25; end
            pr = [mu-2*sd, mu+2*sd];
            if pr(2) <= pr(1), pr(2) = pr(1) + 1e-3; end
        end
        FixedCL = pr;
    end
end

% ---------- figure & handles ----------
fig = figure('Color','w','Units','pixels','Position',[100 100 760 760]);
ax  = axes(fig); 
colormap(ax, parula);

% LaTeX everywhere
set(ax, 'TickLabelInterpreter','latex');
img = imagesc(ax, x*toUnits, y*toUnits, nan(size(RR)));
set(ax, 'YDir', 'normal'); axis(ax, 'image');
ax.XAxis.Exponent = 0; ax.YAxis.Exponent = 0;
xtickformat(ax, '%.1f'); ytickformat(ax, '%.1f');
xlabel(ax, xLabel, 'Interpreter','latex'); 
ylabel(ax, yLabel, 'Interpreter','latex');
xlim(ax, [-L L]*toUnits); ylim(ax, [-L L]*toUnits);

cb = colorbar(ax);
cb.TickLabelInterpreter = 'latex';
cb.Label.Interpreter = 'latex';
if UseLog
    cb.Label.String = '$\log_{10}(\dot{\epsilon}_{\mathrm{eq}})$';
else
    cb.Label.String = '$\dot{\epsilon}_{\mathrm{eq}}~(\mathrm{s}^{-1})$';
end
if LockCB && ~isempty(FixedCL), caxis(ax, FixedCL); end

hold(ax, 'on');
th = linspace(0, 2*pi, 400);
circ = plot(ax, nan(size(th)), nan(size(th)), 'k-', 'LineWidth', 1.1);
hold(ax, 'off');

% ---------- optional video ----------
vw = []; opened = false; framesWritten = 0;
if ~isempty(base)
    prof = chooseProfile();
    vw = VideoWriter([base profileExt(prof)], prof);
    vw.FrameRate = 30;
end

% ---------- main loop ----------
for j = 1:numel(t)
    Rp  = R(j); 
    Rdp = Rdot(j);

    % field on grid and mask inside bubble
    E = eps_eq_fun(Rp, Rdp, RR);
    E(RR < Rp) = NaN;
    V = UseLog * log10(max(E, LogFloor)) + (~UseLog) * E;
    V(~isfinite(V)) = NaN;

    % update image and circle
    set(img,  'CData', V, 'XData', x*toUnits, 'YData', y*toUnits);
    set(circ, 'XData', (Rp*cos(th))*toUnits, 'YData', (Rp*sin(th))*toUnits);

    % color limits
    if LockCB && ~isempty(FixedCL)
        caxis(ax, FixedCL);
    else
        finiteVals = V(isfinite(V));
        [lo, hi] = safeCLim(finiteVals, UseLog);
        caxis(ax, [lo hi]);
    end

    % bounds and title (LaTeX)
    xlim(ax, [-L L]*toUnits); 
    ylim(ax, [-L L]*toUnits);

    t_disp = t(j)*TimeScale + TimeOffset;
    if isempty(TimeUnit)
        ttlTime = sprintf('$t = %.6g$', t_disp);
    else
        if contains(TimeUnit, '\') % already LaTeX-coded (e.g., '\mu s')
            ttlTime = sprintf('$t = %.6g\\,%s$', t_disp, TimeUnit);
        else                       % plain text -> roman
            ttlTime = sprintf('$t = %.6g\\,\\mathrm{%s}$', t_disp, TimeUnit);
        end
    end
    if UseND
        ttlDom = sprintf('(ND domain = $\\pm%.1f$)', NDL);
    else
        ttlDom = sprintf('(domain = $\\pm%.1f~R_{\\max}$)', span);
    end
    title(ax, sprintf('%s  %s', ttlTime, ttlDom), 'Interpreter','latex');

    drawnow;

    % write frame
    if ~isempty(vw)
        if ~opened, open(vw); opened = true; end
        writeVideo(vw, getframe(fig));
        framesWritten = framesWritten + 1;
    end

    % pacing (real Δt, scaled only if requested)
    if j < numel(t)
        dt = (t(j+1) - t(j)) * PaceScale;
        pause(max(dt/spd, 0));
    end
end

if ~isempty(vw) && opened
    close(vw);
    if framesWritten == 0
        warning('No video frames were written.');
    end
end
end

% ---------- helpers ----------
function [lo, hi] = safeCLim(vals, UseLog)
if isempty(vals)
    if UseLog, lo = -10; hi = -9; else, lo = 0; hi = 1; end
    return
end
lo = prctile(vals, 5, 'all');
hi = prctile(vals, 99.5, 'all');
if ~isfinite(lo) || ~isfinite(hi) || hi <= lo
    mu = mean(vals, 'omitnan'); 
    sd = std(vals, 0, 'omitnan');
    if ~isfinite(sd) || sd == 0, sd = 0.25; end
    lo = mu - 2*sd; 
    hi = mu + 2*sd;
    if ~isfinite(lo) || ~isfinite(hi) || hi <= lo
        lo = mu - 1; 
        hi = mu + 1;
        if hi <= lo, hi = lo + 1e-3; end
    end
end
end

function prof = chooseProfile()
try
    names = string({VideoWriter.getProfiles.Name});
    if any(names == "MPEG-4"), prof = 'MPEG-4';
    else, prof = 'Motion JPEG AVI';
    end
catch
    prof = 'Motion JPEG AVI';
end
end

function ext = profileExt(p)
if strcmp(p,'MPEG-4'), ext = '.mp4'; else, ext = '.avi'; end
end