% MATLABScript for obtaning the analytical solution for 2D wave propagation
% in a viscoelastic medium, based on the Appendix B of paper of Carcione et 
% al. (1988)(1), corrected in (2).
%
% The source used is a Ricker wavelet in the y velocity component.
%
% A figure of the source temporal function and its inverse Fourier
% transform are provided to verify the good working of the method. Also
% works for elastic variables by substituting the values of vp and vs by 
% real constant ones (vp_0 & vs_0).
%
% Attenuation is implemented by GMB-EK rheology presented in (4) and
% thoroughfully explained in (3).
%
% In order to make a proper comparison with numerical data it is
% recommended to either use the same rheology in both the numerical and the
% analytical. It is also possible to get an almost constant Q value by 
% increasing a lot the number of mechanisms used in both the numerical and
% the analytical solutions, thus minimizing the effects of the Q fitting in
% the comparison. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  By Josep de la Puente, LMU Geophysics, 2005  %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

home
clear all;
close all;

disp('======================================')
disp('COMPUTATION OF VISCOELASTIC ANALYTICAL')
disp('             SOLUTION                 ')
disp('======================================')
disp('   ((c) Josep de la Puente Alvarez  ')
disp('        LMU Geophysics 2005)')
disp('======================================')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Source parameters (Ricker wavelet):
f0=10;               % cutoff frequency
t0=0.1;              % time shift
nu=0.5;              % weight of exponential (for Carcione's source)
epsilon=1.0;         % weight of cosinus (for Carcione's source)

%Frequency domain:
w_max=1001;          % Maximum frequency (minimum freq = -w_max). Is better to keep an odd number 
w_inc=0.1;           % Increment in frequency

%Position of receiver:
x=0.5;
y=0.5;

%Force of the source:
F=1.0;               % Constant magnitude. Only alters amplitude of the seismograms

%Material parameters:
rho=1.0;             % Density
vp_0=sqrt(3);        % P-wave velocity (elastic)
vs_0=1;              % S_wave velocity (elastic)

%Attenuation Q-Solver parameters
n=3;                 % Number of mechanisms we want to have
QPval=20;            % Desired QP constant value
QSval=10;            % Desired QS constant value
freq=f0;             % Central frequency of the absorption band (in Hertz). Good to center it
                     %  at the source's central frequency
f_ratio=100;         % The ratio between the maximum and minimum frequencies of our bandwidth
                     % (Usually between 10^2 and 10^4)
                        
%Outputting variables
toutmax=1.0;         % Maximum value of t desired as output


disp('--------------------------------------')
disp(' GEOMETRY AND MATERIAL SETTINGS:');
disp('--------------------------------------')
disp(strcat('Density=  ',num2str(rho)));
disp(strcat('VP= ',num2str(vp_0)));
disp(strcat('VS= ',num2str(vs_0)));
disp('Source at coordinates: (0,0)');
disp(strcat('Receiver at coordinates: (',num2str(x),',',num2str(y),')'));
disp(' ');

disp('--------------------------------------')
disp(' ATTENUATION SETTINGS:');
disp('--------------------------------------')
disp(strcat('Number of mechanisms=  ',num2str(n)));
disp(strcat('Q for P-wave=  ',num2str(QPval)));
disp(strcat('Q for S-wave=  ',num2str(QSval)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% PROBLEM SETUP %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Derived quantities initialization
r=sqrt(x^2+y^2);     % Distance to the receiver
w0=f0/(2*pi);        % Conversion to angular velocity
w=[0:w_inc:w_max];
w=2*w-w_max;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% GMB-EK COMPLEX M SOLUTION %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Equation system initialization %%
kmax=2*n-1;             %Number of equations to solve (the system is overdetermined)
AP=zeros(kmax,n);       %Initialization of system matrix (P wave)
AS=zeros(kmax,n);       %Initialization of system matrix (S wave)
QP=ones(kmax,1)/QPval;  %Desired values of Q for each mechanism inverted (P wave)
QS=ones(kmax,1)/QSval;  % " " (S wave)
YP=zeros(n,1);
YS=zeros(n,1);

%% Selection of the logarithmically equispaced frequencies
wmean=2*pi*freq;
wmin_disc=wmean/sqrt(f_ratio); 

for j=1:kmax
    w_disc(j)=exp(log(wmin_disc)+(j-1)/(kmax-1)*log(f_ratio));
end

%% Filling of the linear system matrix %%
for m=1:kmax
    for j=1:n
        AP(m,j)=(w_disc(2*j-1).*w_disc(m)+w_disc(2*j-1).^2/QPval)./(w_disc(2*j-1).^2+w_disc(m).^2);
    end
end

for m=1:kmax
    for j=1:n
        AS(m,j)=(w_disc(2*j-1).*w_disc(m)+w_disc(2*j-1).^2/QSval)./(w_disc(2*j-1).^2+w_disc(m).^2);
    end
end

%% Solving of the system %%
YP=AP\QP;
YS=AS\QS;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% VISUALIZATION OF Q IN FREQUENCY %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Getting values for the continuous representation of Q %%
wmin=wmean/sqrt(f_ratio); 
wmax=wmean*sqrt(f_ratio); 
xfrq=[wmin:(wmax-wmin)/(10000-1):wmax];

%% P-wave Q continuous values
numP=0;
denP=1;
for j=1:n
    numP=numP+(YP(j)*w_disc(2*j-1)*xfrq(:))./(w_disc(2*j-1)^2+xfrq(:).^2);
    denP=denP-(YP(j)*w_disc(2*j-1).^2)./(w_disc(2*j-1)^2+xfrq(:).^2);
end
Q_contP=denP./numP;

%% S-wave Q continuous values
numS=0;
denS=1;
for j=1:n
    numS=numS+(YS(j)*w_disc(2*j-1)*xfrq(:))./(w_disc(2*j-1)^2+xfrq(:).^2);
    denS=denS-(YS(j)*w_disc(2*j-1).^2)./(w_disc(2*j-1)^2+xfrq(:).^2);
end
Q_contS=denS./numS;

%% Computing fitting quality (RMS and maximum difference)
maxPdif=0;
maxSdif=0;

for j=1:length(Q_contP)
    tempP=abs(Q_contP(j)-QPval);
    if tempP >= maxPdif
        maxPdif=tempP;
    end
    tempS=abs(Q_contS(j)-QSval);
    if tempS >= maxSdif
        maxSdif=tempS;
    end
end

subplot(1,2,1),
semilogx(xfrq/(2*pi),Q_contP,xfrq/(2*pi),QPval*ones(length(xfrq),1)),
title('Q for the P-wave'),legend('computed','desired')
xlabel('frequency (Hz)'),ylabel('Q'),axis([wmin/(2*pi) wmax/(2*pi) 0 QPval*1.25])
subplot(1,2,2),
semilogx(xfrq/(2*pi),Q_contS,xfrq/(2*pi),QSval*ones(length(xfrq),1)),
title('Q for the S-wave'),legend('computed','desired'),
xlabel('frequency (Hz)'),ylabel('Q'),axis([wmin/(2*pi) wmax/(2*pi) 0 QSval*1.25])

disp(strcat('Attenuation bandwidth= [',num2str(wmin/(2*pi)),',',num2str(wmax/(2*pi)),'] Hz'));
disp('');
disp(strcat('Maximum QP fitting error=  ',num2str(maxPdif/QPval*100),' %'));
disp('');
disp(strcat('Maximum QS fitting error=  ',num2str(maxSdif/QSval*100),' %'));
disp(' ');
disp('--------------------------------------')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% VELOCITIES COMPUTATION %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Translating into Y_kappa and Y_mu %%
YK=YP*0;
YM=YP*0;
YK(:)=(vp_0^2*YP(:)-4/3*vs_0^2*YS(:))/(vp_0^2-4/3*vs_0^2);
YM(:)=YS(:);

% Complex modulus
mu_0=vs_0^2*rho;
lam_0=vp_0^2*rho-2*mu_0;

MUK=(lam_0+2/3*mu_0);    % Elastic bulk modulus
MUM=2*mu_0;              % Elastic shear modulus

MK=MUK*ones(length(w),1);
MM=MUM*ones(length(w),1);
lambda=zeros(length(w),1);
mu=zeros(length(w),1);
vp=zeros(length(w),1);
vs=zeros(length(w),1);

for j=1:n
    MK(:)=MK(:)-MUK*(YK(j)*(w_disc(2*j-1)./(w_disc(2*j-1)+i*w(:))));
    MM(:)=MM(:)-MUM*(YM(j)*(w_disc(2*j-1)./(w_disc(2*j-1)+i*w(:))));
end

lambda(:)=MK(:)-1/3*MM(:);     
mu(:)=MM(:)/2;                 

% Complex wave velocities (NOTE: Substitute by vp_0 and vs_0 to get the
% elastic solution)
vp(:)=sqrt((lambda(:)+2*mu(:))/rho);
vs(:)=sqrt(mu(:)/rho);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% PROBLEM FUNCTIONS INITIALIZATION %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Green's functions (CORRECTED, as appear in Carcione (2002))
for j=1:length(w)
    if w(j)>0
        G1=-i*pi/2*(1/vp(j)^2*besselh(0,2,(w(j)*r/(vp(j))))+ ...
            1./(w(j)*r*vs(j))*besselh(1,2,(w(j)*r/(vs(j))))- ...
            1./(w(j)*r*vp(j))*besselh(1,2,(w(j)*r/(vp(j)))));
        G2=i*pi/2*(1/vs(j)^2*besselh(0,2,(w(j)*r/(vs(j))))- ...
            1./(w(j)*r*vs(j))*besselh(1,2,(w(j)*r/(vs(j))))+ ...
            1./(w(j)*r*vp(j))*besselh(1,2,(w(j)*r/(vp(j)))));
        u1(j)=F/(2*pi*rho)*(x*y/r^2)*(G1+G2);
        u2(j)=F/(2*pi*rho)*(1/r^2)*(y^2*G1-x^2*G2);
    end
end

%% Geting the symmetric conjugate of the u1 and u2 functions
for j=1:length(w)
        u1(j)=conj(u1(length(w)+1-j));
        u2(j)=conj(u2(length(w)+1-j));
end

%Ricker type source
S=sqrt(pi)*w.^2/(4*(pi*f0)^3).*exp(-i*w*t0).*exp(-w.^2/(4*(pi*f0)^2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% PROBLEM SOLUTION %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Construction of the w-space displacements
PHI_DX=u1.*S;
PHI_DY=u2.*S;

%% Construction of the w-space velocities
PHI_VX=PHI_DX.*i.*w;
PHI_VY=PHI_DY.*i.*w;

%% Time axis scaling by the Nyquist frequency
dt=1/(2*w_max)*2*pi;        

%% FFT of the PHI functions to obtain solution
Sol_x=(ifft(fftshift(PHI_DX)))/dt;  
Sol_y=(ifft(fftshift(PHI_DY)))/dt;

%% Same for the velocities 
Vel_x=(ifft(fftshift(PHI_VX)))/dt;  
Vel_y=(ifft(fftshift(PHI_VY)))/dt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% VARIABLES' SCALING %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=(1:length(w))*dt-dt ;  %The substraction is to have a t=0 value

%% Rescaling of the time variable to adjust to the t-space sampling rate
for j=1:length(w)/2
    t_res(j)=t(2*j-1);
end

%% Rescaling of the sampling of our function to compare with a t-domain
%% function
for j=1:length(w)/2
    SX_res(j)=Sol_x(2*j-1);
    SY_res(j)=Sol_y(2*j-1);
    VX_res(j)=Vel_x(2*j-1);
    VY_res(j)=Vel_y(2*j-1);
end

%% MATLAB's "ifft" routines don't give us a purely real function for the
%% inverse transform of a conjugate symmetric function. Therefore we have
%% to correct the amplitudes with the small complex residues obtained.
SX_res=abs(SX_res).*sign(real(SX_res));
SY_res=abs(SY_res).*sign(real(SY_res));;
VX_res=abs(VX_res).*sign(real(VX_res));
VY_res=abs(VY_res).*sign(real(VY_res));;

%% Looking for the t value closest to our desired tmax
found=0;
for j=1:length(t_res)
    if found==0
    if t_res(j)>=toutmax
        countt=j;
        found=1;
    end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%  PLOTTING     %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Displacements displayed in desired time window
minx=(min(real(SX_res(1:countt)))-eps)*1.2;
maxx=(max(real(SX_res(1:countt)))+eps)*1.2;
miny=(min(real(SY_res(1:countt)))-eps)*1.2;
maxy=(max(real(SY_res(1:countt)))+eps)*1.2;
figure,
subplot(121),plot(t_res(1:countt),SX_res(1:countt),'-r'),xlabel('Time (sec)'),
legend('x component'),
axis([t_res(1) t_res(countt) minx maxx]),ylabel('x-displacement')
title(strcat('Receiver at x= ',num2str(x),',y= ',num2str(y)))
subplot(122),plot(t_res(1:countt),SY_res(1:countt),'-r'),xlabel('Time (sec)'),
legend('y component'),
axis([t_res(1) t_res(countt) miny maxy]),ylabel('y-dispplacement')
title(strcat('Source frequency= ',num2str(f0),' Hz; Source delay= ',num2str(t0),'s'))

% Velocities displayed in desired time window
minx=(min(real(VX_res(1:countt)))-eps)*1.2;
maxx=(max(real(VX_res(1:countt)))+eps)*1.2;
miny=(min(real(VY_res(1:countt)))-eps)*1.2;
maxy=(max(real(VY_res(1:countt)))+eps)*1.2;
figure,
subplot(121),plot(t_res(1:countt),VX_res(1:countt),'-k'),xlabel('Time (sec)'),
legend('x component'),
axis([t_res(1) t_res(countt) minx maxx]),ylabel('u velocity')
title(strcat('Receiver at x= ',num2str(x),',y= ',num2str(y)))
subplot(122),plot(t_res(1:countt),VY_res(1:countt),'-k'),xlabel('Time (sec)'),
legend('y component'),
axis([t_res(1) t_res(countt) miny maxy]),ylabel('v velocity'),
title(strcat('Source frequency= ',num2str(f0),' Hz; Source delay= ',num2str(t0),'s'))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% SOURCE COMPARISON %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% As a check for the scalings chosen, we can compare the ifft of the
%% Ricker source to the analytical expression of the inverse Fourier
%% transform:
B=pi*pi*f0*f0;
 
%% Real time-domain expression of the Ricker wavelet source
SOrig=t*0;
SOrig=(0.5-B*(t-t0).*(t-t0)).*exp(-B*(t-t0).*(t-t0)); 

SS=ifft(fftshift(S))/dt;

for j=1:length(w)/2
    SS_res(j)=SS(2*j-1);
end

SS_res=abs(SS_res).*sign(real(SS_res));
minS=(min(real(SS_res(1:countt)))-eps)*1.2;
maxS=(max(real(SS_res(1:countt)))+eps)*1.2;

figure,
plot(t(1:countt),SOrig(1:countt),'-b',t_res(1:countt),SS_res(1:countt),'.k'),
title('MATLAB Check: ifft of the source')
xlabel('time (s)'),axis([t_res(1) t_res(countt) minS maxS]),
legend('Analytical source','ifft of the Analytical FT of the source')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% END OF THE PROGRAM !!! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%% 
%% REFERENCES %%
%%%%%%%%%%%%%%%%
% (1) "WAVE PROPAGATION SIMULATION IN A LINEAR VISCOELASTIC MEDIUM", Carcione et al. (1988)
%
% (2) "WAVE FIELDS IN REAL MEDIA: WAVE PROPAGSTION IN ANISOTROPIC, ANELASTIC
% AND POROUS MEDIA", Carcione (2002)
%
% (3) "THE FINITE DIFFERENCE METHOD FOR SEISMOLOGISTS: AN INTRODUCTION",
% Moczo et al. (2005)
%
% (4) "INCORPORATION OF ATTENUATION INTO TIME-DOMAIN COMPUTATIONS OF
% SEISMIC WAVE FIELDS", Emmerich & Korn (1987)