% FFT method from slip to stress
% Ref: Fast computation of static stress changes on 2D faults from final slip distribution
% J. Ripperger and P.M. Mai
% GRL 2004

clc
clear all; close all;

%
% load source
% x(km) y(km) strike-slip(cm) dip-slip(cm)
% left-slip is positive and reverse-faulting is positive.
%

load source_stress.dat

source=source_stress;

% option to mirror the source
% this only make sense for vertical strike slip faulting
% the image source is opposite in polarity, symmitric in locations
mirror = 1;

slim = 20;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1 step
% generate arrays of slip  
% slip is a function of x and y 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% space domain sampling interval
% has to be smaller than 1/2*(highest frequency in signal)

%dx=0.2;
%dy=0.2;

%x=0:dx:30;
%y=0:dy:20;

dx=1;
dy=1;

x=0:dx:670;
y=0:dy:260;

ss=source(:,3)/100;
ds=source(:,4)/100;

ss=reshape(ss,length(x),length(y));
ss=flipud(ss');

ds=reshape(ds,length(x),length(y));
ds=flipud(ds');

figure(1)
subplot(3,1,1)
imagesc(x,y,ss), hold on; 
xlabel('along strike distance in km');
ylabel('down dip distance in km');
title('left-slip in meters');
colorbar;
caxis([-10 0])

figure(2)
subplot(3,1,1)
imagesc(x,y,ds), hold on;  
xlabel('along strike distance in km');
ylabel('down dip distance in km');
title('dip-slip in meters');
colorbar;
caxis([0 20])

%return;

if (mirror ==1) 
disp('mirror the source to compute stress drop in half space')
disp('dangrous for dip slip')
y0 = y;
y = [fliplr(-y) y];

ss0 = ss;
ss = [0*ss;ss];
%figure,imagesc(x,y,ss),colorbar

ds0 = ds;
ds = [0*ds;ds];
%figure,imagesc(x,y,ds),colorbar

end

% return;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2 step
% fourier transform of slip function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


NFFTX = 2^(nextpow2(length(x)))
NFFTY = 2^(nextpow2(length(y)))

% truncation length in km has to be large enough to avoid rippling in
% frequency domain

fx=dx*NFFTX    
fy=dy*NFFTY

kx = (-NFFTX/2:NFFTX/2-1)/fx; % wavenumber in km^-1
ky = (-NFFTY/2:NFFTY/2-1)/fy;  

FDs = fftshift(dx*dy*fft2(ss,NFFTY,NFFTX));
FDd = fftshift(dx*dy*fft2(ds,NFFTY,NFFTX));

% figure
% imagesc(kx,ky,(abs((FDs))))
% xlabel('kx in km^-1')
% ylabel('ky in km^-1')
% title('fourier transform of strike slip')
% colorbar;

% figure
% imagesc(kx,ky,(abs((FDd))))
% xlabel('kx in km^-1')
% ylabel('ky in km^-1')
% title('fourier transform of dip slip')
% colorbar;

% check by inverse transform the slip
% tmp=ifft2(fftshift(FDs/dx/dy));
% S=tmp(1:length(y),1:length(x));
% figure
% imagesc(x,y,abs(S))
% xlabel('along strike distance in km');
% ylabel('down dip distance in km');
% title('inverse fourier transform of strike slip')
% colorbar;

% tmp=ifft2(fftshift(FDd/dx/dy));
% S=tmp(1:length(y),1:length(x));
% figure
% imagesc(x,y,abs(S))
% xlabel('along strike distance in km');
% ylabel('down dip distance in km');
% title('inverse fourier transform of dip slip')
% colorbar;

%return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3 step 
% compute the transfer function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[KX,KY]=meshgrid(kx,ky);

% elastic parameters 
poission = 0.25
mu = 40000 % Mpa
lambda = mu*(2*poission/(1-2*poission)) % Mpa

% static stiffness function 

% slip-parallel stress change
% CAUTION
% different for strike-slip and dip-slip as KX and KY is not equalvalent in the following formula
transferXs = -pi*mu*(2*(lambda+mu)/(lambda+2*mu)*KX.^2+KY.^2)./sqrt(KX.^2+KY.^2);
transferXd = -pi*mu*(2*(lambda+mu)/(lambda+2*mu)*KY.^2+KX.^2)./sqrt(KX.^2+KY.^2);

% zero the singularity
transferXs(isnan(transferXs))=0;
transferXd(isnan(transferXd))=0;

% slip-perpendicular stress change
% CAUTION
% same for strike-slip and dip-slip
transferY = -pi*mu*(2*(lambda+mu)/(lambda+2*mu)-1).*KX.*KY./sqrt(KX.^2+KY.^2);

transferY(isnan(transferY))=0;

% figure
% imagesc(kx,ky,transferX)
% title('stiffness function in X')
% colorbar

% figure
% imagesc(kx,ky,transferY)
% title('stiffness function in Y')
% colorbar

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4 step 
% multiplication in frequency fomain = convoluation in space domain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


StressXs=transferXs.*FDs/1000; % Mpa
StressXd=transferXd.*FDd/1000; % Mpa

% figure
% imagesc(kx,ky,(abs(StressXs)))
% title('stiffness function: KX*D')
% colorbar

StressYs=transferY.*FDs/1000;
StressYd=transferY.*FDd/1000;

% figure
% imagesc(kx,ky,log(abs(StressYs)))
% title('stiffness function: KY*D')
% colorbar

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5 step 
% inverse fourier transform
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% get the sign right

signs = - sign(mean(mean(ss)));
signd = - sign(mean(mean(ds)));

tmp = ifft2(fftshift(StressXs/dx/dy));
stressxs = signs * tmp(1:length(y),1:length(x));

tmp = ifft2(fftshift(StressYd/dx/dy));
stressyd = signd * tmp(1:length(y),1:length(x));

tmp = ifft2(fftshift(StressYs/dx/dy));
stressys = signs * tmp(1:length(y),1:length(x));

tmp = ifft2(fftshift(StressXd/dx/dy));
stressxd = signd * tmp(1:length(y),1:length(x));


if (mirror == 1)
stressxs = stressxs+flipud(-stressxs);
stressxs = stressxs(1+length(y0):length(y),:);

stressxd = stressxd+flipud(-stressxd);
stressxd = stressxd(1+length(y0):length(y),:);

stressys = stressys+flipud(-stressys);
stressys = stressys(1+length(y0):length(y),:);

stressyd = stressyd+flipud(-stressyd);
stressyd = stressyd(1+length(y0):length(y),:);

y = y0;
ss = ss0;
ds = ds0;
end

figure(1)
subplot(3,1,2)
imagesc(x,y,real(stressxs))
title('stress drop along strike from strike-slip in Mpa');
colorbar;
caxis([-slim slim])

figure(2)
subplot(3,1,2)
imagesc(x,y,real(stressyd))
title('stress drop along strike from dip-slip in Mpa');
colorbar;
caxis([-slim slim])


figure(1)
subplot(3,1,3)
imagesc(x,y,real(stressys))
title('stress drop along dip from strike-slip in Mpa');
colorbar;
caxis([-slim slim])
%print -depsc2 'stress_drop1'

figure(2)
subplot(3,1,3)
imagesc(x,y,real(stressxd))
title('stress drop along dip from dip-slip in Mpa');
colorbar;
caxis([-slim slim])

%print -depsc2 'stress_drop2'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% step 6 
% combine the dip-slip component and strike-slip component together
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(3)
ts = sqrt(ss.^2+ds.^2);

subplot(3,1,1)
imagesc(x,y,ts), hold on;
xlabel('along strike distance in km');
ylabel('down dip distance in km');
title('total slip magnitude in meters');
colorbar;
caxis([0 20])

tss = real(stressxs+stressyd);
subplot(3,1,2)
imagesc(x,y,tss), hold on;
title('stress drop along strike');
colorbar;
caxis([-20 20])

tsd = real(stressxd+stressys);
subplot(3,1,3)
imagesc(x,y,tsd), hold on;
title('stress drop along dip');
colorbar;
caxis([-20 20])

%subplot(2,1,2)
%imagesc(x,y,sqrt(tsd.^2+tss.^2)), hold on;
%title('stress drop along dip');
%colorbar;
%caxis([-20 20])

% print -depsc2 'stress_drop3'

% save it 
[xx,yy] = meshgrid(x,y);
N=length(x)*length(y);
tmp = [reshape(xx,N,1) reshape(yy,N,1) reshape(tss,N,1) reshape(tsd,N,1)];
save stressdrop.dat tmp -ascii

