! Phase-shift migration and diffraction
!
module phaseshift_mod

  use adj_mod

  implicit none
  real, private		 			:: o1,d1,o2,d2,o3,d3,amax
  real,private,dimension(:,:),allocatable	:: vel
  integer, private				:: n1,n2,n3

  contains

  subroutine phsh_init(vel,n1,d1,o1,n2,d2,o2)!,n3,d3,o3)
    real,dimension(n1_in,n2_in)	:: vel_in
    real 			:: o1_in,o2_in,o3_in,d1_in,d2_in,d3_in
    integer			:: n1_in,n2_in,n3_in

    n1 = n1_in;    n2 = n2_in;  !  n3 = n3_in
    o1 = o1_in;    o2 = o2_in;  !  o3 = o3_in
    d1 = d1_in;    d2 = d2_in;  !  d3 = d3_in
    amax=amax_in
    allocate(vel(n1,n2))
    vel=vel_in
  end subroutine

  function phsh_lop ( adj, add, xx, yy) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    real,dimension(:,:)  :: xx,yy 
!    call adjnull (adj,add,xx,yy )
    call phsh_lop2(adj,add,xx,yy )
    stat=0
  end function phsh_lop

  subroutine phsh_lop2(adj, add, xx, yy)
   integer 		:: ix,iz,it,ih, xstart, xend
   real			:: w,t,z,
   complex		:: xx(:,:), yy(:,:)
   real			:: cosin,a,cosin15,tp
   real 		:: pi
   logical		:: add,adj

  pi=4*atan(1.) 
!  call adjnull( adj, add, xx,n1*n2, yy,n1*n2)



do iw= 1, nw 					! sum each w separately
	w = iw*dw;     				! 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
	do iz= nz,1,-1 				! upward continuation
		wov = w/v(iz)
		sqa = wov**2 - k%all(ik)**2
		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


      if (amp .ne. 0.) then
        if ( adj ) then
            xx(iz,ix )=xx(iz,ix )+yy(it,ix+ih)*amp
        else
            yy(it,ix+ih)=yy(it,ix+ih)+xx(iz,ix )*amp
        end if
      end if
    end do
  end do
  end do
  end subroutine phsh_lop2
end module phaseshift_mod





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

return;	end
