Dear mecej4,
Thank you so much for your quick reply,
In these picture, we will see all files used in this code, in the same time the error I got when I tried to run it.
Here is the mlmcodaq-sn.f file
c f77 mlmcodaq-sn.f butfilter.f four1.f /usr/local/sac/lib/libsac.a /usr/local/sac/lib/libsacio.a -o mlmcodaq-sn
c pgf77 -tp k8-32 mlmcodaq-sn.f butfilter.f four1.f /usr/local/sac/lib/sac.a -o mlmcodaq-sn
c
c ng: geometrical spreading factor (2 for body, 1 for surface)
c number of data
parameter(max=65536,pai=3.1415926,ng=2)
dimension data(2*max)
dimension sactmp(max),sacdum(max),env(max)
character*30 filename,dumfilename
character*2 cmp !UD,NS,EW
real noise
complex ai
c
ai=(0.,1.) !imaginary unit
c
c write(6,*) 'Enter filename,t0,t2s,fl,fh'
read(5,*) dumfilename,t0,t2s,fl,fh,cmp
npole=4
c
cfilename=dumfilename(1:16)//'.UD.sac'
filename=dumfilename(1:16)//'.'//cmp//'.sac'
call rsac1('acc/'//filename,sacdum,ndata,beg,sampt,max,nerr)
c
c * baseline correction *
offset=0. !initialize
do j=1,ndata
offset=offset+sacdum(j)
enddo
offset=offset/ndata
do j=1,ndata
sacdum(j)=sacdum(j)-offset
enddo
c
c ** Integration
sactmp(1)=sacdum(1)*sampt
do j=1,ndata-1
sactmp(j+1)=sactmp(j)+sacdum(j)*sampt
sacdum(j)=0. !initialize
enddo
sacdum(ndata)=0. !initialize
c
c ** Filtering
call butfil(sactmp,sacdum,ndata,fl,fh,npole,sampt)
c
c ** forward transform **
isign=1
c ** input data **
do i=1,ndata
data(2*i-1)=sactmp(i)
data(2*i)=0.
enddo
do i=ndata+1,max
data(2*i-1)=0.
data(2*i)=0.
enddo
c
c ** execute fft subrutine **
call four1(data,max,isign)
c
c ** Analytic signal calculation
do i=1,max/2+1 !positive frequency
omega=2.*pai*(i-1)/(max*sampt)
data(2*i-1)=2.*data(2*i-1)
data(2*i)=2.*data(2*i)
enddo
do i=max/2+2,max !negative frequency
data(2*i-1)=0.
data(2*i)=0.
enddo
c
c ** execute inverse fft subrutine **
isign=-1
call four1(data,max,isign)
c
nshift=int(t0/sampt)
if (nshift.ge.1) then
do i=1,max-nshift
t=(i-1)*sampt
j=i+nshift
env(i)=(data(2*j-1)/max)**2+(data(2*j)/max)**2 !MS envelope
enddo
else if (nshift.le.0) then
do i=1,max+nshift
t=(i-1)*sampt
j=i-nshift
env(j)=(data(2*i-1)/max)**2+(data(2*i)/max)**2 !MS envelope
enddo
endif
c
c
c ** signal and noise estimation
nei=int(77.5/sampt)
ne=int(80./sampt)
tp=t2s/(2.*1.73)
t15s=0.75*t2s
nni=int((tp-3.5)/sampt) !2.5s window before P
nne=int((tp-1.)/sampt)
signal=0. !initialize
noise=0. !initialize
do i=nni,nne
noise=noise+env(i)*sampt
enddo
do i=nei,ne
signal=signal+env(i)*sampt
enddo
if (noise.eq.0.) then
open(11,file='env-data.dat')
write(11,*)
close(11)
open(12,file='env-mlfit.dat')
write(12,*)
close(12)
go to 98
else if (signal/noise.lt.2.) then
open(11,file='env-data.dat')
write(11,*)
close(11)
open(12,file='env-mlfit.dat')
write(12,*)
close(12)
go to 98
else if (-t0.gt.tp-3.5) then
open(11,file='env-data.dat')
write(11,*)
close(11)
open(12,file='env-mlfit.dat')
write(12,*)
close(12)
go to 98
else if (t15s.gt.75.) then
open(11,file='env-data.dat')
write(11,*)
close(11)
open(12,file='env-mlfit.dat')
write(12,*)
close(12)
go to 98
endif
c
open(11,file='env-data.dat')
cdo i=1,ndata
do i=1,max
t=(i-1)*sampt
c write(11,*) t,env(i),data(2*i-1)/max,data(2*i)/max
write(11,*) t,env(i)
enddo
close(11)
c
c ** estimation of coda Q and Nakagami's m-value **
ns=int(t15s/sampt)
c
c *codaQ
betamin=0.
betamax=1.7 !2 is too large
betaacc=1e-5
c
beta=rtbis1(betamin,betamax,betaacc,env,sampt,ns,ne,ng,max)
c
r2sum=0. !initialize
do i=ns,ne
t=(i-1)*sampt
r2sum=r2sum+env(i)*exp(beta*t)*t**ng
enddo
n=ne-ns+1
alpha=r2sum/n
c
c *Nakagami-m
ammin=0.5
ammax=5.
amacc=1e-5
c
am=rtbis2(ammin,ammax,amacc,env,sampt,ns,ne,ng,alpha,beta,max)
c
fmid=(fl+fh)*0.5
codaq=2.*pai*fmid/beta
errm=sqrt((trigamma(am)-1./am)/n) ! 1 standard deviation
cwrite(6,*) filename,beta,am,alpha,codaq
write(6,*) am,errm,codaq,filename
cwrite(6,*) am,errm,codaq,signal/noise,filename
c
open(12,file='env-mlfit.dat')
do i=ns,ne
t=(i-1)*sampt
write(12,*) t,alpha*exp(-beta*t)/t**ng
enddo
close(12)
c
98 continue
stop
end
c
function func1(beta,env,sampt,ns,ne,ng)
real func1,beta,env(*),sampt,tsum,r2sum,r2tsum
integer ns,ne,ng
tsum=0.
r2sum=0.
r2tsum=0.
n=ne-ns+1
do i=ns,ne
t=(i-1)*sampt
tsum=tsum+t
r2sum=r2sum+env(i)*exp(beta*t)*t**ng
r2tsum=r2tsum+env(i)*t*exp(beta*t)*t**ng
enddo
func1=tsum*r2sum/n-r2tsum
cwrite(6,*) 'func1',tsum,r2sum,r2tsum,func1
return
end
c
FUNCTION rtbis1(x1,x2,xacc,env,sampt,ns,ne,ng,nmax)
INTEGER JMAX,ns,ne,ng,nmax
REAL rtbis1,x1,x2,xacc,func1,sampt,env(nmax)
EXTERNAL func1
PARAMETER (JMAX=40)
INTEGER j
REAL dx,f,fmid,xmid
fmid=func1(x2,env,sampt,ns,ne,ng)
f=func1(x1,env,sampt,ns,ne,ng)
if(f*fmid.ge.0.) pause 'root must be bracketed in rtbis'
if(f.lt.0.)then
rtbis1=x1
dx=x2-x1
else
rtbis1=x2
dx=x1-x2
endif
do 11 j=1,JMAX
dx=dx*.5
xmid=rtbis1+dx
fmid=func1(xmid,env,sampt,ns,ne,ng)
if(fmid.le.0.)rtbis1=xmid
if(abs(dx).lt.xacc .or. fmid.eq.0.) return
11 continue
pause 'too many bisections in rtbis'
END
C (C) Copr. 1986-92 Numerical Recipes Software (m41.
c
function func2(am,env,sampt,ns,ne,ng,alpha,beta)
real func2,am,env(*),sampt,tsum,lr2sum,r2sum,alpha,beta
real digamma,ltsum
integer ns,ne,ng
external digamma
tsum=0.
ltsum=0.
r2sum=0.
lr2sum=0.
n=ne-ns+1
do i=ns,ne
t=(i-1)*sampt
tsum=tsum+t
ltsum=ltsum+log(t)
r2sum=r2sum+env(i)*exp(beta*t)*t**ng
lr2sum=lr2sum+log(env(i))
enddo
func2=n*(1.+log(am)-log(alpha)-digamma(am))+ng*ltsum+beta*tsum+
& lr2sum-r2sum/alpha
cwrite(6,*) 'func2',func2,ltsum,tsum,lr2sum,r2sum
return
end
c
FUNCTION rtbis2(x1,x2,xacc,env,sampt,ns,ne,ng,alpha,beta,nmax)
INTEGER JMAX,ns,ne,ng,nmax
REAL rtbis2,x1,x2,xacc,func1,sampt,env(nmax),alpha,beta
EXTERNAL func2,digamma
PARAMETER (JMAX=40)
INTEGER j
REAL dx,f,fmid,xmid
fmid=func2(x2,env,sampt,ns,ne,ng,alpha,beta)
f=func2(x1,env,sampt,ns,ne,ng,alpha,beta)
if(f*fmid.ge.0.) pause 'root must be bracketed in rtbis'
if(f.lt.0.)then
rtbis2=x1
dx=x2-x1
else
rtbis2=x2
dx=x1-x2
endif
do 11 j=1,JMAX
dx=dx*.5
xmid=rtbis2+dx
fmid=func2(xmid,env,sampt,ns,ne,ng,alpha,beta)
cwrite(6,*) j,xmid,fmid,digamma(xmid)
if(fmid.le.0.)rtbis2=xmid
if(abs(dx).lt.xacc .or. fmid.eq.0.) return
11 continue
pause 'too many bisections in rtbis'
END
C (C) Copr. 1986-92 Numerical Recipes Software (m41.
C
c
function digamma ( x0 )
c function digama ( x )
c*********************************************************************72
c
cc DIGAMA calculates DIGAMMA ( X ) = d ( LOG ( GAMMA ( X ) ) ) / dX
c
c Modified:
c
c 18 January 2008
c
c Author:
c
c Jose Bernardo
c Modifications by John Burkardt
c
c Reference:
c
c Jose Bernardo,
c Algorithm AS 103:
c Psi ( Digamma ) Function,
c Applied Statistics,
c Volume 25, Number 3, 1976, pages 315-317.
c
c Parameters:
c
c Input, double precision X, the argument of the digamma function.
c 0 < X.
c
c Output, integer IFAULT, error flag.
c 0, no error.
c 1, X <= 0.
c
c Output, double precision DIGAMA, the value of the digamma function at X.
c
implicit none
real x0,digamma
double precision c
parameter ( c = 8.5D+00 )
double precision d1
parameter ( d1 = -0.5772156649D+00 )
double precision digama
double precision r
double precision s
parameter ( s = 0.00001D+00 )
double precision s3
parameter ( s3 = 0.08333333333D+00 )
double precision s4
parameter ( s4 = 0.0083333333333D+00 )
double precision s5
parameter ( s5 = 0.003968253968D+00 )
double precision x
double precision y
c
c
x=dble(x0)
c Check the input.
c
if ( x .le. 0.0D+00 ) then
digama = 0.0D+00
digamma=sngl(digama)
return
end if
c
c Initialize.
c
y = x
digama = 0.0D+00
c
c Use approximation if argument <= S.
c
if ( y .le. s ) then
digama = d1 - 1.0D+00 / y
digamma=sngl(digama)
return
end if
c
c Reduce to DIGAMA(X + N) where (X + N) >= C.
c
10 continue
if ( y .lt. c ) then
digama = digama - 1.0D+00 / y
y = y + 1.0D+00
go to 10
end if
c
c Use Stirling's (actually de Moivre's) expansion if argument > C.
c
r = 1.0D+00 / y
digama = digama + log ( y ) - 0.5D+00 * r
r = r * r
digama = digama - r * ( s3 - r * ( s4 - r * s5 ) )
digamma=sngl(digama)
return
end
C
C
function trigamma ( x0 )
c function trigam ( x )
c*********************************************************************72
c
cc TRIGAM calculates trigamma(x) = d**2 log(gamma(x)) / dx**2
c
c Modified:
c
c 28 March 1999
c
c Author:
c
c BE Schneider
c Modifications by John Burkardt
c
c Reference:
c
c BE Schneider,
c Algorithm AS 121:
c Trigamma Function,
c Applied Statistics,
c Volume 27, Number 1, pages 97-99, 1978.
c
c Parameters:
c
c Input, double precision X, the argument of the trigamma function.
c 0 < X.
c
c Output, integer IFAULT, error flag.
c 0, no error.
c 1, X <= 0.
c
c Output, double precision TRIGAM, the value of the trigamma function at X.
c
implicit none
real x0,trigamma
double precision a
parameter ( a = 0.0001D+00 )
double precision b
parameter ( b = 5.0D+00 )
double precision b2
parameter ( b2 = 0.1666666667D+00 )
double precision b4
parameter ( b4 = -0.03333333333D+00 )
double precision b6
parameter ( b6 = 0.02380952381D+00 )
double precision b8
parameter ( b8 = -0.03333333333D+00 )
double precision trigam
double precision x
double precision y
double precision z
c
x=dble(x0)
c Check the input.
c
if ( x .le. 0.0D+00 ) then
trigam = 0.0D+00
trigamma=sngl(trigam)
return
end if
z = x
c
c Use small value approximation if X <= A.
c
if ( x .le. a ) then
trigam = 1.0D+00 / x / x
trigamma=sngl(trigam)
return
end if
c
c Increase argument to ( X + I ) >= B.
c
trigam = 0.0D+00
10 continue
if ( z .lt. b ) then
trigam = trigam + 1.0D+00 / z / z
z = z + 1.0D+00
go to 10
end if
c
c Apply asymptotic formula if argument is B or greater.
c
y = 1.0D+00 / z / z
trigam = trigam + 0.5D+00 *
& y + ( 1.0D+00
& + y * ( b2
& + y * ( b4
& + y * ( b6
& + y * b8 )))) / z
trigamma=sngl(trigam)
return
end
and here is the lscodaq-sn.f file
c f77 lscodaq-sn.f butfilter.f four1.f /usr/local/sac/lib/libsac.a /usr/local/sac/lib/libsacio.a -o lscodaq-sn
c pgf77 -tp k8-32 lscodaq-sn.f ../butfilter.f ../four1.f /usr/local/sac/lib/sac.a -o lscodaq-sn
c
c ng: geometrical spreading factor (2 for body, 1 for surface)
c number of data
parameter(max=65536,pai=3.1415926,ng=2)
dimension data(2*max)
dimension sactmp(max),sacdum(max),env(max)
character*30 filename,dumfilename
character*2 cmp !UD,NS,EW
real noise
complex ai
c
ai=(0.,1.) !imaginary unit
c
c write(6,*) 'Enter filename,t0,t2s,fl,fh'
read(5,*) dumfilename,t0,t2s,fl,fh,cmp
npole=4
fmid=(fl+fh)/2.
c
cfilename=dumfilename(1:16)//'.UD.sac'
filename=dumfilename(1:16)//'.'//cmp//'.sac'
call rsac1('acc/'//filename,sacdum,ndata,beg,sampt,max,nerr)
c
c * baseline correction *
offset=0. !initialize
do j=1,ndata
offset=offset+sacdum(j)
enddo
offset=offset/ndata
do j=1,ndata
sacdum(j)=sacdum(j)-offset
enddo
c
c ** Integration
sactmp(1)=sacdum(1)*sampt
do j=1,ndata-1
sactmp(j+1)=sactmp(j)+sacdum(j)*sampt
sacdum(j)=0. !initialize
enddo
sacdum(ndata)=0. !initialize
c
c ** Filtering
call butfil(sactmp,sacdum,ndata,fl,fh,npole,sampt)
c
c ** forward transform **
isign=1
c ** input data **
do i=1,ndata
data(2*i-1)=sactmp(i)
data(2*i)=0.
enddo
do i=ndata+1,max
data(2*i-1)=0.
data(2*i)=0.
enddo
c
c ** execute fft subrutine **
call four1(data,max,isign)
c
c ** Analytic signal calculation
do i=1,max/2+1 !positive frequency
omega=2.*pai*(i-1)/(max*sampt)
data(2*i-1)=2.*data(2*i-1)
data(2*i)=2.*data(2*i)
enddo
do i=max/2+2,max !negative frequency
data(2*i-1)=0.
data(2*i)=0.
enddo
c
c ** execute inverse fft subrutine **
isign=-1
call four1(data,max,isign)
c
nshift=int(t0/sampt)
if (nshift.ge.1) then
do i=1,max-nshift
t=(i-1)*sampt
j=i+nshift
env(i)=(data(2*j-1)/max)**2+(data(2*j)/max)**2 !MS envelope
enddo
else if (nshift.le.0) then
do i=1,max+nshift
t=(i-1)*sampt
j=i-nshift
env(j)=(data(2*i-1)/max)**2+(data(2*i)/max)**2 !MS envelope
enddo
endif
c
c
c ** signal and noise estimation
nei=int(77.5/sampt)
ne=int(80./sampt)
tp=t2s/(2.*1.73)
t15s=0.75*t2s
nni=int((tp-3.5)/sampt) !2.5s window before P
nne=int((tp-1.)/sampt)
signal=0. !initialize
noise=0. !initialize
do i=nni,nne
noise=noise+env(i)*sampt
enddo
do i=nei,ne
signal=signal+env(i)*sampt
enddo
if (noise.eq.0.) then
open(11,file='env-data.dat')
write(11,*)
close(11)
open(12,file='env-lsfit.dat')
write(12,*)
close(12)
go to 98
else if (signal/noise.lt.2.) then
open(11,file='env-data.dat')
write(11,*)
close(11)
open(12,file='env-lsfit.dat')
write(12,*)
close(12)
go to 98
else if (-t0.gt.tp-3.5) then
open(11,file='env-data.dat')
write(11,*)
close(11)
open(12,file='env-lsfit.dat')
write(12,*)
close(12)
go to 98
else if (t15s.gt.75.) then
open(11,file='env-data.dat')
write(11,*)
close(11)
open(12,file='env-lsfit.dat')
write(12,*)
close(12)
go to 98
endif
c
open(11,file='env-data.dat')
cdo i=1,ndata
do i=1,max
t=(i-1)*sampt
c write(11,*) t,env(i),data(2*i-1)/max,data(2*i)/max
write(11,*) t,env(i) !original MS envelope
env(i)=log(env(i)*t**ng) !Ms envelope with correction of t**ng
enddo
close(11)
c
c ** estimation of coda Q and Nakagami's m-value **
ns=int(t15s/sampt)
c
c *codaQ
sum=0.
xsum=0.
ysum=0.
xxsum=0.
xysum=0.
do i=ns,ne
t=(i-1)*sampt
sum=sum+1.
xsum=xsum+t
ysum=ysum+env(i)
xysum=xysum+t*env(i)
xxsum=xxsum+t**2
enddo
delta=xxsum*sum-xsum**2
b=1./delta*(-sum*xysum+xsum*ysum)
a0=1./delta*(xxsum*ysum-xsum*xysum)
codaq=2.*pai*fmid/b
c codaqinv=b/(2.*pai*fmid)
c
** standard deviation **
resid=0. !initialize
do i=ns,ne
t=(i-1)*sampt
resid=resid+(env(i)-(a0-b*t))**2
enddo
resid=resid/(ne-ns+1)
stdcodaqinv=sqrt(resid/delta*sum)/(2.*pai*fmid)
c write(6,*) codaq,b,a0
c write(6,*) codaq,resid,filename
write(6,*) codaq,filename
c
open(12,file='env-lsfit.dat')
do i=ns,ne
t=(i-1)*sampt
c write(12,*) t,exp(a0-b*t)
c write(12,*) t,a0-b*t
write(12,*) t,exp(a0-b*t)/t**ng
enddo
close(12)
c
98 continue
stop
end
Have a great time