Add files via upload

Matlab Relaxation Time/Emmerich.m

@@ -0,0 +1,36 @@
+function [tau_eta, tau_sig]=Emmerich(f_ratio,QPval,wmin,wmax, xfrq, N)
+K=2*N-1;
+
+theta = wmin*(wmax/wmin).^([0:N-1]/(N-1));
+wk=wmin*(wmax/wmin).^([0:K-1]/(K-1));
+QPval_K=interp1(xfrq, QPval, wk);
+over_Q=1./QPval_K;
+%%
+AK=zeros(K,N);
+QB=zeros(K,1);
+QB(:)=over_Q;
+for k=1:K
+ for l=1:N
+ AK(k,l)=(wk(k)*(theta(l)-wk(k)*QB(k)))/(theta(l)^2+wk(k)^2);
+ end
+end
+kl=AK\QB;
+%%
+Q_freq1=zeros(1,size(xfrq,2));
+Q_freq2=zeros(1,size(xfrq,2));
+for l=1:N
+ Q_freq1=Q_freq1+xfrq.*theta(l)*kl(l)./(theta(l)^2+xfrq.^2);
+ Q_freq2=Q_freq2+xfrq.^2*kl(l)./(theta(l)^2+xfrq.^2);
+end
+Q_f=Q_freq1./(1+Q_freq2);
+
+%figure(3),hold on,plot(log10(xfrq/(2*pi)),ones(length(xfrq),1)./QPval)
+%figure(3),semilogx(wk./(2*pi),QB, xfrq/(2*pi),Q_f, '--') 
+%figure(3),hold on,semilogx(log10(xfrq/(2*pi)),Q_f,'b')
+%legend('Q-ref', 'Yang', 'Emmerich')
+%xlabel('log10(f)'); ylabel('1/Q')
+
+%er=sum(1./QPval-Q_f).^2
+%%
+tau_sig=1./theta;
+tau_eta=(1+N.*kl')./theta;

Matlab Relaxation Time/QpDraw.m

@@ -0,0 +1,18 @@
+function Q_f=QpDraw(QPval, xfrq, tau_ep,tau_sig)
+N=length(tau_ep);
+K=4*N;
+wmin=xfrq(1); wmax=xfrq(end);
+wk=wmin*(wmax/wmin).^([0:K-1]/(K-1));
+QPval_K=interp1(xfrq, QPval, wk);
+theta = 1./tau_sig;
+kl=(tau_ep./tau_sig-1)./N;
+Q_freq1=zeros(1,size(xfrq,2));
+Q_freq2=zeros(1,size(xfrq,2));
+
+for l=1:N
+ Q_freq1=Q_freq1+xfrq.*theta(l)*kl(l)./(theta(l)^2+xfrq.^2);
+ Q_freq2=Q_freq2+xfrq.^2*kl(l)./(theta(l)^2+xfrq.^2);
+end
+QB=zeros(K,1);
+QB(:)=1./QPval_K;
+Q_f=Q_freq1./(1+Q_freq2);

Matlab Relaxation Time/RelaxationTime.m

@@ -0,0 +1,71 @@
+% RelaxationTime script calculates the relaxation time of stress and strain
+% for input file: relaxation_time.inp
+% Methodlogy is based on Emmerich and Korn, Incorporation of attenuation
+% into time-domain computations of seismic wave fields, Geophysics, 1987,
+% 52:9, 1252-1264.
+% Strain Relaxation Time: tau_eta_e
+% Stress Relaxation Time: tau_sig_e
+
+clc; clear all
+freq=input('Central Frequency ?');
+f_ratio=101;
+fmin0=1;
+fmin=exp(log(freq)-log(12)/2);
+fmax=12*fmin;
+wmin=2*pi*fmin;%/sqrt(f_ratio);
+wmax=2*pi*fmax;%*sqrt(f_ratio); 
+xfrq=[wmin:(wmax-wmin)/(10000-1):wmax]; %w
+x_axis=xfrq/(2*pi);
+t=[0:10^-4:0.5];
+Legend=cell(1,1); ii=0;
+%%
+%Constant Q
+QPval=zeros(1,length(xfrq));
+Q=input('Attenuation factor Q ?');
+QPval(:)=Q;
+%%
+ns=input('How many springs ?');
+nsn=size(ns,2);
+y2=zeros(size(ns,2),size(t,2));
+ii=0;
+for N=ns
+ ii=ii+1;
+ Legend{ii}=strcat(num2str(N));
+
+ [tau_eta_e, tau_sig_e]=Emmerich(f_ratio,QPval,wmin,wmax, xfrq, N);
+ ex_e=zeros(N,size(t,2));
+ ex_e=exp(-t./tau_sig_e');
+ ex_e=((1-tau_eta_e./tau_sig_e))./N*ex_e; %N ->sum
+ x=t;
+ y=2.25*10^9*(1-ex_e);
+ y2(ii,:)=y;
+
+ figure(2),hold on,plot(t,y)
+ pbaspect([2 1 1])
+end
+
+ylabel('Bulk Modulus (GPa)')
+xlabel('t(s)')
+
+y_log=log(y2);
+theta_e=1./tau_sig_e;
+kl_e=(tau_eta_e./tau_sig_e-1)/N;
+%%
+Qval=zeros(1,length(xfrq));
+Qval(:)=QPval;
+save('Qval.mat','xfrq','Qval')
+x=[sqrt(abs(theta_e)) sqrt(abs(kl_e))];
+
+%%
+Q_f1=QpDraw(QPval, xfrq, tau_eta_e,tau_sig_e);
+
+%%
+figure(4),semilogx(xfrq/(2*pi),1./Qval, 'k')
+hold on
+figure(4),semilogx(xfrq/(2*pi),Q_f1, '--')
+hold on
+
+xlabel('log10(f)'); ylabel('1/Q')
+legend('1/Qref','SolvOpt')
+
+pbaspect([2 1 1])

Matlab Relaxation Time/relaxation_time.inp

@@ -0,0 +1,5 @@
+VISCO | RECURSIVE | Nhete !Fc=15
+1 1 1
+1 3
+3.5774296952513399E-002 6.7481038667282700E-003 1.6313910723490001E-003
+2.2805059493937799E-002 4.7530094027435996E-003 8.9270249012923103E-004