<div class="moz-text-flowed" style="font-family: -moz-fixed">［離散化波数法プログラムの主要部］
parameter ( maxxdv=256, maxkdv=50, maxtim=1024, maxsit=20 )
dimension bnd(maxxdv),rnx(maxxdv),rnz(maxxdv),vs(2),ro(2)
dimension ip(4*maxkdv+2),x(maxsit)
dimension wavei(maxtim),waveo(maxtim,maxsit)
complex a(4*maxkdv+2,maxtim/2),w(4*maxkdv*2)
complex g(4*maxkdv+2,4*maxkdv+2),h(4*maxkdv+2)
complex gw(maxxdv,2,2),hw(maxxdv,2)
complex cwavei(maxtim),cwaveo(maxtim)
c
pai=2.0*acos(0.0)
c
ntim=1024
dt=0.05
f1=0.1
f2=1.0
nkdv=20
c
vs(1)=900.0e2
ro(1)=1.9
vs(2)=2400.0e2
ro(2)=2.3
dmax=1000.0e2
ang=0.0*pai/180.0
c
nxdv=256
dx=100.0e2
do 1 i=1,nxdv
bnd(i)=0.0
1 continue
i1=int(0.25*real(nxdv)+0.5)
i2=int(0.75*real(nxdv)+0.5)
do 2 i=i1,i2
bnd(i)=dmax*0.5*(1.0-cos(2.0*pai*real(i-i1)/real(i2-i1)))
2 continue
c
nsit=20
xl=real(nxdv)*dx
do 3 i=1,nsit
x(i)=xl*real(i-1)/real(nsit-1)
3 continue
c
t0=10.0
fp=0.25
amp=1.0
tp=1.0/fp
call ricker(maxtim,ntim,dt,wavei,t0,tp,amp)
c
call dw2d2l(maxxdv,bnd,rnx,rnz,g,h,gw,hw,nxdv,dx,vs,ro,
& maxtim,maxkdv,a,ntim,dt,nkdv,f1,f2,ang,w,ip)
c
call dwsurf(maxxdv,maxtim,maxkdv,a,wavei,cwavei,nxdv,dx,
& ntim,dt,nkdv,f1,f2,maxsit,x,waveo,cwaveo,nsit,vs,ang)
c
WRITE(*,'(10f7.1)') (waveo(i,5),i=1,ntim)
c
call plots(0.0,0.0,99)
call plot(5.0,17.5,-3)
call plot(0.0,1.2,-3)
call wveplt(maxtim,ntim,wavei,0.0125,0.2)
call plot(0.0,-1.2,-3)
do 4 j=1,nsit
call wveplt(maxtim,ntim,waveo(1,j),0.0125,0.2)
call plot(0.0,-0.8,-3)
4 continue
xwid=real(nsit-1)*0.8
call plot(0.0,0.8*real(nsit),-3)
iplt=3
do 5 i=1,nxdv
xx=-1.0-bnd(i)*0.00001
yy=-1.0*xwid*real(i)/real(nxdv-1)
call plot(xx,yy,iplt)
iplt=2
5 continue
call plot(-1.0,0.0,3)
call plot(-1.0,-1.0*xwid,2)
call plot(0.0,0.0,999)
c
stop
end
［離散化波数法計算サブルーチン部］
c *********************************************************************
subroutine dw2d2l(maxxdv,bnd,rnx,rnz,g,h,gw,hw,nxdv,dx,vs,ro,
& maxtim,maxkdv,a,ntim,dt,nkdv,f1,f2,ang,w,ip)
c *********************************************************************
dimension bnd(maxxdv),rnx(maxxdv),rnz(maxxdv)
dimension vs(2),ro(2),ip(4*maxkdv+2)
complex a(4*maxkdv+2,maxtim/2),w(4*maxkdv+2)
complex g(4*maxkdv+2,4*maxkdv+2),h(4*maxkdv+2)
complex gw(maxxdv,2,2),hw(maxxdv,2)
complex cnu1,cnu2,ccosx,csinx,cexpx,cexpg,en
c
pai=2.0*acos(0.0)
tim=real(ntim)*dt
nf1=int(f1*tim)+1
nf2=int(f2*tim)+1
xl=real(nxdv)*dx
dk=2.0*pai/xl
nk21=2*nkdv+1
c
call bndset(maxxdv,bnd,rnx,rnz,nxdv,dx)
c
do 1 nf=nf1,nf2
c
frq=real(nf-1)/tim
WRITE(*,*) 'Freq(Hz)=',frq
omg=2.0*pai*frq
c
rk1=(omg/vs(1))
rk2=(omg/vs(2))
rk0=rk2*sin(ang)
rnu0=rk2*cos(ang)
rmu1=ro(1)*vs(1)**2.0
rmu2=ro(2)*vs(2)**2.0
c
do 2 n=1,nk21
c
rk=rk0+2.0*pai*real(n-(nkdv+1))/xl
cnu1=csqrt(cmplx(rk1**2.0-rk**2.0,0.0))
cnu2=csqrt(cmplx(rk2**2.0-rk**2.0,0.0))
c
do 3 m=1,nxdv
x=(real(m-1)+0.5)*dx
en=cexp(cmplx(0.0,(rk-rk0)*x))
ccosx=ccos(cnu1*cmplx(bnd(m),0.0))
csinx=csin(cnu1*cmplx(bnd(m),0.0))
cexpx=cexp(cnu2*cmplx(0.0,bnd(m)))
gw(m,1,1)=cmplx(2.0*dk,0.0)*ccosx*en
gw(m,1,2)=cmplx(-1.0*dk,0.0)*cexpx*en
gw(m,2,1)=cmplx(2.0*rmu1*dk,0.0)
& *(cmplx(0.0,rk*rnx(m))*ccosx
& -cnu1*cmplx(rnz(m),0.0)*csinx)*en
gw(m,2,2)=cmplx(-1.0*rmu2*dk,0.0)
& *(cmplx(0.0,rk*rnx(m))
& +cmplx(0.0,rnz(m))*cnu2)*cexpx*en
if(n.eq.1) then
cexpg=cexp(cmplx(0.0,-1.0*rnu0*bnd(m)))
hw(m,1)=cexpg
hw(m,2)=cmplx(rmu2,0.0)
& *(cmplx(0.0,rk0*rnx(m))
& -cmplx(0.0,rnu0*rnz(m)))*cexpg
end if
3 continue
c
call fft(maxxdv,nxdv,gw(1,1,1),-1)
call fft(maxxdv,nxdv,gw(1,1,2),-1)
call fft(maxxdv,nxdv,gw(1,2,1),-1)
call fft(maxxdv,nxdv,gw(1,2,2),-1)
if(n.eq.1) then
call fft(maxxdv,nxdv,hw(1,1),-1)
call fft(maxxdv,nxdv,hw(1,2),-1)
end if
c
do 4 m=nxdv-nkdv+1,nxdv
g(n,m-(nxdv-nkdv))=conjg(gw(m,1,1))
g(nk21+n,m-(nxdv-nkdv))=conjg(gw(m,1,2))
g(n,nk21+m-(nxdv-nkdv))=conjg(gw(m,2,1))
g(nk21+n,nk21+m-(nxdv-nkdv))=conjg(gw(m,2,2))
if(n.eq.1) then
h(m-(nxdv-nkdv))=conjg(hw(m,1))
h(nk21+m-(nxdv-nkdv))=conjg(hw(m,2))
end if
4 continue
do 5 m=1,nkdv+1
g(n,nkdv+m)=conjg(gw(m,1,1))
g(nk21+n,nkdv+m)=conjg(gw(m,1,2))
g(n,nk21+nkdv+m)=conjg(gw(m,2,1))
g(nk21+n,nk21+nkdv+m)=conjg(gw(m,2,2))
if(n.eq.1) then
h(nkdv+m)=conjg(hw(m,1))
h(nk21+nkdv+m)=conjg(hw(m,2))
end if
5 continue
c
2 continue
c
call cdcomp((4*maxkdv+2),(2*nk21),g,w,ip)
call csolve((4*maxkdv+2),(2*nk21),g,h,a(1,nf),ip)
c
1 continue
c
return
end
c *********************************************************************
subroutine dwsurf(maxxdv,maxtim,maxkdv,a,wavei,cwavei,nxdv,dx,
& ntim,dt,nkdv,f1,f2,maxsit,x,waveo,cwaveo,nsit,vs,ang)
c *********************************************************************
dimension x(maxsit),vs(2)
dimension wavei(maxtim),waveo(maxtim,maxsit)
complex a(4*maxkdv+2,maxtim/2),cwavei(maxtim),cwaveo(maxtim)
c
pai=2.0*acos(0.0)
tim=real(ntim)*dt
nf1=int(f1*tim)+1
nf2=int(f2*tim)+1
xl=real(nxdv)*dx
dk=2.0*pai/xl
nk21=2*nkdv+1
c
do 1 i=1,ntim
cwavei(i)=cmplx(wavei(i),0.0)
1 continue
call fft(maxtim,ntim,cwavei,-1)
c
do 2 j=1,nsit
c
do 3 i=1,ntim
cwaveo(i)=(0.0,0.0)
3 continue
c
do 4 nf=nf1,nf2
c
frq=real(nf-1)/tim
omg=2.0*pai*frq
rk0=(omg/vs(2))*sin(ang)
do 5 n=1,nk21
rk=rk0+2.0*pai*real(n-(nkdv+1))/xl
cwaveo(nf)=cwaveo(nf)
& +cmplx(2.0*dk,0.0)*a(n,nf)
& *cexp(cmplx(0.0,rk*x(j)))
5 continue
cwaveo(nf)=cwaveo(nf)*cwavei(nf)
cwaveo(ntim+2-nf)=conjg(cwaveo(nf))
c
4 continue
c
call fft(maxtim,ntim,cwaveo,1)
do 6 i=1,ntim
waveo(i,j)=real(cwaveo(i))
6 continue
c
2 continue
c
return
end
c *********************************************************************
subroutine bndset(maxxdv,bnd,rnx,rnz,nxdv,dx)
c *********************************************************************
dimension bnd(maxxdv),rnx(maxxdv),rnz(maxxdv)
c
do 1 i=2,nxdv-1
dgdx=(bnd(i+1)-bnd(i-1))/(2.0*dx)
sqrtgx=sqrt(1.0+dgdx**2.0)
rnx(i)=-1.0*dgdx/sqrtgx
rnz(i)=1.0/sqrtgx
1 continue
rnx(1)=rnx(2)
rnz(1)=rnz(2)
rnx(nxdv)=rnx(nxdv-1)
rnz(nxdv)=rnz(nxdv-1)
c
return
end
c *********************************************************************
subroutine fft(maxtim,ntim,c,ifft)
c *********************************************************************
complex c(maxtim)
c
call fast(ntim,c,maxtim,ifft)
if(ifft.eq.-1) then
do 1 i=1,ntim
c(i)=c(i)/cmplx(real(ntim),0.0)
1 continue
end if
c
return
end
</div>