Fun_PlotFaultSlip.m

function Fun_PlotFaultSlip(locaSub, depSub, slip, source, grid, epicenter)
% =========================================================================
% Fun_PlotFaultSlip
%
% Plot final slip distribution on:
%   1) Surface projection of the fault plane (Depth = 0 km)
%   2) Actual subsurface finite-fault plane (Depth = depSub)
%
% Color represents final slip amplitude.
%
% Inputs:
%   locaSub   : N-by-2 subfault-center coordinates:
%               [latitude, longitude]
%
%   depSub    : N-by-1 subfault-center depths (km)
%
%   slip      : Final slip distribution.
%               It can be:
%                 - N-by-1 vector, or
%                 - Ndip-by-Nstrike matrix
%
%   source    : Hypocentral subfault index [idip, istrike]
%
%   grid      : Number of subfaults [Ndip, Nstrike]
%
%   epicenter : Surface epicenter [latitude, longitude]
%
% Grid ordering:
%   index = (istrike - 1) * Ndip + idip
%
% Note:
%   Depth increases downward by:
%       set(gca, 'ZDir', 'reverse')
% =========================================================================

%% ------------------------------------------------------------------------
% Input check
% -------------------------------------------------------------------------
grid = round(grid(:)');
source = round(source(:)');

Ndip = grid(1);
Nstrike = grid(2);
nsub = Ndip * Nstrike;

depSub = depSub(:);

if size(locaSub, 1) ~= nsub
    error(['locaSub contains %d points, but grid requires %d ', ...
           'subfaults (%d x %d).'], ...
           size(locaSub, 1), nsub, Ndip, Nstrike);
end

if numel(depSub) ~= nsub
    error('depSub must contain Ndip * Nstrike elements.');
end

if source(1) < 1 || source(1) > Ndip || ...
   source(2) < 1 || source(2) > Nstrike
    error('source = [idip, istrike] is outside the specified grid.');
end

%% ------------------------------------------------------------------------
% Reshape fault geometry
% -------------------------------------------------------------------------

% Rows: dip direction
% Columns: strike direction
latMat = reshape(locaSub(:, 1), Ndip, Nstrike);
lonMat = reshape(locaSub(:, 2), Ndip, Nstrike);
depMat = reshape(depSub,         Ndip, Nstrike);

% Accept either a vector or an Ndip-by-Nstrike slip matrix.
if isvector(slip)

    if numel(slip) ~= nsub
        error('Vector slip must contain Ndip * Nstrike elements.');
    end

    slipMat = reshape(slip(:), Ndip, Nstrike);

elseif isequal(size(slip), [Ndip, Nstrike])

    slipMat = slip;

else
    error(['slip must be either an N-by-1 vector or an ', ...
           'Ndip-by-Nstrike matrix.']);
end

if any(~isfinite(slipMat(:)))
    error('slip contains NaN or Inf values.');
end

%% ------------------------------------------------------------------------
% Hypocenter coordinates
% -------------------------------------------------------------------------

hypoLat = latMat(source(1), source(2));
hypoLon = lonMat(source(1), source(2));
hypoDep = depMat(source(1), source(2));

%% ------------------------------------------------------------------------
% Color limits: use the same range for surface and subsurface planes
% -------------------------------------------------------------------------

slipMax = max(abs(slipMat(:)));

if slipMax <= eps
    slipMax = 1;
end

% For scalar slip amplitude, the color range starts from zero.
cLim = [0, max(slipMat(:))];

% If slip has negative values, use symmetric color limits.
if min(slipMat(:)) < 0
    cLim = [-slipMax, slipMax];
end

%% ------------------------------------------------------------------------
% Create figure
% -------------------------------------------------------------------------
figure('Color', 'w', 'Position', [160, 100, 800, 600]);
tiledlayout('flow','TileSpacing','compact','Padding','loose');
nexttile;
hold on;
box on;
% ------------------------------------------------------------------------
% Surface projection: z = 0
depProjection = zeros(size(depMat));
hProjection = surf(lonMat, latMat, depProjection, slipMat, ...
    'FaceColor', 'interp', ...
    'FaceAlpha', 0.9, ...
    'EdgeColor', [0.30, 0.30, 0.30], ...
    'EdgeAlpha', 0.75, ...
    'LineStyle', '--', ...
    'LineWidth', 0.7, ...
    'DisplayName', 'Slip projection at surface');

% ------------------------------------------------------------------------
% Actual fault plane below the surface
hFault = surf(lonMat, latMat, depMat, slipMat, ...
    'FaceColor', 'interp', ...
    'FaceAlpha', 0.9, ...
    'EdgeColor', [0.20, 0.20, 0.20], ...
    'EdgeAlpha', 0.75, ...
    'LineWidth', 0.7, ...
    'DisplayName', 'Final slip on fault plane');

shading interp;
lighting gouraud;
% camlight('headlight');
%% ------------------------------------------------------------------------
% Mark surface epicenter
plot3(epicenter(2), epicenter(1), 0, 'kp', ...
    'MarkerFaceColor', 'w', ...
    'MarkerEdgeColor', 'k', ...
    'MarkerSize', 18, ...
    'LineWidth', 1.2, ...
    'DisplayName', 'Surface epicenter');

%% ------------------------------------------------------------------------
% Mark hypocentral subfault
plot3(hypoLon, hypoLat, hypoDep, 'rp', ...
    'MarkerFaceColor', [0.8902, 0.1529, 0.4471], ...
    'MarkerEdgeColor', 'k', ...
    'MarkerSize', 18, ...
    'LineWidth', 1.2, ...
    'DisplayName', 'Hypocentral subfault');
% Connect epicenter and hypocenter.
plot3([epicenter(2), hypoLon], ...
      [epicenter(1), hypoLat], ...
      [0, hypoDep], ...
      'k--', ...
      'LineWidth', 1.0, ...
      'HandleVisibility', 'off');

%% ------------------------------------------------------------------------
% Optional: highlight maximum-slip subfault
% -------------------------------------------------------------------------
% [~, maxSlipIndex] = max(abs(slipMat(:)));
% [idipMax, istrikeMax] = ind2sub([Ndip, Nstrike], maxSlipIndex);
% maxSlipLat = latMat(idipMax, istrikeMax);
% maxSlipLon = lonMat(idipMax, istrikeMax);
% maxSlipDep = depMat(idipMax, istrikeMax);
% plot3(maxSlipLon, maxSlipLat, maxSlipDep, 'o', ...
%     'MarkerSize', 10, ...
%     'MarkerFaceColor', 'none', ...
%     'MarkerEdgeColor', 'k', ...
%     'LineWidth', 1.8, ...
%     'DisplayName', 'Maximum-slip subfault');

%% ------------------------------------------------------------------------
% Axis settings
% -------------------------------------------------------------------------
xlabel('Longitude');
ylabel('Latitude');
zlabel('Depth (km)');
% Geophysical convention: depth increases downward.
set(gca, 'ZDir', 'reverse');
% You used this aspect-ratio convention in the original 2-D plot.
% axis vis3d;
view(45, 25);
% Color represents final slip.
colormap(slanCM('viridis', 60));
clim(cLim);
cb = colorbar;
cb.Label.String = 'Final slip';
cb.FontSize = 14;
set(gca, 'FontSize', 18, 'LineWidth', 1.0, 'FontName', 'Times New Roman', ...
    'Projection', 'perspective');
hold off;

end