%w为角频率，lamda，u为弹性参数，fai为孔隙度，e为裂隙密度，tao为时间尺度参数，en为流体粘滞系数，v,miu为泊松比
%r为裂隙纵横比，xiaok为渗透率，cv为裂隙体积,a为裂隙半径，kf为流体体积模量
clear all;
vp=2790;vs=1463;rou=2.08;fai=0.3;tao=2*10^-5;e=0.1;miu=0.3104;v=miu;xiaok=1000;
lamda=rou*vp*vp-2*rou*vs*vs;r=0.1;
u=rou*vs*vs;kf=2E9;
k=rou*(vp^2-4*v*vs/3);
kc= (pi*u*r*(lamda + u))/(kf*(lamda + 2.0*u));
u=rou*vs*vs;
%k=2000;
en=0.01;a=0.1;l=2;
f=5:1:104;
%以下为输入的值；
w=2*pi*f;
cv=4*pi*a*a*a*r/3;
sigma=tao*6*xiaok*l/(en*cv*(1+kc));
kf=(pi*u*u*r*r)/(2*(1-v)*sigma);
%tao=4*en*a*a*a*(1-v)/(9*xiaok*l*u);%裂隙纵横比很小的时候
kp=4*u/(3*kf);
y=(3*pi*(lamda+u)*(1+kp))/(4*(lamda+2*u)*(1+kc));
yyipie=y*(lamda+2*u)/((3*lamda+2*u)*(1+kp));
awfenzi=((ones(1,100)+i*w*y*tao)*(lamda+2*u)/(lamda+u))*(16*e/(27*fai*(1+kp))+(lamda+u)/(3*lamda+2*miu))+i*w*tao*(1/(3*(1+kc))-yyipie);
awfenmu=1+i*w*tao+(1+i*w*y*tao)*(16*e*(1+kc)*(lamda+2*u))/(9*fai*(1+kp)*(lamda+u));
bwfenzi=1/(3*(1+kc))+9*(1+kp)*(lamda+u)/(16*(1+kc)*(3*lamda+2*miu))+i*w*tao*(yyipie-1/(3*(1+kc)))/(1+i*w*pi);
bwfenmu=9*(1+kp)*(lamda+u)/(16*(1+kc)*(lamda+2*u))+(1+i*w*y*tao)/(1+i*w*tao);
aw=awfenzi/awfenmu;
bw=bwfenzi/bwfenmu;
keff=k-e*(4*(3*lamda+2*miu)*(lamda+2*miu)*(1-aw)/(miu*(lamda+miu))-4*pi*aw*r)-fai*((3*lamda+2*miu)*((lamda+2*miu)/(3*lamda+2*miu)+bw)/4*miu-3*bw);
ueff=u*ones(1,100)-((16/45)*e/(1+kc))*(u*(lamda+2*u)/(3*lamda+4*u))*(kc+1./(1+i*w*tao))-ones(1,100)*(32/45)*e*(u*(lamda+2*u)/(3*lamda+4*u))-ones(1,100)*(fai*15*u*(lamda+2*u)/(9*lamda+14*u));
eeff=keff+4*ueff/3;
qp=(-real(eeff))./(imag(eeff));
qs=(-real(ueff))./(imag(ueff));
x=log(abs(f))/log(20);
figure(1);
qpni=1./qp;
plot(x,qpni);
qsni=1./qs;
figure(2);
plot(x,qsni);
