Fun_AutoFaultInversion.m

% =========================================================================
% Automatic finite-fault inversion workflow:
%
%   Stage 1. Preliminary station screening
%   Stage 2. Dense-network station redundancy reduction
%   Stage 3. Iterative finite-fault inversion and fault-plane updating
%   Stage 4. Final fault-plane trimming
%
% Important:
%   This function keeps the same main variables used by the original script:
%
%       ob, g, loca, mm,
%       locasub, locadep, muSubfault,
%       grid, source, index,
%       substf, substfa, syn, obr, res,
%       slip, SeisMoment, Mw
%
% Therefore, the plotting and rupture-output code after this function call
% can remain unchanged
% =========================================================================
function [ob, g, loca, mm, locasub, locadep, muSubfault, grid, source, index, ...
    substf, substfa, syn, obr, res, slip, SeisMoment, Mw]=...
    Fun_AutoFaultInversion(...
    ob, g, loca, mm,locasub, locadep, muSubfault, ...
    PATHg, dataVelFolder, SDR, grid, gridsize, source, ...
    epicenter, srate, IndexEQ, index, depth,fBands, mag, opts)
% User-adjustable inversion settings
if isempty(opts)
    opts.PreScreenIteration = 3;    % Preliminary inversion iterations
    opts.MinStationScreen   = 12;   % Enable initial screening above this number
    opts.DenseStationLimit  = 40;   % Enable redundancy reduction above this number
    opts.MinRemainStation   = 20;   % Try to retain at least this many stations

    opts.MaxOuterIteration  = 5;    % Fault-plane update iteration number
    opts.InnerIteration     = 20;   % Inversion iteration number
    opts.SlipCutoff         = 0.20; % Main rupture threshold: 20% of maximum slip
end
%% ------------------------------------------------------------------------
% Stage 1 and Stage 2: Station selection and consistency checking
Range1=12;
Range2=5;
[ob, g, loca, mm] = Fun_SetInvStations( ...
    ob, g, loca, mm, epicenter, srate, grid, gridsize, source, muSubfault, ...
    Range1,Range2,fBands, mag, opts);

% -------------------------------------------------------------------------
% Stage 3: Iterative inversion and fault-plane updating
[ob, g, loca, ...
 locasub, locadep, muSubfault, grid, source, index, ...
 substf, substfa, syn, obr, res, slip] = Fun_SetInvFault( ...
    ob, g, loca, locasub, locadep, muSubfault, PATHg, dataVelFolder, ...
    SDR, grid, gridsize, source, epicenter, srate, IndexEQ, index, depth, ...
    fBands, mag, opts);

% -------------------------------------------------------------------------
% Stage 4: Final trimming according to the final slip distribution
[g, locasub, locadep, muSubfault, ...
 substf, slip, grid, source] = Fun_SetInvFaultTrim( ...
    g, locasub, locadep, muSubfault, ...
    substf, substfa, grid, source, ...
    gridsize, opts.SlipCutoff);
% Recalculate the final hypocentral linear index
index = (source(2) - 1) * grid(1) + source(1);
% Final seismic moment and moment magnitude
Areas = prod(gridsize) * 1e6;  % km^2 --> m^2
% μ [Pa] * slip [m] * area [m^2] = Moment [N*m]
SeisMoment = sum(slip(:) .* muSubfault(:) .* Areas);
Mw = m2m(SeisMoment);
fprintf('[Result] Final Moment-Magnitude: Mw = %.2f\n', Mw);
fprintf('[Result] Final Grid: [%d, %d]\n', grid(1), grid(2));
fprintf('[Result] Final Source: [%d, %d]\n\n', source(1), source(2));

end