module OW_image_areal
!
! Imaging Condition -- crosscorrelation
!

use OW_type
use OW_parms_areal
!use OW_fftw_areal
!use OW_exec_time

implicit none

contains
!=====================================================================================
subroutine owimage(adj,iw,ws,wr,image,vel,ssou)
real                                              :: vel(:)
logical					:: adj
logical, optional                                 :: ssou
integer					:: ish, ihso, ith, ixw,iw,sign
complex					:: wr(:),ws(:),phase
real					:: image(:,:),hto

phase = cmplx(cos(pi/4.), -sin(pi/4.)) * sqrt(rec%w%all(iw))

if (timeshift) then
   if (adj) then
      do ixw=1,kxm%n
         do ith=1,img%th%n
            hto=(ith-1)*img%th%d+img%th%o
!            if (2*hto/vel(ixw) >= 1.and. 2*hto/vel(ixw) <= img%z%n) then
               if (ssou) then
                  ws(ixw)=ws(ixw)+wr(ixw)*exp(cmplx(0,-2*hto*rec%w%all(iw)))*image(ixw,ith)!*phase
               else
                  wr(ixw)=wr(ixw)+ws(ixw)*exp(cmplx(0,-2*hto*rec%w%all(iw)))*image(ixw,ith)!*phase
               end if
!            end if
         end do
      end do
   else
      do ixw=1,kxm%n
         do ith=1,img%th%n
            hto=(ith-1)*img%th%d+img%th%o
!            if (2*hto/vel(ixw) >= 1.and. 2*hto/vel(ixw) <= img%z%n) then
               image(ixw,ith)=image(ixw,ith)+wr(ixw)*conjg(ws(ixw))*exp(cmplx(0,2*hto*rec%w%all(iw)))!*phase
!            end if
         end do
      end do
   end if
else
   if (adj) then
      do ixw=1,kxm%n
         do ish=1,img%sh%n
            ihso=nint(((ish-1)*img%sh%d+img%sh%o)/(img%xm%d))
            if (ixw-ihso >= 1 .and. ixw-ihso <= kxm%n .and. ixw+ihso >= 1 .and. ixw+ihso <= kxm%n) then
               if (ssou) then
                  ws(ixw-ihso)=ws(ixw-ihso)+wr(ixw+ihso)*image(ixw,ish)!*phase!conjg(phase)!
               else
                  wr(ixw+ihso)=wr(ixw+ihso)+ws(ixw-ihso)*image(ixw,ish)!*phase!conjg(phase)!
               end if
            end if
         end do
      end do
   else
      do ixw=1,kxm%n
         do ish=1,img%sh%n
            ihso=nint(((ish-1)*img%sh%d+img%sh%o)/(img%xm%d))
            if (ixw-ihso >= 1 .and. ixw-ihso <= kxm%n .and. ixw+ihso >= 1 .and. ixw+ihso <= kxm%n) then
               image(ixw,ish)=image(ixw,ish)+conjg(ws(ixw-ihso))*wr(ixw+ihso)!*phase
            end if
         end do
      end do
   end if
end if
end subroutine
!=====================================================================================
end module
