function plot_validity_window()
% =========================================================================
% TAYLOR-ARIS VALIDITY WINDOW (analytical; contains no simulation data)
% Plots the dispersive enhancement kappa_eff/kappa = 1 + Pe(r)^2/210,
% with Pe(r) = H*r*h/kappa, and marks the critical radius
% r_TA = sqrt(210)*kappa/(H*h) for four tracer scenarios.
% Reproduces the validity-window figure of Section 4.
% =========================================================================

r = logspace(-5.2, log10(0.05), 400);      % radial distance (m)

% Scenarios correspond to rows of the design-space table in the text.
K   = [1e-6,   1.4e-7, 4.25e-10, 1.4e-7];  % diffusivities (m^2/s)
Hs  = [0.5,    0.5,    0.5,      0.5];     % expansion rates (1/s)
hs  = [1e-3,   1e-3,   1e-3,     1e-4];    % gap thicknesses (m)
cols = {[0 0.45 0.74], [0.20 0.63 0.17], [0.84 0.15 0.16], [0.20 0.63 0.17]};
sty  = {'-', '-', '-', '--'};
lab  = {'Simulation scalar: $\kappa=10^{-6}$, $h=1$ mm', ...
        'Thermal, $h=1$ mm', ...
        'Fluorescein, $h=1$ mm', ...
        'Thermal, $h=100\,\mu$m (design point)'};

figure('Color', 'w'); hold on;
for i = 1:4
    keff = 1 + (Hs(i) .* r .* hs(i) ./ K(i)).^2 / 210;
    plot(r*1000, keff, sty{i}, 'Color', cols{i}, 'LineWidth', 2);
end

% keff = 2*kappa criterion (defines r_TA) and the cell rim
plot([6e-3, 80], [2, 2], 'k--', 'LineWidth', 1.2);
plot([50, 50], [0.8, 3e6], ':', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);

for i = 1:4
    r_TA = sqrt(210) * K(i) / (Hs(i) * hs(i));
    if r_TA < 0.05
        plot(r_TA*1000, 2, 'o', 'MarkerFaceColor', cols{i}, ...
             'MarkerEdgeColor', 'k', 'MarkerSize', 7);
    end
    fprintf('r_TA scenario %d: %.4g mm\n', i, r_TA*1000);
end

set(gca, 'XScale', 'log', 'YScale', 'log');
xlim([6e-3, 80]); ylim([0.8, 3e6]); grid on;
xlabel('Radial distance, $r$ (mm)', 'Interpreter', 'latex');
ylabel('Dispersive enhancement, $\kappa_{\mathrm{eff}}/\kappa$', ...
       'Interpreter', 'latex');
title('Taylor--Aris validity window: $\kappa_{\mathrm{eff}}/\kappa = 1 + Pe(r)^2/210$', ...
      'Interpreter', 'latex');
legend(lab, 'Interpreter', 'latex', 'Location', 'northwest');
end