Fun_PlotFault.m

function Fun_PlotFault(locaSub, depSub, source, grid, epicenter)
% =========================================================================
% Plot a 3-D filled finite-fault surface and hypocentral subfault.
%
% Input:
%   locaSub   : Subfault-center coordinates [latitude, longitude], N x 2
%   depSub    : Subfault-center depths [km], N x 1
%   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 through:
%       set(gca, 'ZDir', 'reverse')
% =========================================================================
if size(locaSub, 2) ~= 2
    error('locaSub must be an N-by-2 matrix: [latitude, longitude].');
end

if numel(grid) ~= 2 || numel(source) ~= 2
    error('grid and source must both be two-element vectors.');
end

grid = round(grid(:)');
source = round(source(:)');
depSub = depSub(:);

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

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

if numel(depSub) ~= nsub
    error('The length of depSub must equal prod(grid).');
end

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

% Convert [idip, istrike] to linear index.
index = (source(2) - 1) * Ndip + source(1);

%% ------------------------------------------------------------------------
% Reshape subfault-center vectors into a 2-D fault-plane grid
% -------------------------------------------------------------------------

% Each column corresponds to one strike-direction position;
% each row corresponds to one dip-direction position.
latMat = reshape(locaSub(:, 1), Ndip, Nstrike);
lonMat = reshape(locaSub(:, 2), Ndip, Nstrike);
depMat = reshape(depSub,         Ndip, Nstrike);

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

%% ------------------------------------------------------------------------
% Create figure
% -------------------------------------------------------------------------
figure('Color', 'w', 'Position', [160, 100, 800, 600]);
tiledlayout('flow','TileSpacing','compact','Padding','loose');
nexttile;
hold on;
box on;
%% ------------------------------------------------------------------------
% Plot the fault-plane projection at the surface (Depth = 0 km)
% -------------------------------------------------------------------------
deppro = zeros(size(depMat));
hProjection = surf(lonMat, latMat, deppro, depMat, ...
    'FaceColor', 'interp', ...
    'FaceAlpha', 0.7, ...
    'EdgeColor', [0.25, 0.25, 0.25], ...
    'EdgeAlpha', 0.7, ...
    'LineStyle', '--', ...
    'LineWidth', 0.7, ...
    'DisplayName', 'Fault projection at surface');

%% ------------------------------------------------------------------------
% Plot the filled fault surface
% -------------------------------------------------------------------------
% CData = depMat means that the fault color represents depth.
hFault = surf(lonMat, latMat, depMat, depMat, ...
    'FaceColor', 'interp', ...
    'FaceAlpha', 0.7, ...
    'EdgeColor', [0.35, 0.35, 0.35], ...
    'EdgeAlpha', 0.7, ...
    'LineWidth', 0.7, ...
    'DisplayName', 'Finite-fault plane');
% Improve the visual three-dimensional appearance.
shading interp;
lighting gouraud;
% camlight('headlight');

%% ------------------------------------------------------------------------
% Plot surface epicenter
% -------------------------------------------------------------------------
plot3(epicenter(2), epicenter(1), 0, 'kp', ...
    'MarkerFaceColor', [0.8902    0.1529    0.4471], ...
    'MarkerEdgeColor', 'k', ...
    'MarkerSize', 16, ...
    'LineWidth', 1.2, ...
    'DisplayName', 'Surface epicenter');

%% ------------------------------------------------------------------------
% Plot hypocentral subfault
% -------------------------------------------------------------------------
plot3(hypoLon, hypoLat, hypoDep, 'rp', ...
    'MarkerFaceColor', [0.8902    0.1529    0.4471], ...
    'MarkerEdgeColor', 'k', ...
    'MarkerSize', 18, ...
    'LineWidth', 1, ...
    'DisplayName', 'Hypocentral subfault');
% Dashed line between surface epicenter and hypocentral subfault.
plot3([epicenter(2), hypoLon], ...
      [epicenter(1), hypoLat], ...
      [0, hypoDep], 'k--', ...
      'LineWidth', 1, ...
      'HandleVisibility', 'off');


%% ------------------------------------------------------------------------
% Axis and appearance settings
% -------------------------------------------------------------------------
xlabel('Longitude');
ylabel('Latitude');
zlabel('Depth (km)');
% Seismological/geophysical convention: depth grows downward.
set(gca, 'ZDir', 'reverse');
colormap(slanCM('viridis', 60));
% A suitable initial viewing angle
view(45, 25);
axis vis3d;
cb = colorbar;
cb.Label.String = 'Depth';
cb.FontSize = 14;
set(gca, ...
    'FontSize', 16, ...
    'LineWidth', 1, ...
    'Projection', 'perspective','fontname','Times New Roman');
hold off;

end