function plot_reynolds_figure()
% =========================================================================
% FIGURE 6 - FINITE-REYNOLDS BEHAVIOUR OF THE ELLIPTICAL CELL
%
% (a) anisotropy ratio Hx/Hy versus Re_squeeze, five fixed-domain points
%     plus three moving-mesh (ALE) points, against the five-point
%     least-squares law  Hx/Hy = (be^2/ae^2)(1 + a1 Re + a2 Re^2).
% (b) relaxation of the transverse profile in the ALE run: peak location
%     migrating to the mid-plane, and the profile deformation NORMALISED
%     by the instantaneous Reynolds number.
%
% All data traceable to HeleShaw_CFD_results.xlsx, sheets 'Runs',
% 'Analysis', 'ALE_Transient'.
%
% REPLACES the earlier version of this figure, which plotted only three
% fixed-domain points, used the superseded two-point fit (0.151, -0.034),
% and carried two label artefacts.
%
% NOTE ON PANEL (a): Eq. (24) describes the GRID-CONVERGED solution, so the
% fit is drawn through grid-corrected values. Raw measurements are shown as
% open markers so the size of the discretization correction is visible.
%
% NOTE ON PANEL (b): the deformation is plotted as RMS/Re, not raw RMS.
% Re falls from 0.998 to 0.69 during the run as the gap closes, so raw RMS
% falls partly for that reason alone. Normalised, the deformation collapses
% by a factor of ~3.5 (the physical relaxation) and then plateaus near 1.8,
% which is the honest statement and is what Sect. 7.3 claims.
% =========================================================================

% ---- edit this if the equation renumbers -------------------------------
EQNUM = 26;   % current number of Eq.(Refit); check the compiled PDF

% ---- theory and fit -----------------------------------------------------
ratio_theory = 0.25;            % be^2/ae^2 for the 2:1 cell
a1 =  0.1454;                   % five-point least-squares, Analysis sheet
a2 = -0.0275;
grid_err_pct = 0.0692;          % Richardson coarse-grid error, in dev %

% ---- fixed-domain series (Runs sheet; Re=1 uses the fine mesh) ----------
Re_fd    = [0.005   0.25     0.5      1.0      2.0    ];
ratio_fd = [0.250395 0.259134 0.266890 0.279322 0.295431];

% grid-corrected values, which is what Eq. (EQNUM) describes
ratio_fd_corr = ratio_fd - ratio_theory*grid_err_pct/100;

% ---- moving-mesh (ALE) series (ALE_Transient sheet) ---------------------
Re_ale    = [0.99     0.89     0.6905  ];
ratio_ale = [0.279343 0.277753 0.272826];

% ---- ALE profile relaxation (ALE_Transient sheet) ----------------------
t_ale    = [0.001  0.002  0.003  0.004  0.005  0.055  0.15474];
peak_ale = [0.4620 0.4630 0.4639 0.4648 0.4658 0.4962 0.5020 ];
rms_ale  = [6.3166 6.0432 5.7856 5.5421 5.3116 1.5963 1.2715 ];
h_ale    = 0.001 - 0.002*t_ale;          % constant plate speed 2 mm/s
Re_t     = 0.002*h_ale/(2e-6);           % Re_squeeze = |hdot| h / 2 nu
rms_norm = rms_ale ./ Re_t;              % deformation per unit Re

RMS_G       = 1.3195;   % RMS of Paper B's sextic g, per unit Re
RMS_FD_ARTE = 6.9;      % fixed-domain velocity-inlet artefact, per unit Re

% ---- colours (matching the other figures in the paper) -----------------
cBlue = [0.00 0.45 0.74];
cRed  = [0.84 0.15 0.16];
cGold = [0.93 0.69 0.13];
cGrey = [0.50 0.50 0.50];
cGrn  = [0.16 0.58 0.24];

figure('Color','w','Position',[60 120 1250 460]);

% =========================== PANEL (a) ==================================
subplot(1,2,1); hold on; box on;

Refine = linspace(0, 2.15, 400);
fit    = ratio_theory*(1 + a1*Refine + a2*Refine.^2);

h1 = plot(Refine, fit, '-', 'Color', [0.2 0.2 0.2], 'LineWidth', 1.8);
h2 = plot([0 2.15], ratio_theory*[1 1], '--', 'Color', cRed, 'LineWidth', 1.5);
h3 = plot(Re_fd, ratio_fd_corr, 'o', 'MarkerSize', 8, ...
          'MarkerFaceColor', cBlue, 'MarkerEdgeColor', 'k');
h4 = plot(Re_fd, ratio_fd, 'o', 'MarkerSize', 8, ...
          'MarkerFaceColor', 'none', 'MarkerEdgeColor', cBlue, 'LineWidth', 1);
h5 = plot(Re_ale, ratio_ale, 's', 'MarkerSize', 9, ...
          'MarkerFaceColor', cGold, 'MarkerEdgeColor', 'k');

xlabel('$Re_{\mathrm{squeeze}}$', 'Interpreter', 'latex');
ylabel('$H_x/H_y$', 'Interpreter', 'latex');
title('(a) anisotropy ratio versus inertia', 'Interpreter', 'latex');
xlim([0 2.15]); ylim([0.245 0.300]); grid on;

legend([h1 h2 h3 h4 h5], ...
  {sprintf('least-squares fit, Eq. (%d)', EQNUM), ...
   'lubrication theory $b_e^2/a_e^2 = 0.25$', ...
   'CFD, fixed domain (grid-corrected)', ...
   'CFD, fixed domain (raw)', ...
   'CFD, moving mesh (ALE)'}, ...
  'Interpreter', 'latex', 'Location', 'northwest', 'FontSize', 8);

% =========================== PANEL (b) ==================================
subplot(1,2,2); box on;

yyaxis left
pk = plot(t_ale, peak_ale, 'o-', 'Color', cBlue, 'LineWidth', 1.6, ...
          'MarkerSize', 7, 'MarkerFaceColor', cBlue);
hold on
mp = plot([8e-4 0.2], [0.5 0.5], '--', 'Color', cRed, 'LineWidth', 1.5);
fd = plot(1e-3, 0.4611, 'v', 'MarkerSize', 9, 'MarkerFaceColor', cGrey, ...
          'MarkerEdgeColor', 'k');
ylabel('peak location $\zeta$', 'Interpreter', 'latex');
ylim([0.455 0.508]);
set(gca, 'YColor', cBlue);

yyaxis right
rn = plot(t_ale, rms_norm, '^--', 'Color', cGrn, 'LineWidth', 1.5, ...
          'MarkerSize', 7, 'MarkerFaceColor', cGrn);
hold on
gl = plot([8e-4 0.2], RMS_G*[1 1], ':', 'Color', cGrn, 'LineWidth', 1.8);
al = plot([8e-4 0.2], RMS_FD_ARTE*[1 1], '-.', 'Color', cGrey, 'LineWidth', 1.4);
ylabel('RMS departure per unit $Re$ (\%)', 'Interpreter', 'latex');
ylim([0 7.6]);
set(gca, 'YColor', cGrn);

set(gca, 'XScale', 'log');
xlim([8e-4 0.2]); grid on;
xlabel('$t$ (s)', 'Interpreter', 'latex');
title('(b) moving mesh restores mid-plane symmetry', 'Interpreter', 'latex');

legend([pk mp fd rn gl al], ...
  {'peak $\zeta$ (ALE) --- left axis', ...
   'mid-plane $\zeta = 1/2$', ...
   'fixed-domain steady value', ...
   'RMS$/Re$ (ALE) --- right axis', ...
   'first-order theory, $1.32$', ...
   'fixed-domain artefact, $\simeq 6.9$'}, ...
  'Interpreter', 'latex', 'Location', 'southoutside', ...
  'NumColumns', 2, 'FontSize', 8);

% ---- export ------------------------------------------------------------
exportgraphics(gcf, 'fig_reynolds.png', 'Resolution', 600);
% For the journal, also produce a vector copy:
% exportgraphics(gcf, 'Fig6.eps', 'ContentType', 'vector');

% ---- console cross-check against the manuscript ------------------------
fprintf('\n  Re      raw       grid-corr   Eq.(%d)     residual\n', EQNUM);
for i = 1:numel(Re_fd)
    f = ratio_theory*(1 + a1*Re_fd(i) + a2*Re_fd(i)^2);
    fprintf('  %-6.3f  %.6f  %.6f   %.6f   %+.4f %%\n', Re_fd(i), ...
        ratio_fd(i), ratio_fd_corr(i), f, ...
        100*(ratio_fd_corr(i)-f)/ratio_theory);
end
fprintf('\n  ALE profile relaxation:\n');
fprintf('   t (s)     Re      peak zeta   RMS %%    RMS/Re\n');
for i = 1:numel(t_ale)
    fprintf('  %7.5f  %6.4f   %.4f    %6.4f   %5.3f\n', ...
        t_ale(i), Re_t(i), peak_ale(i), rms_ale(i), rms_norm(i));
end
fprintf('\n  Paper B first-order value: %.4f %%/Re\n', RMS_G);
fprintf('  ALE plateau:               %.3f %%/Re  (%.2fx)\n', ...
        rms_norm(end), rms_norm(end)/RMS_G);
end