module OW_image_areal1
!
! Imaging Condition -- crosscorrelation
!

use OW_type
use OW_parms_areal
use OW_fftw_areal

implicit none

contains
!----------------------------------------------------------
subroutine owimage(ws,wr,GG) 

integer					:: ish, ihso, ixw,iw,iz
complex					:: GG(:,:),wr(:),ws(:)

do ish=1,img%sh%n
   ihso=nint(((ish-1)*img%sh%d+img%sh%o)/(img%xm%d))
   do ixw=1,rec%h%n
      if (ixw-ihso > 1 .and. ixw-ihso < rec%h%n .and. ixw+ihso > 1 .and. ixw+ihso < rec%h%n) then
          GG(ish,ixw)=GG(ish,ixw)+conjg(ws(ixw-ihso))*(wr(ixw+ihso))!*(cmplx(0,rec%w%all(iw)))
      end if
   end do
end do

end subroutine
!----------------------------------------------------------
subroutine owimage1(ws,wr,GG) 

integer					:: ish, ihso, ixw,iw,iz
complex					:: GG(:,:,:),wr(:,:),ws(:,:)

do ish=1,img%sh%n
   ihso=nint(((ish-1)*img%sh%d+img%sh%o)/(img%xm%d))
   do ixw=1,rec%h%n
      if (ixw-ihso > 1 .and. ixw-ihso < rec%h%n .and. ixw+ihso > 1 .and. ixw+ihso < rec%h%n) then
          do iz=1,img%z%n
             GG(iz,ish,ixw)=GG(iz,ish,ixw)+conjg(ws(iz,ixw-ihso))*(wr(iz,ixw+ihso))!*(cmplx(0,rec%w%all(iw)))
          end do
      end if
   end do
end do

end subroutine
!----------------------------------------------------------
subroutine owimage2(ws,wr,GG) 

integer					:: ish, ihso, ixw,iw,iz
real					:: GG(:,:,:),wr(:,:),ws(:,:)

do ish=1,img%sh%n
   ihso=nint(((ish-1)*img%sh%d+img%sh%o)/(img%xm%d))
   do ixw=1,rec%h%n
      if (ixw-ihso > 1 .and. ixw-ihso < rec%h%n .and. ixw+ihso > 1 .and. ixw+ihso < rec%h%n) then
          do iz=1,img%z%n
             GG(iz,ish,ixw)=GG(iz,ish,ixw)+ws(iz,ixw-ihso)*wr(iz,ixw+ihso)
          end do
      end if
   end do
end do

end subroutine
!       phase=cmplx(cos(pi/4.), -sin(pi/4.)) / sqrt(rec%w%all(iwa))           
!----------------------------------------------------------
subroutine owimage_stereo(ws,wr,GG) 

integer				:: ish, ihso, ixw,iw,iz,ikx,ikz,ikxr
integer				:: begx1,begx2,endx1,endx2,begz1,begz2,endz1,endz2
real				:: GG(:,:,:),wr(:,:),ws(:,:),hso,den,kxr,rs
complex,allocatable			:: cwr(:,:),cws(:,:)
complex,allocatable			:: tfldst(:,:)

allocate(tfldst(azn,axn),cwr(azn,axn),cws(azn,axn))
call owfftw_stereo(tfldst)

do ish=1,img%sh%n
   ihso=nint(((ish-1)*img%sh%d+img%sh%o)/(img%xm%d))
   hso=(ish-1)*img%sh%d+img%sh%o
   do ixw=1,rec%h%n
      if (ixw-int(axn/2.) .le. 0) then
         begx1=1
         begx2=int(axn/2.)-ixw+1
      else
         begx1=ixw-int(axn/2.)
         begx2=1
      end if
      if (ixw+int(axn/2.) .ge. img%xm%n) then
         endx1=img%xm%n
         endx2=int(axn/2.)+img%xm%n-ixw+1
      else
         endx1=ixw+int(axn/2.)
         endx2=axn
      end if
      if (ixw-ihso > 1 .and. ixw-ihso < rec%h%n .and. ixw+ihso > 1 .and. ixw+ihso < rec%h%n) then
          do iz=1,img%z%n
             cwr=0.;cws=0.
             if (iz-int(azn/2.) .le. 0) then
                begz1=1
                begz2=int(azn/2.)-iz+1
             else
                begz1=iz-int(azn/2.)
                begz2=1
             end if
             if (iz+int(azn/2.) .ge. img%z%n) then
                endz1=img%z%n
                endz2=int(azn/2.)+img%z%n-iz+1
             else
                endz1=iz+int(azn/2.)
                endz2=azn
             end if
             tfldst=0.;tfldst(begz2:endz2,begx2:endx2)=wr(begz1:endz1,begx1:endx1); 
             call sfftw_execute(planstF); cwr=tfldst*sqrt(1./(1.*azn*axn))
             tfldst=0.;tfldst(begz2:endz2,begx2:endx2)=ws(begz1:endz1,begx1:endx1);
             call sfftw_execute(planstF); cws=tfldst*sqrt(1./(1.*azn*axn))
             do ikx=1,kix%n
                do ikz=1,kiz%n
                   if (kiz%all(ikz) .ne. 0.) then
                      den1 = 1 + tan(2*(geodip-ddip)*pi/180.)*kix%all(ikx)/kiz%all(ikz)
                      den2 = 1 + tan(2*(geodip+ddip)*pi/180.)*kix%all(ikx)/kiz%all(ikz)
                      if (den1 .gt. tiny(real(1)) .and. den2 .gt. tiny(real(1))) then
                         rs1 = (tan(2*(geodip-ddip)*pi/180.) - kix%all(ikx)/kiz%all(ikz))/den
                         rs2 = (tan(2*(geodip+ddip)*pi/180.) - kix%all(ikx)/kiz%all(ikz))/den
                         kxr1 = - kiz%all(ikz) * rs1
                         kxr2 = - kiz%all(ikz) * rs2
                      end if
                      if (kxr1 .ge. 0. .and. kxr1 .le. maxval(kix%all) & 
                          kxr2 .ge. 0. .and. kxr2 .le. maxval(kix%all)) then
                         ikxr=int(kxr/kix%d) + 1.5
                         if (ikxr .ge. 1) then
                            cws(ikz,ikx)=cws(ikz,ikx)*exp(cmplx(0.,-2*kix%all(ikx)*hso))
                            cwr(ikz,ikxr)=cwr(ikz,ikxr)*exp(cmplx(0.,kxr*hso))
                            GG(iz,ish,ixw)=GG(iz,ish,ixw)+real(cws(ikz,ikx)*cwr(ikz,ikxr))
                         end if
                      else if (kxr .ge. minval(kix%all) .and. kxr .le. kix%all(kix%n-1)) then
                         ikxr=int((kxr+minval(kix%all))/kix%d)+int(kix%n/2.)+2
                         if (ikxr .ge. 1) then
                            cws(ikz,ikx)=cws(ikz,ikx)*exp(cmplx(0.,-2*kix%all(ikx)*hso))
                            cwr(ikz,ikxr)=cwr(ikz,ikxr)*exp(cmplx(0.,kxr*hso))
                            GG(iz,ish,ixw)=GG(iz,ish,ixw)+real(cws(ikz,ikx)*cwr(ikz,ikxr))
                         end if
                      end if
                   end if
                end do
             end do
          end do
      end if
   end do
end do

!                   if (kiz%all(ikz) .ne. 0. .and. &
!                      (-kix%all(ikx)/kiz%all(ikz) .le. tan((geodip+ddip)*pi/180.) .and. &
!                      (-kix%all(ikx)/kiz%all(ikz) .ge. tan((geodip-ddip)*pi/180)))) then

deallocate(tfldst,cwr,cws)

end subroutine

end module
