module OW_image_areal
!
! Imaging Condition -- crosscorrelation
!

use OW_type
use OW_parms_areal

implicit none

contains
!=====================================================================================
subroutine owimage(adj,ww,ws,wr,image,ssou)

logical					:: adj
logical, optional                                 :: ssou
integer					:: ihx, ihxo, ixw,ihy, ihyo, iyw,iw,sign
complex					:: wr(:,:),ws(:,:),phase
real					:: image(:,:,:,:),ww

phase = ww**2 !cmplx(cos(pi/4.), -sin(pi/4.)) * sqrt(ww)

      if (adj) then
         do ihy=1,img%shy%n
            ihyo=nint(((ihy-1)*img%shy%d+img%shy%o)/(img%ym%d))
            do iyw=1,kym%n
               do ihx=1,img%shx%n
                  ihxo=nint(((ihx-1)*img%shx%d+img%shx%o)/(img%xm%d))
                  do ixw=1,kxm%n
                     if (ixw-ihxo >= 1 .and. ixw-ihxo <= kxm%n .and. ixw+ihxo >= 1 .and. ixw+ihxo <= kxm%n .and. &
                         iyw-ihyo >= 1 .and. iyw-ihyo <= kym%n .and. iyw+ihyo >= 1 .and. iyw+ihyo <= kym%n ) then
                        if (ssou) then
                           ws(ixw-ihxo,iyw-ihyo)=ws(ixw-ihxo,iyw-ihyo)+wr(ixw+ihxo,iyw+ihyo)*image(ixw,ihx,iyw,ihy)*phase
                        else
                           wr(ixw+ihxo,iyw+ihyo)=wr(ixw+ihxo,iyw+ihyo)+ws(ixw-ihxo,iyw-ihyo)*image(ixw,ihx,iyw,ihy)*phase
                        end if
                     end if
                  end do
               end do
            end do
         end do
      else
         do ihy=1,img%shy%n
            ihyo=nint(((ihy-1)*img%shy%d+img%shy%o)/(img%ym%d))
            do iyw=1,kym%n
               do ihx=1,img%shx%n
                  ihxo=nint(((ihx-1)*img%shx%d+img%shx%o)/(img%xm%d))
                  do ixw=1,kxm%n
                     if (ixw-ihxo >= 1 .and. ixw-ihxo <= kxm%n .and. ixw+ihxo >= 1 .and. ixw+ihxo <= kxm%n .and. &
                         iyw-ihyo >= 1 .and. iyw-ihyo <= kym%n .and. iyw+ihyo >= 1 .and. iyw+ihyo <= kym%n ) then
                         image(ixw,ihx,iyw,ihy)=image(ixw,ihx,iyw,ihy)+conjg(ws(ixw-ihxo,iyw-ihyo))*wr(ixw+ihxo,iyw+ihyo)*phase
                     end if
                  end do
               end do
            end do
         end do
      end if
end subroutine
!=============================================
 real function rms1(w)
   real              :: w(:)
   rms1=sqrt(sum(w*w)/size(w))
   if (rms1.lt.tiny(1.)) rms1=1.
 end function
!=============================================
 real function rms2(w)
   real              :: w(:,:)
   integer           :: i
   rms2=0.
   do i=1,size(w,2)
      rms2=rms2+sqrt(sum(w(:,i)*w(:,i))/(size(w,1)*size(w,2)))
   end do
   if (rms2.lt.tiny(1.)) rms2=1.
 end function
!=============================================
end module
