clc
clear
MkSize=12;
LdWidth=1.5;
FtSize=20;
FtName='Arial';
FtSize2=16;
%%%%%% Fontain sandstone laboratory data
fm_3_5MPa=[0.00697 0.00999 0.01969 0.0396 0.0679 0.09728 0.111878 0.19966 0.39385 0.7033 0.98807];
fm_5MPa=[0.00697 0.00999 0.01969 0.0396 0.0679 0.09728 0.111878 0.19966 0.39385 0.7033 0.98807];
fm_10MPa=[0.00695 0.00994 0.00196 0.0395 0.0691544 0.1073 0.1995 0.3937];
fm_20MPa=[0.006947 0.00955 0.019619 0.03949 0.069 0.097];
K_3_5MPa=[19.56 19.78 21.43 22.5275 24.5055 26.1538 26.7033 30.5495 34.3956 37.5824 39.3407]*1e9;
K_5MPa=[30.769 31.97 34.0659 36.4835 38.2418 38.9011 39.6308 40.329 40.879 40.3297 40]*1e9;
K_10MPa=[37.4725 42.8571 0 42.6374 42.63 43.5165 42.1978 41.4286]*1e9;
K_20MPa=[41.32 44.62 44.18 43.74 43.956 43.516]*1e9;
QK_3_5MPa=[0.1105 0.1526 0.1539 0.2 0.2368 0.2605 0.2578 0.2605 0.215789 0.14868 0.071];
QK_5MPa=[0.133 0.0899 0.129 0.1213 0.1081 0.0939 0.08678 0.04984 -0.021 -0.0789 -0.07894];


%% %%% Basic parameters for modelling
Kgr =    77e9;
mu_gr =  33.4e9;
fai_s0 = 0.139;   %
etagly = 1;
Kf =   4.36e9;
rho_f =  1261;
k0 =     5*10e-17;
at = 1.6;
K_drs =  37.35e9;
mu_drs = 21.65e9;
lamda_drs = K_drs-2*mu_drs/3;
L_drs = lamda_drs+mu_drs;
vdrs = (3*K_drs-2*mu_drs)/(6*K_drs+2*mu_drs);
j = sqrt(-1);
fai_c0 = 1.8e-4;
theta_c = 4070;
Cdrs = 1/K_drs;
Cgr = 1/Kgr;
rho_total = 2191;

theta_cK = -150000e9;      %% control vp_fast
theta_cmu = -2500e9;  
% theta_cL = -700e9;
theta_cpasai1s = 0e9;
theta_cpsai2s = 0e9;
psai_1s = -765e9;
psai_2s = -710e9;       
gamma = -440e9;
gamma1 = 0e9;               %%作用不大
gamma2 = 150e9;
gamma3 =  45e9;

yy = 1;
for f = -2:0.01:20
    
    w = 2*pi*10^(f-1);
    
    arf = 0.0001;
    
    i = 1;
    for pc = (10:5:95)*1e6
        pf = 5e6;
        p = pc-pf;
        
%% Kf_shear relaxation
%  Gx = j*w*etagly;
%  r_f = sqrt((36*Gx)/(3*Kf+4*Gx));  
%  Kfx = Kf+4/3*Gx-(Kf-2/3*Gx)^2./((Kf+4/3*Gx)*r_f/(2*arf)...
%             *besseli(0,r_f/(2*arf))/(2*besseli(1,r_f/(2*arf)))-Gx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Kf____squirt flow+shear relaxation
ka = 1/arf*sqrt(3*j*w*etagly/Kf);
Kf_x = (1-2*besseli(1,ka)./(ka*besseli(0,ka)))*Kf;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        rho_wave = at.*rho_f./(fai_s0-p.*(Cdrs-Cgr))-j*etagly./(k0.*w);
        
        varib = 1/(K_drs+theta_cK.*fai_c0*exp(-theta_c.*p.*Cdrs))-1/K_drs;

        Kmf(yy,i) = 1./(1/K_drs+1./(1./varib...
            +1./((1/Kf_x-1/Kgr)*fai_c0*exp(-theta_c*p*Cdrs))));
        
        mu_mf(yy,i) = 1./(1/(mu_drs+theta_cmu.*fai_c0*exp(-theta_c.*p.*Cdrs))-(4/15)*(1/(K_drs+theta_cK.*fai_c0*exp(-theta_c.*p.*Cdrs))-1/Kmf(yy,i)));

        
        lamda_mf(yy,i) = Kmf(yy,i)-(2/3)*mu_mf(yy,i);
        
        n = 1-Kmf(yy,i)/Kgr;
        
        M  = 1./((n-fai_s0+p.*(Cdrs-Cgr))./Kmf(yy,i)+(fai_s0-p.*(Cdrs-Cgr))./Kf);
        es = (pc-(fai_s0-(pc-pf).*(Cdrs-Cgr))*pf)./((3-3*fai_s0+3*(pc-pf)*(Cdrs-Cgr)).*Kgr);
        theta = (pc-(1-Kmf(yy,i)./Kgr).*pf)./Kmf(yy,i);
        ef = (theta-3*(1-fai_s0+(pc-pf).*(Cdrs-Cgr)).*es)./(3*fai_s0-3.*(pc-pf).*(Cdrs-Cgr));
        
        E1313 = mu_mf(yy,i)-(n*M+gamma/2)*ef...
            +(3*lamda_mf(yy,i)+4*mu_mf(yy,i)+psai_1s+theta_cpasai1s*fai_c0...
            *exp(-theta_c*p*Cdrs)+3*n^2*M)*es;
        E3333 = lamda_mf(yy,i)+2*mu_mf(yy,i)+n.^2.*M...
            -(n.*M-2*gamma2).*ef+(5*lamda_mf(yy,i)+6*mu_mf(yy,i)+psai_2s...
            +theta_cpsai2s.*fai_c0.*exp(-theta_c.*p.*Cdrs)+5*n^2.*M).*es;
        R33 = n.*M+(n.*M-2*gamma-6*gamma2).*es-(2*gamma3+n.*M).*ef;
        R = M+(6*gamma3+3*n.*M).*es-(3*M-6*gamma1).*ef;
        
        vp_fast(yy,i) = sqrt((rho_total.*R+E3333.*rho_wave-2*rho_f.*R33)./(2*(rho_total.*rho_wave-rho_f.^2))...
            .*(1+sqrt(1-(4*(rho_total.*rho_wave-rho_f.^2).*(R.*E3333-R33.^2))./(rho_total.*R...
            +rho_wave.*E3333-2*rho_f.*R33).^2)));
        vs(yy,i) = sqrt(E1313./(rho_total-rho_f.^2/rho_wave));
        i = i+1;
        
    end
    yy = yy+1;
end

Kw_pre = rho_total.*(vp_fast.^2-(4/3)*vs.^2);
muw_pre = rho_total.*vs.^2;

aK = imag(Kw_pre)./(real(Kw_pre));  %%% 1/Q
aG = imag(muw_pre)./(real(muw_pre));  %%% 1/Q


f = -3:0.01:19;
figure(1)
set(gcf,'Position',[500 150 800 600])
plot(f,Kw_pre(:,1)/1e9,'k-','linewidth',1);
hold on
plot(f,Kw_pre(:,2)/1e9,'k--','linewidth',LdWidth);
hold on
plot(f,Kw_pre(:,4)/1e9,'k-.','linewidth',LdWidth);
hold on
plot(log10(fm_5MPa),K_5MPa/1e9,'k^','MarkerSize',MkSize,'MarkerEdgecolor','k','MarkerFacecolor','k');
hold on
plot(log10(fm_10MPa),K_10MPa/1e9,'ks','MarkerSize',MkSize+2,'MarkerEdgecolor','k');
hold on
plot(log10(fm_20MPa),K_20MPa/1e9,'kd','MarkerSize',MkSize,'MarkerEdgecolor','k','MarkerFacecolor','k');
xlabel('Log[frequency (Hz)]','FontName',FtName,'FontSize',FtSize);
ylabel('Bulk modulus (GPa)','FontName',FtName,'FontSize',FtSize);
axis([-3 1 10 50])
set(gca,'xtick',-3:1:12,'ytick',10:10:50,'FontName',FtName,'FontSize',FtSize)
lgd=legend('Predict,  5  MPa','Predict, 10 MPa','Predict, 20 MPa',...
    'Measure,  5  MPa','Measure, 10 MPa','Measure, 20 MPa','Location','Southeast','NumColumns',2);
legend('boxoff')
title(lgd,'Indiana limestone \phi_{\its} = 13.9%','FontName',FtName,'FontSize',FtSize,'FontWeight','Normal')


figure(2)
set(gcf,'Position',[500 150 800 600])
plot(f,aK(:,1),'k-','linewidth',1);
hold on
plot(f,aK(:,2),'k--','linewidth',LdWidth);
hold on
plot(f,aK(:,4),'k-.','linewidth',LdWidth);
hold on
plot(log10(fm_5MPa),QK_5MPa,'k^','MarkerSize',MkSize,'MarkerEdgecolor','k','MarkerFacecolor','k');
xlabel('Log[frequency (Hz)]','FontName',FtName,'FontSize',FtSize);
ylabel('Bulk mudulus attenuation','FontName',FtName,'FontSize',FtSize);
axis([-3 1 0 0.2])
set(gca,'xtick',-3:1:12,'ytick',0:0.05:0.2,'FontName',FtName,'FontSize',FtSize)
lgd=legend('Predict,   5  MPa','Predict,  10 MPa','Predict,  20 MPa',...
    'Measure, 5 MPa','Location','Northeast');
legend('boxoff')
title(lgd,'Indiana limestone \phi_{\its} = 13.9%','FontName',FtName,'FontSize',FtSize,'FontWeight','Normal')
