Fun_SetGreenModel2.m

function [ob,loca,g,mm,muSubfault]=Fun_SetGreenModel2(dataVelFolder,locadep,loca,ob,g,mm,epicenter,grid)
%% ------------------------------------------------------------------------
% Velocity model: [depth Vp Vs rho]
DataEarth=load([dataVelFolder,'\DataEarth.mat']);
DataEarth=DataEarth.DataEarth;
% Calculate the shear modulus (=rho*Vs^2) for each velocity-model layer
muEarth=DataEarth(:,4) .* DataEarth(:,3).^2 .* 1e9;
% Assign a shear modulus to each subfault
[~,Layer]=histc(locadep(:), DataEarth(:,1));
% Check all subfault
if any(Layer==0)
    warning(['Some subfault depths are outside the range of the velocity ', ...
             'model. The nearest valid velocity-model layer will be used.']);
    Layer(Layer == 0 & locadep(:) < DataEarth(1,1)) = 1;
    Layer(Layer == 0 & locadep(:) >= DataEarth(end,1)) = size(DataEarth,1);
end
muSubfault=muEarth(Layer);

% Correct the Green's-function Amplitude for each subfault
muRef=3e10;      % Reference shear modulus (Pa)
nsub=prod(grid);
for iSub=1:nsub
    % The correction is applied to all stations and components;
    % This makes the smoothing constraint act on slip rather than seismic
    % moment when the shear modulus varies with depth
    g(:,:,iSub)=g(:,:,iSub).*(muSubfault(iSub)/muRef);
end


%% ------------------------------------------------------------------------
% Screen stations according to distance or others
% You can change the parameters
MaxDist=400;
MinDist=5;
MinNsta=12;
MinAzi=0.01;
len1=size(loca,1);
% Do screen
[ob,loca,g,mm]=dist_azim_del(ob,loca,g,mm,epicenter,MaxDist,MinNsta,MinDist,MinAzi);
len2=size(loca,1);

%% ------------------------------------------------------------------------
% Update the number of stations after station screening
nsta = size(loca, 1);
if size(ob, 2) ~= 3 * nsta
    error(['The number of waveform channels is inconsistent with the ', ...
           'number of retained stations.']);
end
if size(g, 2) ~= 3 * nsta
    error(['The second dimension of g is inconsistent with the number ', ...
           'of retained stations.']);
end
% Select a common waveform length
% Normalize each waveform by its own maximum absolute amplitude
oba = abs(ob);
obMax = max(oba, [], 1);
obMax(obMax==0) = 1; % Avoid division by zero
oba = oba./obMax;

% Sum normalized amplitudes over all channels
cumoba = sum(oba, 2);
% Accumulate the amplitude over time
cumoba = cumsum(cumoba);
% Normalize the cumulative curve to the range [0,1]
if cumoba(end)==0
    error('All observation waveforms are zero after station screening');
end
cumoba = cumoba./cumoba(end);
% Retain the waveform until 90 percent of the cumulative normalized
% amplitude has been reached
CCrit = 0.9;
[~, Iscut]=min(abs(cumoba - CCrit));
ob = ob(1:Iscut,:);
g  = g(1:Iscut,:,:);

fprintf('\n[Setting] Manually removed %d stations\n',len2-len1);
fprintf('    90%%-Interval reduces the wavelength to %d/%d \n\n',Iscut,length(cumoba));

end