Fault-Inversion-Auto/Fun_SetSourceDep.m at main · Yuxiao-git2025/Fault-Inversion-Auto · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
% =========================================================================
% Adjust the designed finite-fault grid to the real earthquake depth
function [source, index, locaSub, depSub, depthInfo] = ...
    Fun_SetSourceDep(SDR, grid, gridsize, source, loca, epicenter, realDepth)
if numel(SDR) < 2
    error('SDR must contain at least strike and dip.');
end
if numel(grid) ~= 2
    error('grid must be a two-element vector.');
end
if numel(gridsize) ~= 2
    error('gridsize must be a two-element vector.');
end
if numel(source) ~= 2
    error('source must be [idip, istrike].');
end
if numel(epicenter) ~= 2
    error('epicenter must be [latitude, longitude].');
end
if ~isscalar(realDepth) || ~isfinite(realDepth)
    error('realDepth must be a finite scalar.');
end
grid = round(grid(:)');
gridsize = double(gridsize(:)');
source = round(source(:)');
% Main-script convention:
%   grid   = [Ndip, Nstrike]
%   source = [idip, istrike]
%
% get_subloca convention:
%   grid   = [Nstrike, Ndip]
%   source = [segmentNumber, istrike, idip]
gridsub = grid([2, 1]);
sourcesub = [1, source([2, 1])];

%% ------------------------------------------------------------------------
% Calculate the initial subfault locations and depths
% -------------------------------------------------------------------------
[locaSub, depInitial] = get_subloca( ...
    SDR, ...
    gridsub, ...
    gridsize, ...
    sourcesub, ...
    loca, ...
    epicenter);
depInitial = depInitial(:);
numberOfSubfaults = prod(grid);

if numel(depInitial) ~= numberOfSubfaults
    error(['The number of calculated subfault depths (%d) does not ', ...
           'match prod(grid) (%d).'], ...
           numel(depInitial), numberOfSubfaults);
end

%% ------------------------------------------------------------------------
% Find the designed hypocentral subfault
% -------------------------------------------------------------------------
indexInitial=(source(2) - 1) * grid(1) + source(1);
if indexInitial < 1 || indexInitial > numberOfSubfaults
    error('The initial source index is outside the subfault grid.');
end

%% ------------------------------------------------------------------------
% Shift the whole fault so that the hypocentral subfault reaches real depth
% -------------------------------------------------------------------------
depthShift = realDepth - depInitial(indexInitial);
depSub = depInitial + depthShift;

%% ------------------------------------------------------------------------
% If the upper edge is above the surface, move the source downward
% -------------------------------------------------------------------------
Counts = sum(depSub < 0);
if Counts > 0
    % Each dip-direction row contains grid(2) subfaults
    % Use ceil so that the source remains on an integer grid index and
    % moves sufficiently far downward
    Nshift = ceil(Counts / grid(2));
    source(1) = source(1) + Nshift;
    % Keep the source inside the dip-direction grid.
    source(1) = min(max(source(1), 1), grid(1));

    % Recalculate the grid location and depths after changing source.
    sourcesub = [1, source([2, 1])];
    [locaSub, depInitial] = get_subloca( ...
        SDR, ...
        gridsub, ...
        gridsize, ...
        sourcesub, ...
        loca, ...
        epicenter);

    depInitial = depInitial(:);

    indexInitial = ...
        (source(2) - 1) * grid(1) + source(1);

    depthShift = realDepth - depInitial(indexInitial);
    depSub = depInitial + depthShift;
end

%% ------------------------------------------------------------------------
% Final source index
% -------------------------------------------------------------------------
index = (source(2) - 1) * grid(1) + source(1);
if depSub(index) < 0
    warning(['The corrected hypocentral subfault is still above the ', ...
             'surface. Please check grid, source, and fault geometry.']);
end

%% ------------------------------------------------------------------------
% Store diagnostic information
% -------------------------------------------------------------------------
depthInfo.depInitial = depInitial;
depthInfo.depFinal = depSub;
depthInfo.depthShift = depthShift;
depthInfo.indexInitial = indexInitial;
depthInfo.indexFinal = index;
depthInfo.numberOfNegativeDepth = sum(depSub < 0);
depthInfo.minimumDepth = min(depSub);
depthInfo.maximumDepth = max(depSub);
depthInfo.hypocenterDepth = depSub(index);


fprintf('\n[Setting] Subfault-Depth correction:\n');
fprintf('   Source index    : %d\n', depthInfo.indexFinal);
fprintf('   Shift           : %.2f km\n', depthInfo.depthShift);
fprintf('   Hypocenter      : %.2f km\n', depthInfo.hypocenterDepth);
fprintf('   Minimum subfault: %.2f km\n', depthInfo.minimumDepth);
fprintf('   Maximum subfault: %.2f km\n', depthInfo.maximumDepth);
end