module OW_image_areal
!
! Imaging Condition -- crosscorrelation
!

use OW_type
use OW_parms_areal
use OW_fftw_areal
use OW_exec_time
use OW_slantstk

implicit none

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

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

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)+real(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					:: wr(:,:),ws(:,:)
real					:: GG(:,:,:)

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)+real(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_tmsh(iw,ws,wr,GG) 

integer					:: ith,ixw,iw,iz
real					:: GG(:,:), hto
complex					:: wr(:),ws(:)

do ith=1,img%th%n
   hto=(ith-1)*img%th%d+img%th%o
   do ixw=1,rec%h%n
      GG(ith,ixw)=GG(ith,ixw)+real(conjg(ws(ixw)*exp(cmplx(0,-hto*rec%w%all(iw))))*wr(ixw))
   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, finish, start
complex,allocatable			:: cwr(:,:),cws(:,:)
complex,allocatable			:: tfldst(:,:)

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

do ixw=1,rec%h%n
   call owexec_time(start)
   if (ixw .le. int(axn/2.)) then
      begx1=1
      begx2=int(axn/2.)-ixw+2
   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
   do iz=1,img%z%n
       cwr=0.;cws=0.
       if (iz .le. int(azn/2.)) then
          begz1=1
          begz2=int(azn/2.)-iz+2
       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 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	
          if (ixw-ihso > 1 .and. ixw-ihso < rec%h%n .and. ixw+ihso > 1 .and. ixw+ihso < rec%h%n) then
             do ikx=1,kix%n
                do ikz=1,kiz%n
                   if (kiz%all(ikz) .ne. 0.) then
                      den = 1 + tan(2*geodip*pi/180.)*kix%all(ikx)/kiz%all(ikz)
                      if (den .gt. tiny(real(1))) then
                         rs = (tan(2*geodip*pi/180.) - kix%all(ikx)/kiz%all(ikz))/den
                         kxr = - kiz%all(ikz) * rs
                      end if
                      if (kxr .ge. 0. .and. kxr .le. maxval(kix%all)) then
                         ikxr=int(kxr/kix%d) + 1.5
                         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))
                      else if (kxr .ge. minval(kix%all) .and. kxr .le. kix%all(kix%n)) then
                         ikxr=int((minval(kix%all)-kxr)/kix%d)+int(kix%n/2.)+1.5
                         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 do
             end do	
          end if
       end do
   end do	
   call owexec_time(finish)
   if (mod(ixw,10).eq.1) write(0,*) "ixw = ",ixw, " in ", 10*(finish-start), " s"
end do

deallocate(tfldst,cwr,cws)

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

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

allocate(cwr(azn,axn),cws(azn,axn),cwrs(azn,ndip),cwss(azn,ndip))

call owexec_time(start)
do ixw=1,rec%h%n
   if (ixw .le. int(axn/2.)) then
      begx1=1
      begx2=int(axn/2.)-ixw+2
   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
   do iz=1,img%z%n
       cwr=0.;cws=0.;cwrs=0.;cwss=0.;
       if (iz .le. int(azn/2.)) then
          begz1=1
          begz2=int(azn/2.)-iz+2
       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

       if(ixw .eq. 1 .or. ixw .eq. 101) write(2,222) begz1,endz1,begx1,endx1,begz2,endz2,begx2,endx2
222    format(8I5)
!       cwr(begz2:endz2,begx2:endx2)=wr(begz1:endz1,begx1:endx1); 
!       cws(begz2:endz2,begx2:endx2)=ws(begz1:endz1,begx1:endx1);
       cwr=wr(begz1:endz1,begx1:endx1); 
       cws=ws(begz1:endz1,begx1:endx1);
       call owslant(cwr,cwrs)
       call owslant(cws,cwss)

       if (ixw .eq. int(rec%h%n/2.)) then
          call srite("cwss",cwss,4*azn*ndip)
       end if

       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	
          if (ixw-ihso > 1 .and. ixw-ihso < rec%h%n .and. ixw+ihso > 1 .and. ixw+ihso < rec%h%n) then
             do ikx=1,ndip
                do ikz=1,azn
                   den = 1 + tan(2*geodip*pi/180.)*dips(ikx)
                   if (den .gt. tiny(real(1))) then
                      rs = (tan(2*geodip*pi/180.) - dips(ikx))/den
                      if (abs(rs).le.amax) then
                         irs=int((rs+amax)/ddip) + 1.5
!                         GG(iz,ish,ixw)=GG(iz,ish,ixw)+conjg(cwss(ikz,ikx))*cwrs(ikz,irs)
                         GG(iz,ish,ixw)=GG(iz,ish,ixw)+cwss(ikz,ikx)*cwrs(ikz,irs)
                      end if
                   end if
                end do
             end do	
          end if
       end do
   end do	
end do
call owexec_time(finish)
write(0,*) "ixw = ",ixw, " in ", (finish-start), " s"

deallocate(cwr,cws,cwrs,cwss)

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

integer				:: ish, ihso, ixw,iw1,iw2,iz
real				:: GG(:,:,:)
complex				:: 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
             do iw1=1,wnw
                do iw2=1,wnw
                   GG(iz,ish,ixw)=GG(iz,ish,ixw)+real(conjg(ws(iz,ixw,iw2))*wr(iz,ixw,iw1))
                end do
             end do
          end do
      end if
   end do
end do

end subroutine

end module
