7-(w-x-z)域内的偏移程序-p76.txt

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

# 公式见P75
# Migration in the (omega,x,z) domain

real q(48,64),pi2,alpha,dt,dtau,dw
complex cq(48,64),cd(48),ce(48),cf(48),aa,a,b,c,cshift
integer ix,nx,iz,nz,iw,nw,it,nt

open(4,file='plot36',status='new',access='direct',form='unformatted',recl=1)

nt=64; nz=nt; nx=48; pi2=2*3.1415926
dt=1.; dtau=1.; dw=pi2/(dt*nt); nw=nt/2	!dtau=dz/v
alpha=.25;			# alpha=v*v*dtau/(4*dx*dz)

do iz=1,nz
	do ix=1,nx
	{	q(ix,iz)=0.; cq(ix,iz)=0.	}

do it=nt/3,nt,nt/4
	do ix=1,4		# Broadened impluse source
	{	cq(ix,it)=(5.-ix); cq(ix,it+1)=(5.-ix)	}

call rowcc(nx,nt,cq,+1.,+1.)	# F.T.over time

do iz=1,nz			# iz and iw loops interchangeble
{
	do iw=2,nw
	{
		aa=-alpha/((0.,-1.)*(iw-1)*dw)
		a=-aa; b=1.+2.*aa; c=-aa
		do ix=z,nx-1
			cd(ix)=aa*cq(ix+1,iw)+(1.-2.*aa)*cq(ix,iw)+aa*cq(ix-1,iw)
		cd(1)=0.; cd(nx)=0.
		call ctris(nx,-a,a,b,c,-c,cd,cq(1,iw),ce,cf)
		cshift=cexp(cmplx(0.,-(iw-1)*dw*dtau))
		do ix=1,nx
			cq(ix,iw)=cq(ix,iw)*cshift
		do ix=1,nx
			q(ix,iz)=q(ix,iz)+cq(ix,iw)	# q(t=0)=Q(w)求和
	}
}
write(4,rec=1)((q(ix,iz),iz=1,nz),ix=1,nx)
stop
end