% MATLABScript for obtaning the analytical solution for wave propagation
% in a homogenous poroelastic medium, based on the solution described in 
% the paper by Dai et al. (1995).
%
% The source used is a point explosion with time component Ricker or gaussian first derivative. 
%
% The solution is computed in the frequency domain. Only non-viscous fluids
% are considered.
%
% The resulting seismograms are the radial components of the solid and
% liquid particle velocities.
%
% Material parameters can be customized if "input_type" option 2 is chosen
% (lines 100-110).
%
% The problem setup can be modified by changing lines 42-54.
%
% Different source time functions can be added or modified in lines 205-210. 
%
%FIXIT: The gradient computation might have associated some error
%FIXIT: Signs are not verified!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  By Josep de la Puente, LMU Geophysics, 2006  %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
home
clear all;
close all;

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

% Time and space problem setup
dt=0.0001;                        %Timestep
toutmax=0.25;                     %Maximum outputting time
t=[0:dt:toutmax];                 %Time vector
r=100;                            %Distance source-receiver

% Source parameters
t0=0.04;                         %Source delay
f0=30;                           %Source frequency
Bcoef=pi*pi*f0*f0;

% Frequency domain setup
w_max=10000;                     %Maximum frequency used 
dw=1;                            %Frequency increment

disp(' ')
disp('---------------------------------------------------------------')
disp('The model setup is as follows:')
disp([' *Source-receiver offset = ',num2str(r),' m'])
disp([' *Max. output time       = ',num2str(toutmax),' s'])
disp([' *Timestep               = ',num2str(dt),' s'])
disp([' *Max. frequency in FFT  = ',num2str(w_max),' rad/s'])
disp([' *Freq. step in FFT      = ',num2str(dw),' rad/s'])
disp(' ')
disp('---------------------------------------------------------------')
disp('Which time function do you want to use for ')
sou_time=input('     the source (1:Ricker, 2:first gaussian der.)?');
inputtype=input('Which input do you want to use (1:default, 2:personal)?');
disp(' ')
if sou_time==1
    disp('The source used is a Ricker wavelet with the following characteristics:')
elseif sou_time==2
    disp('The source used is a Gaussian derivative with the following characteristics:')
end
disp([' *Central frequency = ',num2str(f0),' Hz'])
disp([' *Time delay = ',num2str(t0),' s'])
disp(' ')
disp('---------------------------------------------------------------')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% INPUT MATERIAL PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if inputtype==1 %After (Dai et al. 1995). Note that they use a different time source!
    %Poroelastic parameters
    Poro=0.1;
    q=0.1;

    %Biot's elastic constants
    P=2.0332e10;
    Q=0.0953e10;
    N=0.684e10;
    R=0.0331e10;

    % Phases inertia
    rho(1,1)=2167;
    rho(2,2)=191;
    rho(1,2)=-83;
elseif inputtype==2    
    %Poroelastic parameters (S.I. units)
    Poro=0.5;                           %Porosity
    nu=0;                               %Viscosity (never use non-zero values!)
    kappa=300e-15;                      %Permeability
    rhoS=2500;                          %Solid's density
    rhoF=1040;                          %Fluid's density
    Tor=2;                              %Tortuosity
    lam=19.6e9;                         %Elastic lambda of the solid matrix
    mu= 26.1e9;                         %Elastic mu of the solid matrix
    KS=80e9;                            %Solid's bulk modulus
    KF=2.5e9;                           %Fluid's bulk modulus
    
    % Derived parameters
    alpha1=1-(3*lam+2*mu)/(3*KS);       %Biot's effective stress coefficient
    KK=lam+2/3*mu;                      %Generalized drained modulus
    MM=KS/((1-KK/KS)-Poro*(1-KS/KF));   %Fluid-solid coupling modulus
    q=0.5;                              %Ratio of force distribution
 
    %Biot's elastic constants
    P=((1-Poro)*(1-Poro-KK/KS)+Poro*KK/KF)*MM+4/3*mu;
    Q=(1-Poro-KK/KS)*Poro*MM;
    N=mu;
    R=(Poro^2*KS)/(1-Poro-KK/KS+Poro*KS/KF);

    % Phases inertia
    rho(2,2)=Tor*Poro*rhoF;
    rho(1,2)=-Poro*rhoF*(Tor-1);       
    rho(1,1)=(1-Poro)*rhoS-rho(1,2);

end

disp('The material chosen has the following Biot.s elastic constants:')
disp(['  P = ',num2str(P)])
disp(['  Q = ',num2str(Q)])
disp(['  N = ',num2str(N)])
disp(['  R = ',num2str(R)])
disp(' ')
disp('and the following phases inertia:')
disp(['  rho_{11} = ',num2str(rho(1,1))])
disp(['  rho_{22} = ',num2str(rho(2,2))])
disp(['  rho_{12] = ',num2str(rho(1,2))])
disp(' ')
disp('---------------------------------------------------------------')
disp(' ')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% DERIVED PARAMETER COMPUTATION AND INITIALIZATION %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% D
clear i;
for k=1:2
    for j =k:2
        D(k,j)=rho(k,j); 
    end
end

% Damping
b= 0; 

% Velocity computation
A=D(1,1)*D(2,2)-D(1,2)^2;
B=D(1,1)*R+D(2,2)*P-2*D(1,2)*Q;
C=P*R-Q^2;

VS=sqrt((B-sqrt(B^2-4*A*C))/(2*A));
VF=sqrt((B+sqrt(B^2-4*A*C))/(2*A));

AF=-(D(1,2)*VF^2-Q)/(D(2,2)*VF^2-R);
AS=-(D(1,2)*VS^2-Q)/(D(2,2)*VS^2-R);

Alpha=((1-q)*(Q+R*AS)-q*(P+Q*AS))/((P*R-Q^2)*(AF-AS));
Beta= ((q)*(P+Q*AF)-(1-q)*(Q+R*AF))/((P*R-Q^2)*(AF-AS));

%Potential calculation
PHIS=0*t;
PHIF=0*t;
PHIS_min=0*t;
PHIF_min=0*t;
PHIS_max=0*t;
PHIF_max=0*t;
DS=0*t;
DF=0*t;
VelS=0*t;
VelF=0*t;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% FREQUENCY DOMAIN COMPUTATION %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

w=[-w_max:dw:w_max];
for j=1:length(w)
    if w(j)==0
        zero_w=j;
        count=1;
    end
end

kF=w/VF;
kS=w/VS;
clear i;
dt2=1/(2*w_max)*2*pi;

% Source choice
if sou_time==1 %% Ricker (second derivative of gaussian)
    SFREQ=sqrt(pi)*w.^2/(4*(pi*f0)^3).*exp(-i*w*t0).*exp(-w.^2/(4*(pi*f0)^2));
elseif sou_time==2 %% First derivative of gaussian
    SFREQ=sqrt(pi)*w/(4*(pi*f0)^3).*exp(-i*w*t0).*exp(-w.^2/(4*(pi*f0)^2)); 
end

% Gradient computation initialization
rinc=0.00001*r;
r_min=r-rinc;
r_max=r+rinc;
rvec=[r_min (r-r_min)/2 (r_max-r)/2 r_max];    

StoreS=zeros(3,length(w));
StoreF=zeros(3,length(w));

for j=1:4
    rloc=rvec(j);
    PHISF=SFREQ.*(Alpha*exp(-i*kF*rloc)+Beta*exp(-i*kS*rloc))/(4*pi*rloc);
    PHIFF=SFREQ.*(Alpha*AF*exp(-i*kF*rloc)+Beta*AS*exp(-i*kS*rloc))/(4*pi*rloc);
    StoreS(j,:)=ifft(fftshift(PHISF.*i.*w))/dt;
    StoreF(j,:)=ifft(fftshift(PHIFF.*i.*w))/dt;
end

%Velocity computation through gradient (Taylor operator, 4-points)
velS=(StoreS(1,:)/6-StoreS(2,:)*4/3+StoreS(3,:)*4/3-StoreS(4,:)/6)/(rinc);
velF=(StoreF(1,:)/6-StoreF(2,:)*4/3+StoreF(3,:)*4/3-StoreF(4,:)/6)/(rinc);

for j=1:length(w)/2
    VS_res(j)=velS(2*j-1);
    VF_res(j)=velF(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.
VS_res=abs(VS_res).*sign(real(VS_res));
VF_res=abs(VF_res).*sign(real(VF_res));;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% TIME DOMAIN SOLUTION PLOTTING %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t2=(1:length(w))*dt2-dt2 ;  

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

%% 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

% Fluid velocities w.r.t to the solid
VR_res=Poro*(VF_res-VS_res);

% Velocities displayed in desired time window
minx=(min(real(VS_res(1:countt))))*1.2;
maxx=(max(real(VS_res(1:countt))))*1.2;
miny=(min(real(VR_res(1:countt))))*1.2;
maxy=(max(real(VR_res(1:countt))))*1.2;
figure,
subplot(121),plot(t_res(1:countt),VS_res(1:countt),'-k'),xlabel('Time (sec)'),
legend('Solid velocity'),
axis([t_res(1) t_res(countt) minx maxx]),ylabel('radial velocity (m/s)')
title(strcat('Receiver at r= ',num2str(r)))
subplot(122),plot(t_res(1:countt),VR_res(1:countt),'-k'),xlabel('Time (sec)'),
legend('Fluid velocity (w.r.t solid)'),
axis([t_res(1) t_res(countt) miny maxy]),ylabel('radial velocity (m/s)'),
title(strcat('Source frequency= ',num2str(f0),' Hz; Source delay= ',num2str(t0),'s'))


disp(' ')
disp('Computation successfully finished!')