-
-
194.05 KB, 下载次数: 1
parameter (nn=500000)
dimension x1(nn),x2(nn),x3(nn)
dimension xns(nn),xew(nn),xud(nn)
dimension f(nn),f1(nn),f2(nn),f3(nn),period(nn)
character title*32,dim*3,idsn*40,odsn*40,idsn_ctl*40
c
write(*,*) 'input control file name (*.ctl)'
read(*,*) idsn_ctl
c
open(10,file=idsn_ctl)
read(10,*)
read(10,*) idsn
read(10,*)
read(10,*) sf1,nd
read(10,*)
read(10,*) band
read(10,*)
read(10,*) odsn
c
c
call wvread_csv(x1,x2,x3,n,nd,title,hz,dim,idsn)
c
c
nd2=nd/2+1
c
dt=1.0/hz
t=float(nd)*dt
c
ls=int(sf1*hz)
do 110 l=1,nd
ll=l+ls
x1(l)=x1(ll)
x2(l)=x2(ll)
x3(l)=x3(ll)
110 continue
c
call fpac(nd,x1,nd,dt,100,f,g,g,nd2,nfold,df)
do l=1,nfold
f1(l)=f(l)
enddo
call swin(nfold,f1,g,nd2,100,df,band)
do m=1,nd2
xns(m)=f1(m)
enddo
c
call fpac(nd,x2,nd,dt,100,f,g,g,nd2,nfold,df)
do l=1,nfold
f2(l)=f(l)
enddo
call swin(nfold,f2,g,nd2,100,df,band)
do m=1,nd2
xew(m)=f2(m)
enddo
c
call fpac(nd,x3,nd,dt,100,f,g,g,nd2,nfold,df)
do l=1,nfold
f3(l)=f(l)
enddo
call swin(nfold,f3,g,nd2,100,df,band)
do m=1,nd2
xud(m)=f3(m)
enddo
c
c
do k=1,nd2
period(k)=t/k
enddo
c
open(70,file=odsn)
write(70,200)
200 format(2x,'Period(s)',3x,'Amp_NS',6x,'Amp_EW',6x,'Amp_UD')
do j=1,nd2
write(70,220) period(j),xns(j),xew(j),xud(j)
220 format(7e12.4)
enddo
c
close(70)
stop
end
c---------------------------------------------------------------------
c subroutine for fourier,power spectra and autocorrelation
c---------------------------------------------------------------------
c
c coded by y.ohasaki
c
subroutine fpac(n,x,nd1,dt,ind,f,g,r,nd2,nfold,df)
c
complex c(8192)
dimension x(nd1),f(nd2),g(nd2),r(nd2)
c
c initialization
do 110 m=1,n
c(m)=cmplx(x(m),0.)
110 continue
nt=2
120 if(nt.ge.n) go to 130
nt=nt*2
go to 120
130 if(nt.eq.n) go to 150
do 140 m=n+1,nt
c(m)=(0.,0.)
140 continue
150 nfold=nt/2+1
t=real(nt)*dt
df=1./t
c
c
c fourier transform
c
call fast(nt,c,8192,-1)
c
c fourier spectrum
c
if(ind.eq.001) goto 180
do 151 k=1,nt
c(k)=c(k)/nt
151 continue
do 160 k=1,nfold
f(k)=cabs(c(k))*dt
160 continue
if(ind.eq.100) return
c
c power spectrum
c
if(ind.eq.101) go to 180
g(1)=f(1)**2/t
do 170 k=2,nfold-1
g(k)=2.*f(k)**2/t
170 continue
g(nfold)=f(nfold)**2/t
if(mod(ind,10).eq.0) return
c
c autocorrelation
c
180 do 190 k=1,nt
c(k)=c(k)*conjg(c(k))
190 continue
call fast(nt,c,8192,+1)
r0=real(c(1))
do 200 j=1,nfold
r(j)=real(c(j))/r0
200 continue
c
return
end
c
c---------------------------------------------------------------------
c subroutine for fast fourier transform
c---------------------------------------------------------------------
subroutine fast(n,x,nd1,ind)
complex x(nd1),temp,theta
j=1
do 140 i=1,n
if(i.ge.j) go to 110
temp=x(j)
x(j)=x(i)
x(i)=temp
110 m=n/2
120 if(j.le.m) go to 130
j=j-m
m=m/2
if(m.ge.2) go to 120
130 j=j+m
140 continue
kmax=1
150 if(kmax.ge.n) return
istep=kmax*2
do 170 k=1,kmax
theta=cmplx(0.0,3.141593*float(ind*(k-1))/float(kmax))
do 160 i=k,n,istep
j=i+kmax
temp=x(j)*cexp(theta)
x(j)=x(i)-temp
x(i)=x(i)+temp
160 continue
170 continue
kmax=istep
go to 150
end
c
c
c---------------------------------------------------------------------
c subroutine for smoothed spectra by parzen's spectral window
c---------------------------------------------------------------------
c coded by y.ohsaki
c
subroutine swin(nfold,f,g,nd,ind,df,band)
c
dimension f(nd),g(nd),w(101),g1(4497),g2(4497)
c
c initialization
c
t=1./df
if(band.eq.0.) go to 120
udf=1.854305/band*df
if(udf.gt.0.5) go to 230
lmax=ifix(2./udf)+1
if(lmax.gt.101) go to 240
c
c spectral window
c
w(1)=0.75*udf
do 110 l=2,lmax
dif=1.570796*real(l-1)*udf
w(l)=w(1)*(sin(dif)/dif)**4
110 continue
c
c conversion from fourier to power spectrum
c
120 if(ind.ne.100) go to 140
g(1)=f(1)**2/t
do 130 k=2,nfold-1
g(k)=2.*f(k)**2/t
130 continue
g(nfold)=f(nfold)**2/t
c
c smoothing of power spectrum
c
140 if(band.eq.0.) go to 210
ll=lmax*2-1
ln=ll-1+nfold
lt=(ll-1)*2+nfold
le=lt-lmax+1
do 150 k=1,lt
g1(k)=0.
150 continue
do 160 k=1,nfold
g1(ll-1+k)=g(k)
160 continue
do 180 k=lmax,le
s=w(1)*g1(k)
do 170 l=2,lmax
s=s+w(l)*(g1(k-l+1)+g1(k+l-1))
170 continue
g2(k)=s
180 continue
do 190 l=2,lmax
g2(ll+l-1)=g2(ll+l-1)+g2(ll-l+1)
g2(ln-l+1)=g2(ln-l+1)+g2(ln+l-1)
190 continue
do 200 k=1,nfold
g(k)=g2(ll-1+k)
200 continue
c
c smoothed fourier spectrum
c
210 f(1)=sqrt(g(1)*t)
do 220 k=2,nfold-1
f(k)=sqrt(g(k)*t/2.)
220 continue
f(nfold)=sqrt(g(nfold)*t)
return
c
230 write(30,601)
stop
240 write(30,602)
stop
c
c format statements
c
601 format('bandwidth is too narrow')
602 format('bandwidth is too wide')
end
c---- 攇宍僨乕僞傪撉傒崬傓僒僽儖乕僠儞 ----------------------------
subroutine wvread_csv(x1,x2,x3,n,nd1,title,hz,dim,idsn)
dimension x1(nd1),x2(nd1),x3(nd1)
character title*32,dim*3,idsn*40,stcode*6,date*22,dim0*19
c
open(60,file=idsn)
read(60,601) stcode
601 format(11x,a6)
read(60,602) n
602 format(16x,i5)
read(60,603) hz
603 format(15x,f5.1)
read(60,*) dim0
if(dim0.eq.'UNIT = gal(cm/s/s)') dim='ACC'
if(dim0.eq.'UNIT = kine(cm/s) ') dim='VEL'
if(dim0.eq.'UNIT = cm ') dim='DIS'
read(60,604) date
604 format(15x,a22)
title=stcode(1:6)//' , '//date(1:22)
read(60,*)
do i=1,n
read(60,*) temp,x1(i),x2(i),x3(i)
enddo
close(60)
c
return
end
c
194.05 KB, 下载次数: 1
vvt 发表于 2019-6-5 14:48
你的这个代码应该是要求 NUMBER OF DATA= 14900 这里是 2的整幂次方的。
而且必须要求和 ctl 文件里的一致 ...
要哭了| 欢迎光临 Fortran Coder (http://bbs.fcode.cn/) | Powered by Discuz! X3.2 |