module OW_pspi_plane

use OW_type
use OW_parms_areal
use OW_wei_plane

implicit none

!public						:: owpspi_init, owpspi
private						:: ni,ixs,nwb,node,vsal,raypar

integer						:: ni,ixs,nwb,node
real						:: wb,raypar

contains

!---------------------------------------------------------
subroutine owpspi_plane_init(ni_in,ixs_in,wb_in,node_in,raypar_in)

integer						:: ni_in,ixs_in,node_in
real						:: wb_in,raypar_in

ni=ni_in
ixs=ixs_in
node=node_in
wb=wb_in
wb=floor(wb/img%z%d)*img%z%d
nwb=floor(wb/img%z%d)+1
raypar=raypar_in
end subroutine
!---------------------------------------------------------
subroutine owpspi_plane(taper,recv,source,vsal,vel,G1)

integer					:: iw,inode,ino,i1,nout,iowrk,ih,ixphase
real					:: vel(:,:),vsal,taper(:),owrk,xmed,xphase
complex					:: recv(:,:,:),source(:,:,:),G1(:,:,:)
complex, allocatable			:: GG(:,:,:),synshot(:),synsou(:),rwr(:,:),rws(:,:)
complex, allocatable			:: synshot1(:),synsou1(:)
integer, external				:: omp_get_thread_num

if(exist_file("wr")) allocate(rwr(img%z%n,nxcut))
if(exist_file("ws")) allocate(rws(img%z%n,nxcut))

if (raypar < 0) then 
   xmed = img%xm%all(img%xm%n)
else if (raypar > 0) then 
   xmed = img%xm%all(1)
else 
   xmed=0
end if

xmed=(img%xm%all(1)+img%xm%all(img%xm%n))/2.

! $OMP parallel default(shared) private(inode) ! num_threads(8)

       inode=omp_get_thread_num()+1;

! $OMP do private(iw)
do iw = 1,rec%w%n
!  write(0,*) iw , rec%w%all(iw)/(2*pi) , abs(1/(4*raypar*img%xm%d))
   if (rec%w%all(iw)/(2*pi) <= abs(1/(4*raypar*img%xm%d))) then

   if (allocated(synsou)) deallocate(synsou); allocate(synsou(img%xm%n));    synsou=0
   if (allocated(synshot)) deallocate(synshot); allocate(synshot(img%xm%n)); synshot=0
   synsou=source(1,iw,1)*exp(cmplx(0,rec%w%all(iw)*raypar*(img%xm%all-xmed)))
   do ixs=1,rec%x%n
      owrk=(ixs-1)*rec%x%d+rec%x%o+rec%h%o
      iowrk=int((owrk-img%xm%all(1))/rec%h%d)+1.5
      ixphase=int((rec%x%all(ixs)-img%xm%all(1))/img%xm%d)+1.5
      xphase=img%xm%all(ixphase)-xmed
      do ih=1,rec%h%n
         synshot(iowrk+ih-1)=synshot(iowrk+ih-1)+recv(ih,iw,ixs)*exp(cmplx(0,rec%w%all(iw)*raypar*xphase))*csqrt(cmplx(0,rec%w%all(iw)))
      end do
   end do

   if (allocated(synsou1)) deallocate(synsou1); allocate(synsou1(nxcut));    synsou1=0
   if (allocated(synshot1)) deallocate(synshot1); allocate(synshot1(nxcut)); synshot1=0
   
   synsou1=synsou(ixcut:nxcut+ixcut-1)
   synshot1=synshot(ixcut:nxcut+ixcut-1)

   if (exist_file("planeshot")) call srite("planeshot",synsou1,8*nxcut)

      if(allocated(GG)) deallocate(GG); allocate(GG(img%z%n,img%sh%n,nxcut)); GG=0 !,node))
      call owwei_init(ni,iw,wb,nwb,node,inode)
      call owwei(taper,synshot1,synsou1,vsal,vel,GG,rwr,rws)! (:,:,:,inode),rwr)

      if(exist_file("wr")) call srite("wr",rwr,8*img%z%n*nxcut)
      if(exist_file("ws")) call srite("ws",rws,8*img%z%n*nxcut)

! $OMP critical
!    G1(:,:,:)=G1(:,:,:)+GG(:,:,:,inode)
      G1=G1+GG
! $OMP end critical   
      if (mod(iw,10)==1) write(0,5) " FREQUENCY # ",iw," OUT OF ",rec%w%n," FROM NODE: ",inode
   end if
end do
! $OMP End do
! $OMP end parallel

5 format(A13,i3,A8,i3,A12,i3)

end subroutine 
end module
