[離散化波数法プログラムの主要部] 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