Fun_SetInvFaultBoundary.m

function [isStable, NormSlip] = Fun_SetInvFaultBoundary( ...
    slip, locadep, gridsize, CCritSlip)
% =========================================================================
% Determine whether low-slip margins exist around the main rupture area
% =========================================================================

maxSlip = max(slip(:));

if isempty(maxSlip) || ~isfinite(maxSlip) || maxSlip <= 0
    warning('The inverted slip is zero, negative, or invalid.');
    isStable = true;
    NormSlip = zeros(size(slip));
    return;
end

NormSlip = slip ./ maxSlip;

% Retain only the main rupture region.
NormSlip(NormSlip <= CCritSlip) = 0;

% Number of zero-slip rows from the top edge.
topCumulative = cumsum(sum(NormSlip, 2));
nZeroTop = sum(topCumulative == 0);

% Number of zero-slip rows from the bottom edge.
bottomCumulative = cumsum(flipud(sum(NormSlip, 2)));
nZeroBottom = sum(bottomCumulative == 0);

% Number of zero-slip columns from the left edge.
leftCumulative = cumsum(sum(NormSlip, 1));
nZeroLeft = sum(leftCumulative == 0);

% Number of zero-slip columns from the right edge.
rightCumulative = cumsum(fliplr(sum(NormSlip, 1)));
nZeroRight = sum(rightCumulative == 0);

% If the top row is near or above the surface, it is accepted directly;
% otherwise, require at least one low-slip row at the top edge.
Cond1 = (locadep(1) < gridsize(1)) || (nZeroTop >= 1);

% Require low-slip margins at the lower, left, and right edges.
Cond2 = (nZeroBottom >= 1) && ...
        (nZeroLeft >= 1) && ...
        (nZeroRight >= 1);

isStable = Cond1 && Cond2;

fprintf(['   Zero-slip margins [Top Bottom Left Right] = ','[%d %d %d %d]\n'], ...
    nZeroTop, nZeroBottom, nZeroLeft, nZeroRight);

end