6-产生单频波之和的电影程序-p71.txt

地球物理大师Claebout的几个有用的程序

!公式见P70
! wave field extrapoiation program
!computer program to make a movie of a sum of monochromatic waves.(Lynn,Gonzalez,JFC,Hale)

!implicit undefined (n-z)

complex cd(48),ce(48),cf(48),q(48),aa,a,b,c,cshift
real p(96,48,12),phase,pi2,dz,dx,v,z0,x0,dt,dw,lambda,w,wov,x
integer ix,nx,iz,nz,iw,nw,it,nt
character*4   fname

nx=48; nz=96; nt=12; dx=2; dz=1; pi2=2*3.1415926
v=1; lambda=nz*dz/4; dw=v*pi2/lambda; dt=pi2/(nt*dw); nw=2

do iz=1,nz
	do ix=1,nx
		do it=1,nt
			p(iz,ix,it)=0.0 
		end do
	end do
end do


do iw=1,nw				!superimpose nw frequencies
	w=iw*dw; wov=w/v		! fequency/velocity
	x0=nx*dx/3;z0=nz*dz/3		! (x0,z0)-the source coordinate ??

	do ix=1,nx			! intial conditions for a collapsing spherical wave
		x=ix*dx-x0
		phase=-wov*sqrt(z0**2+x**2)
		q(ix)=cexp(cmplx(0.,phase)) ! the wave field value on the surface (iz=0) ??
	end do

	aa=dz/(4.*(0.,-1.)*wov*dx**2)   ! tridiagonal matrix coefficients
	a=-aa; b=1.+2.*aa; c=-aa

	do iz=1,nz			! extrapolation in depth
		do ix=2,nx-1		! cd(1,nx) restores the right terms of the equation while z==iz
			cd(ix)=aa*q(ix+1)+(1.-2.*aa)*q(ix)+aa*q(ix-1)	! the diffraction term
		end do
		cd(1)=0.; cd(nx)=0.	! 零斜率边界条件

		call ctris(nx,-a,a,b,c,-c,cd,q,ce,cf)
				! "ctris"solves complex tridiagonal equations, i.e."rtris", with complex variable
				! the returned value is the left terms of the equation while z==iz+1

		cshift=cexp(cmplx(0.,wov*dz))

		do ix=1,nx			! the shifting term
			q(ix)=q(ix)*cshift
		end do

		do it=1,nt			! evolution in time
			cshift=cexp(cmplx(0.,-w*it*dt))
			do ix=1,nx
				p(iz,ix,it)=p(iz,ix,it)+q(ix)*cshift
			end do
		end do
	end do
end do

do it=1,nt
	WRITE(fname,'(I2)')it
	open(3,file='E:\2D_Viscoelastic_W\movie'//fname//'.grd')

	WRITE(3,'(A4)') 'DSAA'
	WRITE(3,'(2I6)') nx+1,nz+1
	WRITE(3,'(2I6)') 0.0,nx
	WRITE(3,'(2I6)') 0.0,nz
	WRITE(3,'(2E12.4)') -1.,1.

	write(3,*)((p(iz,ix,it),iz=1,nz),ix=1,nx)
end do
stop
end

subroutine ctris(n,endl,a,b,c,endr,d,q,e,f)
complex q(n),d(n),f(n),e(n),a,b,c,den,endl,endr
e(1)=-a/endl; f(1)=d(1)/endl
do i=2,n-1
	den=b+c*e(i-1)
	e(i)=-a/den
	f(i)=(d(i)-c*f(i-1))/den
end do
q(n)=(d(n)-c*f(n-1))/(endr+c*e(n-1))

do i=n-1,1,-1
	q(i)=e(i)*q(i+1)+f(i)
end do
return
end
! q0=bl*q(1); qnxpl=br*q(nx)	
! cd(1)=aa*q(2)+(1.-2.*aa)*q(1)+aa*q0
! cd(nx)=aa*q(nx-1)+(1.-2.*aa)*q(nx)+aa*qnxpl
! endl=c*bl+b; endr=a*br+b
! call ctris(nx,endl,a,b,c,endr,cd,q,ce,cf)	

! 上面代码是Dave Hale编制的程序,替换上面程序中的：

! cd(1)=0.; cd(nx)=0.	
! call ctris(nx,-a,a,b,c,-c,cd,q,ce,cf)
! 可得到快速而节省内存的零边界条件
! 注意：对于零值边界应bl=br=0;对于零斜率边界则应bl=br=1.令bl,br为复数，则得到4.4节导出的吸收边界条件