clc
clear
MkSize=12;
LdWidth=2;
FtSize=20;
FtName='Arial';
FtSize2=16;
%%%%%% Vos sandstone%%%%%%%%%%
v = xlsread('C:\Users\DELL\Desktop\桌面文献及材料\Fontainebleau-Vosges sandstone data\Vosges\Vosges_data.xlsx');
vpw_linear = xlsread('vp_Vos_linear.xls');
vs_linear =  xlsread('vs_Vos_linear.xls');
Kw_linear =  xlsread('Kw_Vos_linear.xls');
muw_linear = xlsread('Muw_Vos_linear.xls');

Kw_DZ =  xlsread('K_Vos_DZ_model.xls');
muw_DZ = xlsread('Mu_Vos_DZ_model.xls');

p = v(:,1);     %%% pressure
%%%%%%%%%%%  vp    %%%%%%%%%%%%
vpdry_lab = v(:,3);  
vpw_lab = v(:,2);
%%%%%%%%%%%  vs    %%%%%%%%%%%%
vsdry_lab = v(:,5); 
vsw_lab = v(:,4);
%%%%%%%%% density %%%%%%%%%%%
den_d = 1950;
rho_total = 2130;
%% DZ_velocity
vsw_DZ = sqrt(muw_DZ(:,121)./rho_total);
vpw_DZ = sqrt((Kw_DZ(:,121)+(4/3).*muw_DZ(:,121))./rho_total); 

%%%%%%%%% moduli_dry_g_w %%%%%%%%%%%
Kdry_lab = den_d.*(vpdry_lab.^2-(4/3)*vsdry_lab.^2);
Kw_lab = rho_total.*(vpw_lab.^2-(4/3)*vsw_lab.^2);
mudry_lab = den_d.*vsdry_lab.^2;
muw_lab = rho_total.*vsw_lab.^2;


%% %%% Basic parameters for modelling
Kgr = 39e9;
mu_gr = 24e9;
fai_s0 =0.25;   %
etagly = 1e-3;
Kf = 2.24e9;
rho_f = 1000;
k0 = 1*10e-16;
at = 1.8;
K_drs = 13.2e9;
mu_drs = 10.6e9;
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);
pos_drs = (3*K_drs-2*mu_drs)./(6*K_drs+2*mu_drs);
j = sqrt(-1);
fai_c0 = 0.00125;
theta_c = 1070;
Cdrs = 1/K_drs;
Cgr = 1/Kgr;
f = 1e6;
w = 2*pi*f;
rho_total = 2130;

theta_cK = -7300e9;      %% control vp_fast
theta_cmu = -5500e9;  
theta_cpasai1s = 0e9;
theta_cpsai2s = 0e9;
psai_1s = -850e9;
psai_2s = -700e9;       
gamma = -360e9;
gamma1 = 0e9;               %%no effect
gamma2 = 130e9;
gamma3 =  75e9;

f = 1e6;
w = 2*pi*f;
p_close = 120e6;
arf0 = (4-4*pos_drs^2)./(3*pi*K_drs*(1-2*pos_drs))*p_close;

i = 1;
for pc = (5:5:115)*1e6
pf = 5e6;
p = pc-pf;
arf(i) = arf0-(4-4*pos_drs^2)./(3*pi*K_drs*(1-2*pos_drs))*p;

rho_wave = at.*rho_f./(fai_s0-p.*(Cdrs-Cgr))-j*etagly./(k0.*w);
% x = 1/K_drs*(1+theta_cK.*fai_c0*exp(-theta_c.*p.*Cdrs))

varib = 1/(K_drs+theta_cK.*fai_c0*exp(-theta_c.*p.*Cdrs))-1/K_drs;

ka = sqrt(3*j*w*etagly/Kf)/arf0;
Kf_x = (1-2*besseli(1,ka)/(ka*besseli(0,ka)))*Kf;

Kmf(i) = 1./(1./K_drs+1./(1/varib...
    +1/((1/Kf_x-1/Kgr)*fai_c0*exp(-theta_c*p*Cdrs))));
mu_mf(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(i)));

% kk = K_drs+theta_cK.*fai_c0*exp(-theta_c.*p.*Cdrs);
% cc = mu_drs+theta_cmu.*fai_c0*exp(-theta_c.*p.*Cdrs);

lamda_mf(i) = Kmf(i)-(2/3)*mu_mf(i);
n = 1-Kmf(i)/Kgr;

M  = 1./((n-fai_s0+p.*(Cdrs-Cgr))./Kmf(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(i)./Kgr).*pf)./Kmf(i);
ef = (theta-3*(1-fai_s0+(pc-pf).*(Cdrs-Cgr)).*es)./(3*fai_s0-3.*(pc-pf).*(Cdrs-Cgr));

E1313 = mu_mf(i)-(n*M+gamma/2)*ef...
    +(3*lamda_mf(i)+4*mu_mf(i)+psai_1s+theta_cpasai1s*fai_c0...
    *exp(-theta_c*p*Cdrs)+3*n^2*M)*es;
E3333 = lamda_mf(i)+2*mu_mf(i)+n.^2.*M...
    -(n.*M-2*gamma2).*ef+(5*lamda_mf(i)+6*mu_mf(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(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(i) = sqrt(E1313./(rho_total-rho_f.^2/rho_wave));
i = i+1;

end

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

p1 = 0:5:110;
pp = 0:0.1:120;
figure(1)
set(gcf,'Position',[500 150 800 600])
plot(p1,vpw_lab/1e3,'ko','MarkerSize',10,'MarkerFaceColor','k');
hold on
plot(p1,vpdry_lab/1e3,'ko','MarkerSize',10);
hold on
plot(p1,vpw_linear/1e3,'k--','linewidth',LdWidth);
hold on
plot(pp,vpw_DZ/1e3,'m--','linewidth',LdWidth);
hold on
plot(p1,vp_fast/1e3,'k-','linewidth',LdWidth);
xlabel('Effective pressure (MPa)','FontName',FtName,'FontSize',FtSize);
ylabel('P-wave velocity (km/s)','FontName',FtName,'FontSize',FtSize);
axis([0 110 2.2 4])
set(gca,'xtick',0:10:110,'ytick',2.2:0.3:4,'FontName',FtName,'FontSize',FtSize)
lgd=legend(' Measured Vp, saturated',' Measured Vp, dry',' Predicted  Vp, saturated, Liu model',...
    ' Predicted  Vp, saturated, DZ model',' Predicted  Vp, saturated, new model','Location','Southeast');
legend('boxoff')
title(lgd,'Vosges \phi_{\its} = 25%','FontName',FtName,'FontSize',FtSize,'FontWeight','Normal')

figure(2)
set(gcf,'Position',[500 150 800 600])
plot(p1,vsw_lab/1e3,'ko','MarkerSize',10,'MarkerFaceColor','k');
hold on
plot(p1,vsdry_lab/1e3,'ko','MarkerSize',10);
hold on
plot(p1,vs_linear/1e3,'k--','linewidth',LdWidth);
hold on
plot(pp,vsw_DZ/1e3,'m--','linewidth',LdWidth);
hold on
plot(p1,vs/1e3,'k-','linewidth',LdWidth);
xlabel('Effective pressure (MPa)','FontName',FtName,'FontSize',FtSize);
ylabel('S-wave velocity (km/s)','FontName',FtName,'FontSize',FtSize);
axis([0 110 1.4 2.4])
set(gca,'xtick',0:10:110,'ytick',1.4:0.2:2.4,'FontName',FtName,'FontSize',FtSize)
lgd=legend(' Measured Vs, saturated',' Measured Vs, dry',' Predicted  Vs, saturated, Liu model',...
    ' Predicted  Vs, saturated, DZ model',' Predicted  Vs, saturated, new model','Location','Southeast');
legend('boxoff')
title(lgd,'Vosges \phi_{\its} = 25%','FontName',FtName,'FontSize',FtSize,'FontWeight','Normal')


figure(3)
set(gcf,'Position',[500 150 800 600])
plot(p1,Kw_lab/1e9,'ko','MarkerSize',10,'MarkerFaceColor','k');
hold on
plot(p1,Kdry_lab/1e9,'ko','MarkerSize',10);
hold on
plot(p1,Kw_linear/1e9,'k--','linewidth',LdWidth);
hold on
plot(pp,Kw_DZ(:,121)/1e9,'m--','linewidth',LdWidth);
hold on
plot(p1,Kw_pre/1e9,'k-','linewidth',LdWidth);
xlabel('Effective pressure (MPa)','FontName',FtName,'FontSize',FtSize);
ylabel('Bulk modulus (GPa)','FontName',FtName,'FontSize',FtSize);
axis([0 110 4 18])
set(gca,'xtick',0:10:110,'ytick',4:2:18,'FontName',FtName,'FontSize',FtSize)
lgd=legend(' Measured K, saturated',' Measured K, dry',' Predicted  K, saturated, Liu model',...
    ' Predicted  K, saturated, DZ model',' Predicted  K, saturated, new model','Location','Southeast');
legend('boxoff')
title(lgd,'Vosges \phi_{\its} = 25%','FontName',FtName,'FontSize',FtSize,'FontWeight','Normal')


figure(4)
set(gcf,'Position',[500 150 800 600])
plot(p1,muw_lab/1e9,'ko','MarkerSize',10,'MarkerFaceColor','k');
hold on
plot(p1,mudry_lab/1e9,'ko','MarkerSize',10);
hold on
plot(p1,muw_linear/1e9,'k--','linewidth',LdWidth);
hold on
plot(pp,muw_DZ(:,121)/1e9,'m--','linewidth',LdWidth);
hold on
plot(p1,muw_pre/1e9,'k-','linewidth',LdWidth);
xlabel('Effective pressure (MPa)','FontName',FtName,'FontSize',FtSize);
ylabel('Shear modulus (GPa)','FontName',FtName,'FontSize',FtSize);
axis([0 110 4 12])
set(gca,'xtick',0:10:110,'ytick',4:2:12,'FontName',FtName,'FontSize',FtSize)
lgd=legend(' Measured \mu, saturated',' Measured \mu, dry',' Predicted  \mu, saturated, Liu model',...
    ' Predicted  \mu, saturated, DZ model',' Predicted  \mu, saturated, new model','Location','Southeast');
legend('boxoff')
title(lgd,'Vosges \phi_{\its} = 25%','FontName',FtName,'FontSize',FtSize,'FontWeight','Normal')



figure(5)
set(gcf,'Position',[500 150 800 600])
plot(p1,arf,'k--','MarkerSize',10,'MarkerFaceColor','k');
xlabel('Effective pressure (MPa)','FontName',FtName,'FontSize',FtSize);
ylabel('Crack aspect ratio (-)','FontName',FtName,'FontSize',FtSize);
axis([0 110 0 6e-3])
set(gca,'xtick',0:10:110,'ytick',0:1e-3:6e-3,'FontName',FtName,'FontSize',FtSize)
lgd=legend(' ','Location','northeast');
legend('boxoff')
title(lgd,'Fontainebleau \phi_{\its} = 4%','FontName',FtName,'FontSize',FtSize,'FontWeight','Normal')


