!
!	This is the same as the wavemovie.rt in the book,
!       but (slightly) more readable
!	-- Jun ( 5-24-93 )
!
subroutine wavemovie( nz,nx,nt,nw, dz,dx, 	lambda, v, ppp)

integer			:: nz,nx,nt,nw, iz,ix,it,iw
real			:: z0,x0,dz,dx,dt,dw, lambda, v, ppp(nz,nx,nt)
real	  		:: pi2, phase,	       	       x,    w, wov
complex   		:: aa, a, b, c, cshift
temporary complex	:: cd(nx), pp(nx)

pi2 = 8.*atan(1.);	dw = v*pi2/lambda;	 x0 = nx*dx/3
			dt =   pi2/(nt*dw);	 z0 = nz*dz/3

call null( ppp, nz*nx*nt )			! zero ppp

do iw= 1, nw 					! sum each w separately
	w = iw*dw;	    wov = w/v		! freq / vel

	do ix= 1, nx 				! initial conditions for a
		x = ix*dx - x0			! downward propagating plane wave 
		phase  = -wov*(-x*sin(pi2/24)-z0*cos(pi2/24))
	pp(ix) = cexp(cmplx(0.,phase))
	end do

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

	do iz= 1, nz 				! extrapolation in depth

		do ix= 2, nx-1 			! diffraction term
			cd(ix) = aa*pp(ix+1) + (1.-2.*aa)*pp(ix) + aa*pp(ix-1)
		end do
		cd( 1) = 0.
		cd(nx) = 0.
		call ctris( nx, -a, a, b, c, -c, cd, pp)

		cshift = cexp( cmplx( 0., wov*dz))
		do ix= 1, nx 			! shifting term
			pp(ix) = pp(ix) * cshift
		end do

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