Fun_SetFaultModel.m

function [fModel, fGrid] = Fun_SetFaultModel(mag, Type)
% =========================================================================
% Estimate rupture dimensions, subfault grid, source position, and
% inversion frequency band.
%
% This function follows the original logic based on Wells and Coppersmith
% (1994).
%
% Input:
%   mag       : Moment magnitude
%   Type      : Fault type
%               1 = reverse / thrust fault
%               2 = strike-slip fault
%               3 = normal fault
%
% Output:
%   faultModel.Ruplen : Rupture length [km]
%   faultModel.Rupwid  : Rupture width  [km]
%   faultModel.fBands : Frequency band [Hz]
%
%   sourceGrid.gridsize : Subfault size [dip-direction, strike-direction]
%   sourceGrid.grid     : Number of subfaults [dip-direction, strike-direction]
%   sourceGrid.source   : Initial source position [dip-direction, strike-direction]
% =========================================================================
if nargin < 2
    error('Both are required');
end
if ~isscalar(mag) || ~isfinite(mag)
    error('mag must be a finite scalar');
end
if ~isscalar(Type) || ~ismember(Type, [1, 2, 3])
    error('faultType must be 1, 2, or 3');
end

%% ------------------------------------------------------------------------
% Estimate rupture length and width
% -------------------------------------------------------------------------
switch Type
    case 1
        % Reverse / thrust fault
        Ruplen = 10^(-2.42 + 0.58 * mag);
        Rupwid  = 10^(-1.61 + 0.41 * mag);
    case 2
        % Strike-slip fault
        Ruplen = 10^(-2.57 + 0.62 * mag);
        Rupwid  = 10^(-0.76 + 0.27 * mag);
    case 3
        % Normal fault
        Ruplen = 10^(-1.88 + 0.50 * mag);
        Rupwid  = 10^(-1.14 + 0.35 * mag);
end

%% ------------------------------------------------------------------------
% Select subfault size
%
% gridsize = [dip-direction, strike-direction]
% grid     = [number along dip, number along strike]
% source   = initial source location in the subfault grid
% -------------------------------------------------------------------------
if mag < 6.5
    gridsize = [2, 2];
elseif mag >= 6.5 && mag < 7.5
    gridsize = [5, 5];
elseif mag >= 7.5 && mag < 8.5
    gridsize = [10, 10];
else
    gridsize = [20, 20];
end

% Keep the original grid construction logic:
% Make the number of subfaults odd and put the initial source near center.
grid = [ ...
    ceil(ceil(Rupwid  / gridsize(1)) / 2) * 2 + 1, ...
    ceil(ceil(Ruplen / gridsize(2)) / 2) * 2 + 1];
source = [ ...
    ceil(ceil(Rupwid  / gridsize(1)) / 2) + 1, ...
    ceil(ceil(Ruplen / gridsize(2)) / 2) + 1];

%% ------------------------------------------------------------------------
% Select frequency band
% -------------------------------------------------------------------------
if mag >= 8.5
    fBands = [0.02, 0.05];
elseif mag >= 7.5 && mag < 8.5
    fBands = [0.02, 0.10];
elseif mag >= 6.5 && mag < 7.5
    fBands = [0.02, 0.20];
else
    fBands = [0.02, 0.50];
end

%% ------------------------------------------------------------------------
% Store outputs
% -------------------------------------------------------------------------
fModel.Ruplen = Ruplen;
fModel.Rupwid = Rupwid;
fModel.fBands = fBands;

fGrid.gridsize = gridsize;
fGrid.grid = grid;
fGrid.source = source;

fprintf('\n[Setting] Fault model parameters:\n');
fprintf('   Rupture length : %.2f km\n', Ruplen);
fprintf('   Rupture width  : %.2f km\n', Rupwid);
fprintf('   Subfault size  : [%.1f, %.1f] km\n', gridsize(1), gridsize(2));
fprintf('   Grid size      : [%d, %d]\n', grid(1), grid(2));
fprintf('   Initial source : [%d, %d]\n',source(1), source(2));
fprintf('   Frequency band : %.2f-%.2f Hz\n',fBands(1), fBands(2));
end