subroutine phasemig( u, nt, nx, dt, dx, image, ntau, dtau, v)
integer	nt,nw,nx,ntau, iw,ikx,itau
real	dt,dx, w,w0,dw, kx,kx0,dkx,dtau, kzkz, v, sig1,sig2,pi, signum
complex	u(nt,nx), image(ntau,nx), c
			pi = 3.14159265;	sig1 =+1.;  sig2 =-1.
call ft1axis( 0, sig1, nt, nx, u)   
call ft2axis( 0, sig2, nt, nx, u)

w0  = -pi/dt;	dw = 2.*pi/(nt*dt);	nw = nt
kx0 = -pi/dx;	dkx= 2.*pi/(nx*dx)

call zero( ntau*nx, image)
do itau = 1, ntau {
do ikx  = 2, nx   {        kx = kx0 + (ikx-1) * dkx
do iw = 2, nw     {        w  = w0  + (iw-1)  * dw
	    kzkz = w*w - v*v * kx*kx
	    if (kzkz >0) { c = cexp(cmplx(0., -signum(w)*dtau*sqrt(kzkz))) }
	    else         { c = 0.                                          }
	    u(iw,ikx) = u(iw,ikx) * c
	    image(itau,ikx) = image(itau,ikx) + u(iw,ikx)
	    }}}
call ft2axis( 1, sig2, ntau, nx, image)
return; end