Fun_SetGreenModel4.m

function [g, muSubfault]=Fun_SetGreenModel4(g, locadep, dataVelFolder, grid)
% =========================================================================
% Correct Green's-function amplitudes using depth-dependent shear modulus.
%
% Input:
%   g             : Green's functions [Nt, 3*Nsta, Nsub]
%   locadep       : Subfault depths [Nsub, 1], km
%   dataVelFolder : Folder containing DataEarth.mat
%   grid          : [Ndip, Nstrike]
%
% Output:
%   g             : Shear-modulus-corrected Green's functions
%   muSubfault    : Shear modulus of each subfault, Pa
% =========================================================================
File=fullfile(dataVelFolder, 'DataEarth.mat');
DataEarth=load(File, 'DataEarth');
DataEarth=DataEarth.DataEarth;
% DataEarth columns:
%   1: depth (km)
%   2: Vp (km/s)
%   3: Vs (km/s)
%   4: density (g/cm^3)
muEarth=DataEarth(:,4) .* DataEarth(:,3).^2 .* 1e9;
% Assign each subfault depth to a velocity-model layer
[~, layer]=histc(locadep(:), DataEarth(:,1));
% Handle depths outside the velocity-model range
layer(layer == 0 & locadep(:) < DataEarth(1,1)) = 1;
layer(layer == 0 & locadep(:) >= DataEarth(end,1)) = ...
    size(DataEarth,1);
muSubfault = muEarth(layer);

nsub=prod(grid);
if numel(muSubfault) ~= nsub
    error('The number of shear-modulus values does not match Nsub.');
end
if size(g,3) ~= nsub
    error('The third dimension of g does not match Nsub.');
end
% Reference shear modulus
muRef=3e10;
for iSub = 1:nsub
    g(:,:,iSub)=g(:,:,iSub).*(muSubfault(iSub)./muRef);
end

end