MainInv.m

%% ========================================================================
%                    >> Rupture inversion process <<
% =========================================================================
% Data:
%   1. TELESEISMIC
%   2. STRONG-MOTION
%   3. GPS
%   4. INSAR
% Refer to:
%   Zhang et al., The 2009 L'Aquila MW 6.3 earthquake: a new technique to
%   locate the hypocentre in the joint inversion of earthquake rupture process
% =========================================================================
% 1. The Green's function library calculations are "universal": a precomputed 
% database of theoretical waveforms for unit moment tensor sources, 
% generated for 1D Earth model (IASP91), covering specific ranges of 
% source depth, epicentral distance, frequency, and sampling rate
% 2. ob_tele.mat and ob_sm.mat are the input files for this study; note 
% that the data format is velocity (which will later be converted to 
% displacement for final inversion)
% 3. The Tele-data contains at least: [ob,loca,da,mm,srate] or [leng,leng_before]
%   where da is [distance (km), azimuth]; mm is station message; ob only
%   contains the Z-component
% 4. The Sm-data is similar to that, but ob contains [Nsample,Component=3,Nsta]
% =========================================================================
% Parameters of fault plane (1)
clc;
EveLoc = [-42.725,173.065];  % [latitude,longitude]
depth  = 19;
Focal  = [219,38,128];   % S/D/R
NGrid  = [8, 25];        % Along dip and strike
LGrid  = [10,10];        % Length of each grid in variable 'NGrid' (km)
SGrid  = [3, 21];        % Location of event in defined grids
% Note that the actual along-strike position of the source is 0, so the
% entire strike-scale is determined jointly by SGrid, LGrid, and NGrid
IsDisp = 1;              % Whether invert to Displacement 1
Files  = 'Results.mat';  % Save file


% Parameters of types of records (2)
% ==================
% >> Teleseismic recording
UseTe = 1;     % Use or not
WayTe = 0;     % Green's function calculation: exact=1, approx=0
DoSTF = 0;     % Whether inverse the apparent source time function
WeightSTF = 0.1;        % Recommended to choose a smaller value, e.g., 0.01
BandTe = [0.01,0.2];    % Filter bandwidth
DataTe = 'Data\ob_tele.mat';  % Loading the initial tele-data (must exists)
DatagTe ='Data\obg_tele.mat'; % Loading/Creating the tele-green-data
PathTe = 'D:\Matlab\MyMechanismFunc\FaultInverse\Green_tele\green_func\';
ShiftTe = 1;            % Whether wave shift after an initial inversion
DoGreenTe = 0;          % Calculate the green function or not (if exists)
% Using DoGreenTe, a new file (obg_tele) will be regenerated subsequently,
% with an additional variable 'g' representing the Green's function
% ==================
% >> Strong-motion recording
UseSm = 1;              % See above for more details
BandSm = [0.02,0.2];
DataSm = 'Data\ob_sm.mat';
DatagSm = 'Data\obg_sm.mat';
PathSm = 'D:\Matlab\MyMechanismFunc\FaultInverse\Green_sm\green_func\';
ShiftSm = 1;
DoGreenSm = 0;
% ==================
% >> GPS recording
UseGPS = 0;
DataGPS = 'Data\ob_gps.mat';
% ==================
% >> INSAR recording
UseINSAR = 0;
Datainsar = 'Data\ob_insar.mat';


% Parameters of the inversion calculation (3)
Tspan= 110;        % Total source time of subfaults
Tmax = 30;         % Maximum rupture time
Vmax = 3;          % Maximum rupture velocity
lambda1 = 50;      % Spatial-smooth factor
lambda2 = 100;     % Temporal-smooth factor
lambda3 = 0.01;    % Minimum scalar moment factor
Rakevary= 1;       % If set to 1, rake is allowed to vary within 45¡ã
% Note that Rakevary must be set to 0 when performing STF inversion
IsEdge  = 1;        % Edge displacement control
%   0 for setting no edge to 0;
%   1 for setting all edges except for the upper edge to 0;
%   2 for setting all edges to 0

fprintf('# [1] Parameters are ready\n\n');

%% ========================================================================
%                       >> Calculation Prepare <<
% =========================================================================
% The Green's function corresponding to a specific source fault will be 
% calculated and expressed as:
% 1. The teleseismic Green function:   gt=gt(:,i,k,j)
%   The j-th teleseismic record, generated by the k-th sub-fault, 
%   corresponding to the time series of the Green's function for the i-th 
%   unit moment tensor component
% 2. The Strong-Motion Green function: gsm=gsm(:,i,k,j)
%   g(:,i,k,3*(j-1)+1) means E-component of the j-th station
%   g(:,i,k,3*(j-1)+2) means N-component of the j-th station
%   g(:,i,k,3*(j-1)+3) means U-component of the j-th station
% The velocity waveform will be integrated and filtered in related function
% then the data length will be standardized and merged
% =========================================================================
% Calculate Shear Mudulus at each subfaults
load('earth.mat');
[miu1,miun]=getmiu(depth,Focal,NGrid,LGrid,SGrid,EveLoc,earth);

%                   >> 1. Teleseismic processing <<
% Notice:
% During this process, the data will be saved as obg_tele or obg_sm. When
% needed, you can use the precomputed values instead of recalculating each
% time, referring to the logical values of DoGreenTe and DoGreenSm
% (Just start by doing the calculation directly)
% Indeed:
% The saved data has not been filtered, so filtering will be 
% performed again during function <getobg_tele*> and <getobg_sm>
if DoSTF==1 % Approx solution & Cal STF
    [obt,gt,locat,mmt,s1,leng,t,t0]=Fun_getobg_tele4(UseTe,DoGreenTe,DataTe,DatagTe,...
        EveLoc,depth,Focal,NGrid,LGrid,SGrid,PathTe,miun,BandTe,IsDisp);
    tdt=t-repmat(t0',size(t,1),1);
    if WayTe==1 % Exact solution & Cal STF
        [obt,gt,locat,mmt,s1,leng]=Fun_getobg_tele2(UseTe,DoGreenTe,DataTe,DatagTe,...
            EveLoc,depth,Focal,NGrid,LGrid,SGrid,PathTe,miun,BandTe,IsDisp);
    end
elseif WayTe==0 % Approx solution & No STF (**Recommendation**)
    [obt,gt,locat,mmt,s1,leng]=Fun_getobg_tele3(UseTe,DoGreenTe,DataTe,DatagTe,...
        EveLoc,depth,Focal,NGrid,LGrid,SGrid,PathTe,miun,BandTe,IsDisp);
elseif WayTe==1 % Exact solution & No STF
    [obt,gt,locat,mmt,s1,leng]=Fun_getobg_tele2(UseTe,DoGreenTe,DataTe,DatagTe,...
        EveLoc,depth,Focal,NGrid,LGrid,SGrid,PathTe,miun,BandTe,IsDisp);
end

%                  >> 2. Strong-Motion processing <<
[obsm,gsm,locasm,mmsm,s2]=Fun_getobg_sm(UseSm,DoGreenSm,DataSm,DatagSm,...
    EveLoc,depth,Focal,NGrid,LGrid,SGrid,PathSm,miun,BandSm,IsDisp);

% Easy processing and zero padding
[obt,gt]=length_unification(obt,gt,obsm);
% Verify uniform sampling rate
srate=set_srate(s1,s2);
% Combine structure
ob=[obt,obsm];
g=cat(4,gt,gsm);
loca=[locat;locasm];
mm=[mmt;mmsm];
clear gt gsm; % clear Green function
% Weighted average
[ob,g,mob]=normalize_obg(ob,g,obt,obsm,UseTe,UseSm);
% Construct finite fault matrix
[G,sGDT0,hdata,meangreen,sumob,trup_con,locasub,locadep,trup]=get_Gmatrix(...
    ob,g,loca,EveLoc,depth,Focal,NGrid,LGrid,SGrid,Tspan,Tmax,Vmax,srate,...
    lambda1,lambda2,lambda3,Rakevary,IsEdge);
if DoSTF==1
    obtele=ob(:,1:size(locat,1));
    gtele=g(:,:,:,1:size(locat,1));
    [astftele,synastf]=getastftele(obtele,gtele,SGrid,NGrid,Focal,Tspan,BandTe,s1);
    Gastf=get_stfinvG(NGrid,trup_con,tdt,Tspan*srate);
    astftele=astftele(:)/srate;
    mobstf=max(max(Gastf))/WeightSTF/max(max(G))*2;
    astftele=astftele/mobstf;
    Gastf=Gastf/mobstf;
    grida=NGrid;
    save('Data\ResultSTF.mat','gtele','SGrid','grida','Focal','Tspan','BandTe','s1','Rakevary','mobstf');
else
    Gastf=[];
    astftele=[];
end

%                       >> 3. GPS processing <<
[slipg,Ggps,weit_gps,locagps] = getobg_gps2(UseGPS,DataGPS,sumob,...
    locasub,locadep,Focal,NGrid,LGrid,trup_con,Rakevary);

%                       >> 4. INSAR processing <<
[slipsar,Gsar,weit_sar,locasar,nsar,dr] = getobg_insar2(UseINSAR,Datainsar,...
    sumob,locasub,locadep,Focal,NGrid,LGrid,trup_con,Rakevary);

% Construct the Green's function matrix
G_jt = [G; Ggps*weit_gps; Gsar*weit_sar];
sGDT=[Gastf;sGDT0];
% Construct the observation record matrix
ObsUsing = [ob(:)/meangreen;slipg(:)*weit_gps;slipsar*weit_sar;astftele(:);zeros(...
    size(sGDT0,1),1)];
fprintf('# [2] Data-Input are ready\n\n');
clear g G Ggps Gsar;



%% ========================================================================
%                       >> Calculation Start <<
% =========================================================================
Ccrit=10^(-5); % Threshold
MaxIter=80; % 200 is enough
st1=datetime('now'); fprintf('# [Now] %s \n',st1);
[Xsolve,rsq]=Fun_MainInverse(G_jt,sGDT,ObsUsing,Ccrit,MaxIter);
% req: sum of squared residuals (SSE)
st2=datetime('now');
fprintf('# [Process] %s time \n',st2-st1);
% =========================================================================
%               >> Recalculate again when shifting finished <<
% =========================================================================
st1=datetime('now');
syn=reshape(G_jt(1:numel(ob),:)*(Xsolve*meangreen),size(ob)); % Synthesis Record
[Xsolve,ob0,syn,astftele2,synastf2]=Fun_MainInverseShift(ShiftTe*UseTe,ShiftSm*...
    UseSm,Xsolve,ob,syn,obt,obsm,ObsUsing,meangreen,G_jt,sGDT,srate,DoSTF,Ccrit,MaxIter);
st2=datetime('now');
fprintf('# [Re-Process] %s time \n',st2-st1);
fprintf('# [Final] %s \n',st2);
% =========================================================================
%                           >> Do More Info <<
% =========================================================================
if DoSTF==1
    synstf=sGDT(1:size(obt,2)*Tspan*srate,:)*Xsolve*mobstf*srate;
end
misfta0=sum(sum((ob0-syn).^2))/sum(sum(ob0.^2));
[X,ob,syn] = get_obgX(Xsolve,ob0,syn,mob,Tspan,trup_con,prod(NGrid),srate,Rakevary);
obsm = ob(:,size(obt,2)+1:end);
synsm = syn(:,size(obt,2)+1:end);
obsm = reshape(obsm,[size(ob,1),size(obsm,2)/3,3]);
obsm = permute(obsm,[1,3,2]);
synsm = reshape(synsm,[size(ob,1),size(obsm,3),3]);
synsm = permute(synsm,[1,3,2]);
syng = G_jt(numel(ob)+1:numel(ob)+numel(slipg),:)/weit_gps*Xsolve;
synsar = G_jt(numel(ob)+numel(slipg)+1:end,:)/weit_sar*Xsolve;
% ==============
% slip: total slip amount for each subfault  
% stf: overall source time function  
% slipx: slip component in one local coordinate direction of a fault  
% slipy: slip component in another local coordinate direction of a fault  
% slipe: eastward slip component  
% slipn: northward slip component  
% zslip: vertical slip component  
% stfsub_out: slip time function for each subfault
% ==============
% [X,ob,syn] = get_obgX(Xsolve,ob,syn,mob,lensub,trup_con,prod(grid),srate);
% [slip,stf,slipx,slipy,slipe,slipn,zslip] = get_result(X,lensub,grid,...
%    gridsize,fault,miun,srate);
[slip,stf,slipx,slipy,slipe,slipn,zslip,stfsub_out] = get_result2(X,Tspan,NGrid,...
    LGrid,Focal,miun,srate);
stfsub_out=stfsub_out*srate;
% M=¦ÌDS
Mw=m2m(sum(slip(:).*miun(:)*prod(LGrid)*1e6));
fprintf('\n# Moment magnitude= %.2f\n',Mw);
msftj=sum(sum((ob-syn).^2))/sum(sum((ob).^2));
if UseTe
    obtelef=ob(:,1:size(obt,2));
    syntele=syn(:,1:size(obt,2));
    msftt=sum(sum((obtelef-syntele).^2))/sum(sum((obtelef).^2));
    msftti=sum((obtelef-syntele).^2)./sum((obtelef).^2);
end
if UseSm
    msftsm=sum(sum(sum((obsm-synsm).^2)))/sum(sum(sum((obsm).^2)));
    msftsmi(:,:)=sum((obsm-synsm).^2)./sum((obsm).^2);
end

fprintf('# [3] Finished\n\n');
% clear G_jt sGDT sGDT0;
% save(Files); % Save results (need to exclude G_jt and other large)