kv_uniqueness_from_stress_integrals.m

function kv_uniqueness_from_stress_integrals(opts)
% KV uniqueness maps from stress integrals (Table 1) using expDataPrep outputs.
% One figure: manifold over (mu, G) with uniqueness U_* = (balance)*(magnitude).
%
% Balance:  B = 1 - |S_N - S_H| / S_K            (rewards simultaneous viscous+elastic)
% Magnitude: M = S_K / max(S_K(:))               (rewards strong overall KV stress)
% Final:     U_* = B .* M
%
% Table 1 (non-dimensional stress integrals):
%   S*_v  = -(4/Re) * (Rdot*/R*)                          (Newtonian part)
%   S*_NH = (1/(2*Ca)) * [4*lambda^-1 + lambda^-4 - 5]    (Neo-Hookean part)
%   S*_KV = S*_v + S*_NH
% where
%   Re_k(mu) = Rmax_k * sqrt(rho*P_inf) / mu
%   Ca(G)    = P_inf / G
% and lambda = R*/Req*.
%
% Dependencies: expDataPrep must define Rmatrix, tmatrix, Rdotmatrix, R_d, tc (or tc_each),
% Req_each (or will be estimated), and Rmax_each or Rmax_range.

% ---------------- Options ----------------
if nargin < 1
    opts = struct();
end

if isfield(opts,'material') && ~isempty(opts.material)
    material = opts.material;
else
    material = 0;   % default: synthetic data
end

if isfield(opts,'mu_vec') && ~isempty(opts.mu_vec)
    mu_vec = opts.mu_vec;                % Pa·s
else
    mu_vec = logspace(-4,0,64);
end

if isfield(opts,'G_vec') && ~isempty(opts.G_vec)
    G_vec = opts.G_vec;                  % Pa
else
    G_vec = logspace(2,log10(5e5),64);
end

if isfield(opts,'P_inf') && ~isempty(opts.P_inf)
    P_inf = opts.P_inf;                  % Pa
else
    P_inf = 101325;
end

if isfield(opts,'rho') && ~isempty(opts.rho)
    rho = opts.rho;                      % kg/m^3
else
    rho = 1064;
end

% ---------------- Load data via expDataPrep ----------------
% expDataPrep is a script; 'material' will be visible to it from this workspace.
expDataPrep;

% Pull variables to local scope (handle alternate naming)
if exist('tc','var'), tc_each = tc; end
if ~exist('tc_each','var'), error('tc_each not found (expDataPrep must define tc or tc_each).'); end
if ~exist('Rmatrix','var') || ~exist('tmatrix','var') || ~exist('Rdotmatrix','var')
    error('Rmatrix/tmatrix/Rdotmatrix not found. Ensure expDataPrep ran successfully.');
end
if ~exist('Rmax_range','var') && exist('Rmax_each','var')
    Rmax_range = Rmax_each;
end
if ~exist('Rmax_range','var')
    Rmax_range = max(R_d,[],1); % fallback
end
if ~exist('Req_each','var')
    % Estimate Req from dimensional tail if missing (consistent with your prep)
    last10 = max(1, ceil(0.10*size(R_d,1)));
    Req_each = mean(R_d(end-last10+1:end, :), 1);
end

% ---------------- Sort trials by Rmax for clean vertical ordering ----------------
[T, N] = size(Rmatrix);
[Rs, order] = sort(Rmax_range(:));
Rmax_sorted = Rs(:).';
permute_col = @(X) X(:,order);

Rmatrix     = permute_col(Rmatrix);
tmatrix     = permute_col(tmatrix);
Rdotmatrix  = permute_col(Rdotmatrix);
if exist('R_d','var'), R_d = permute_col(R_d); end
Req_each    = Req_each(:,order);
tc_each     = tc_each(:,order);

% Per-trial nondimensional equilibrium radius
Req_nd_each = Req_each ./ Rmax_sorted;          % 1 x N
Req_nd_each = max(Req_nd_each, 1e-12);

% ---------------- Per-trial integrals using Table 1 (non-dimensional S*) ----------------
% S*_v  = -(4/Re_k) * (Rdot*/R*)
% S*_NH = (1/(2*Ca)) * [4*lambda^-1 + lambda^-4 - 5]
% where lambda = R*/Req*.
%
% We accumulate purely kinematic sums so (mu,G) can be inserted analytically:
%   Qv2_sum = Σ_k ∫ (Rdot*/R*)^2 dτ / Rmax_k^2
%   Qh2_sum = Σ_k ∫ [4 λ^-1 + λ^-4 - 5]^2 dτ
%   Qvh_sum = Σ_k ∫ (Rdot*/R*) * [4 λ^-1 + λ^-4 - 5] dτ / Rmax_k
%
% Then, with:
%   A_k(μ) = -4/Re_k = -4*μ/(Rmax_k*sqrt(rho*P_inf))
%   B(G)   =  1/(2*Ca)=  G/(2 P_inf)
%
% we have:
%   S_N^2 = Σ_k ∫ (A_k q_v)^2 dτ = (16 μ^2/(rho*P_inf)) * Qv2_sum
%   S_H^2 = Σ_k ∫ (B q_h)^2 dτ   = (G^2/(4 P_inf^2))     * Qh2_sum
%   S_K^2 = S_N^2 + S_H^2 + 2 Σ_k ∫ (A_k q_v)(B q_h) dτ
%         = S_N^2 + S_H^2 + [ 2*B * Σ_k A_k ∫ q_v q_h dτ ]
%         = S_N^2 + S_H^2 - (4 μ G /(P_inf*sqrt(rho*P_inf))) * Qvh_sum

Qv2_sum = 0.0;
Qh2_sum = 0.0;
Qvh_sum = 0.0;

for k = 1:N
    rstar  = max(Rmatrix(:,k), 1e-12);
    rEq    = max(Req_nd_each(k), 1e-12);
    lambda = rstar ./ rEq;                     % λ = R*/Req*
    q_v    = Rdotmatrix(:,k) ./ rstar;         % (Rdot*/R*)
    q_h    = 4./lambda + lambda.^(-4) - 5;     % NH bracket (no B here)

    tau  = tmatrix(:,k);
    if numel(tau) < 2
        continue; % safeguard
    end
    dtau = [diff(tau); tau(end) - tau(end-1)];

    Qv2_sum = Qv2_sum + sum( (q_v.^2) .* dtau ) / (Rmax_sorted(k)^2);
    Qh2_sum = Qh2_sum + sum( (q_h.^2) .* dtau );
    Qvh_sum = Qvh_sum + sum( (q_v .* q_h) .* dtau ) /  Rmax_sorted(k);
end

% ---------------- Closed-form insertion of μ and G via Re(μ) and Ca(G) ----------------
[MU, GG] = meshgrid(mu_vec, G_vec);

SN2 = (16 .* (MU.^2) ./ (rho .* P_inf)) .* Qv2_sum;
SH2 = ((GG.^2) ./ (4 .* P_inf.^2))      .* Qh2_sum;
SK2 = SN2 + SH2 + (4 .* MU .* GG ./ (P_inf .* sqrt(rho .* P_inf))) .* Qvh_sum;

SN  = sqrt(max(SN2, realmin));
SH  = sqrt(max(SH2, realmin));
SK  = sqrt(max(SK2, realmin));

% ---------------- Uniqueness metric ----------------
% Balanced + magnitude-aware (peaks top-right when SN≈SH and SK is large)
B = 1 - abs(SN - SH) ./ SK;      % balance term in [0,1]
B = max(min(B,1),0);
M = SK ./ max(SK(:));            % magnitude term in [0,1]
Ustar = B .* M;                   % final uniqueness

% If you prefer the geometric-mean uniqueness instead, uncomment:
% U_geo = 1 - sqrt( (SN./SK) .* (SH./SK) );
% U_geo = max(min(U_geo,1),0);

% ---------------- Single manifold (μ–G → U*) ----------------
hfig = figure('Color','w', 'Name','KV uniqueness manifold: mu–G → U_*', ...
              'WindowStyle','normal', 'DockControls','on', ...
              'Visible','on', 'Position',[100 100 880 660]);
movegui(hfig,'center'); drawnow;

s = surf(MU, GG, Ustar, 'EdgeColor','none', 'FaceAlpha',0.96);
set(gca,'XScale','log','YScale','log');
colormap(parula); caxis([0 1]); grid on; box on; view(46,28);
hold on; contour3(MU, GG, Ustar, 10, 'k:');

% ---- Labeling (LaTeX-safe) ----
ax = gca;
set(ax,'TickLabelInterpreter','latex');
xlabel(ax,'$\mu\,[\mathrm{Pa\cdot s}]$','Interpreter','latex');
ylabel(ax,'$G\,[\mathrm{Pa}]$','Interpreter','latex');
zlabel(ax,'$U_\star=(1-|S_N-S_H|/S_K)\,(S_K/\max S_K)$','Interpreter','latex');

% Colorbar with label (older MATLAB pattern)
cbar = colorbar;
ylabel(cbar,'Uniqueness (0..1)','Interpreter','latex');

% Lighting without clashing with your variable named "material"
lighting gouraud; camlight headlight; feval('material','dull');

end