# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module image_mod
  !! imaging condition
  use down_mod
  use R_mod

  !---------------
  implicit none

  complex, dimension(:,:,:,:,:), allocatable,private:: dsliceR,dsliceM,dsliceN
  real,dimension(:,:),allocatable,save,private::smoother
  integer,dimension(:,:),allocatable,save,private::jxd,jyd,jxu,jyu
  integer,dimension(:),allocatable,save,private::jhx
  integer,private,save::cy,cx
contains
!------------------------------------------------------------------------------------------
  subroutine SPimg(izs,iws,ith,wfld_u,wfld_d,X,oper)
    type(imag)::X
    integer, intent(in) :: izs,iws!,ith
    integer::i1,i2,i3,i4,i5
    complex, dimension(:,:,:), allocatable :: wfld_u,wfld_d
    integer :: st,iw,ith



    interface !!-------------------------
       integer function oper(wfld_u,wfld_d,iws,ith) result(st)
         complex, dimension(:,:,:), allocatable :: wfld_u,wfld_d
         integer,intent(in):: iws,ith
       end function oper
    end interface !!---------------------



    call zeroimage()
    !! image all offsets for all frequencies
    ith=1

!$OMP PARALLEL PRIVATE(ith)
       ith=1+omp_get_thread_num()
!$OMP DO PRIVATE(iw)

    do iw=1,down%nws
       st = oper(wfld_u,wfld_d,iw,ith)
    end do

!$OMP END DO
!$OMP END PARALLEL

    call addimagethrds(izs,X)




  end subroutine SPimg
!----------------------------------------------------------------
!! offset domain Shot-profile imaging
  integer function Ig_init() result(st)
    integer::i1,i2,i3,c,d,er
    real :: a




    !! imaging wavefield, size of imaging space (x,y,hx,hy)
    if(down%image) then
       allocate(dsliceR(down%amx%n,down%amy%n,down%ahx%n,down%ahy%n,run%o_nth),stat=er)
       if(er/=0) call seperr('cant allocate dslice array (image.F90)')
       call imagenull(dsliceR)
    end if
    if(down%mult) then
       allocate(dsliceM(down%amx%n,down%amy%n,down%ahx%n,down%ahy%n,run%o_nth),stat=er)
       if(er/=0) call seperr('cant allocate dslice array (image.F90)')
       call imagenull(dsliceM)
    end if
    if(down%normd) then
       allocate(dsliceN(down%amx%n,down%amy%n,down%ahx%n,down%ahy%n,run%o_nth),stat=er)
       if(er/=0) call seperr('cant allocate dslice array (image.F90)')
       call imagenull(dsliceN)
    end if

    a=8*(size(dsliceR)+size(dsliceM)+size(dsliceN))
    down%mem=down%mem+a




    !! ADDED BY GA
    !! precompute 1-D arrays of indicies necessary to update the
    !! wavefield at the imaging step (subroutines imaging, imaging_decon,
    !! and normfiled).
    !! Hopefully this will
    !! cut the time it takes to do imaging when computing subsurface
    !! offsets with 3D data
    call compute_wfld_indices()

    if(.not.down%smoothme) return
    c=down%rect1*2
    if(mod(c,2)==0)c=c+1
    down%rect1=c/2
    cx=c/2+1

    if(down%amy%n.gt.1) then
       d=down%rect2*2
       if(mod(d,2)==0)d=d+1
       down%rect2=d/2
       cy=d/2+1
    else
       d=1;down%rect2=1;cy=1
    end if

    allocate(smoother(c,d))
    do i2=1,d
       do i1=1,c
          smoother(i1,i2)=sqrt(real(cx-i1)**2 + real(cy-i2)**2)+1.
       end do
    end do
    smoother=1./smoother

    a=4*size(smoother)
# 130 "<stdin>"
    st=0
end function Ig_init

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !! This subroutine added by GA.
    !! Extracted from imaging
    subroutine compute_wfld_indices()

    !! local variables
    integer :: ihy,ihx,imx,imy,iyd,iyu,ixd,ixu,jhy
    real :: hy,hx,a


    allocate(jyd(down%amy%n,down%ahy%n),jyu(down%amy%n,down%ahy%n))
    allocate(jxd(down%amx%n,down%ahx%n),jxu(down%amx%n,down%ahx%n))
    allocate(jhx(down%ahx%n))

    a =4*(2*(size(jxd)+size(jyd)))





    do ihy=1, down%ahy%n
       hy = down%ahy%o + (ihy-1)*down%ahy%d
       jhy = nint(hy/down%amy%d)

       do ihx=1, down%ahx%n
          hx = down%ahx%o + (ihx-1)*down%ahx%d
          jhx(ihx) = nint(hx/down%amx%d)

          do imy=1, down%amy%n
             iyd=imy - jhy
             iyu=imy + jhy
             if ( iyd.lt.1 ) iyd=1
             if ( iyd.gt.down%amy%n ) iyd=down%amy%n
             if ( iyu.lt.1 ) iyu=1
             if ( iyu.gt.down%amy%n ) iyu=down%amy%n
             jyd(imy,ihy) = iyd
             jyu(imy,ihy) = iyu

             do imx=1, down%amx%n
                ixd=imx - jhx(ihx)
                ixu=imx + jhx(ihx)
                if ( ixd.lt.1 ) ixd=1
                if ( ixd.gt.down%amx%n ) ixd=down%amx%n
                if ( ixu.lt.1 ) ixu=1
                if ( ixu.gt.down%amx%n ) ixu=down%amx%n
                jxd(imx,ihx) = ixd
                jxu(imx,ihx) = ixu
            end do
          end do
       end do
    end do

write(0,*)'jxu=',size(jxu,1),size(jxu,2)
end subroutine compute_wfld_indices

!----------------------------------------------------------------
!! Conventional offset generating imaging: sum_w(UD*)
!! wfld_d conjugated when initialized
  subroutine imaging(wfld_u, wfld_d,iws,ith)
    complex, dimension(:,:,:), allocatable:: wfld_u, wfld_d
    integer, intent(in)::iws,ith
    integer :: ihx, ihy, imx, imy,iw
    integer :: jjhx, jjhy, jjxd, jjxu, jjyd, jjyu,nn
    !!real :: hx, hy





    nn=down%amx%n

    do ihy=1, down%ahy%n
       !!hy = down%ahy%o + (ihy-1)*down%ahy%d
       !!jhy = nint(hy/down%amy%d)

       do ihx=1, down%ahx%n
          !!hx = down%ahx%o + (ihx-1)*down%ahx%d
          !!jhx = nint(hx/down%amx%d)
          jjhx=jhx(ihx)

          do imy=down%ay%b, down%ay%n+down%ay%b-1
             !!jjyd = imy - jjhy
             !!jjyu = imy + jjhy
             jjyd=jyd(imy,ihy)
             jjyu=jyu(imy,ihy)

             do imx=down%ax%b, down%ax%n+down%ax%b-1
                jjxd = imx - jjhx
                jjxu = imx + jjhx
                !jjxd = jxd(imx,ihx)
                !jjxu = jxu(imx,ihx)
                dsliceR(imx-down%ax%b+1, imy-down%ay%b+1, ihx, ihy,ith) = &
                  & dsliceR(imx-down%ax%b+1, imy-down%ay%b+1, ihx, ihy,ith) + &
                  & wfld_d(jjxd, jjyd, iws) * wfld_u(jjxu, jjyu, iws)
            end do
          end do
       end do
    end do



  end subroutine imaging
!----------------------------------------------------------------
!! Conventional offset generating imaging: sum_w(UD*)
  integer function offIgop(wfld_u,wfld_d,iws,ith) result(st)
    complex, dimension(:,:,:), allocatable :: wfld_u,wfld_d
    integer,intent(in)::iws,ith




    st=1

    if(.not.down%dcnv) then
       call imaging(wfld_u,wfld_d,iws,ith)
    else
       call imaging_decon(wfld_u,wfld_d,iws,ith)
    end if




    st=0
  end function offIgop
!---------------------------------------------------------------
 subroutine imaging_decon(wfld_u,wfld_d,iws,ith)
   complex, dimension(:,:,:) , allocatable:: wfld_u, wfld_d
   integer, intent(in)::iws,ith
   integer :: ihx, ihy, imx, imy
   integer :: jhx, jhy, jxd, jxu, jyd, jyu
   real :: hx, hy,denom
   integer :: wxe=1,wxb=1,wye=1,wyb=1,sxb=1,sxe=1,sye=1,syb=1




   do ihy=1, down%ahy%n
      hy = down%ahy%o + (ihy-1)*down%ahy%d
      jhy = nint(hy/down%amy%d)

      do ihx=1, down%ahx%n
         hx = down%ahx%o + (ihx-1)*down%ahx%d
         jhx = nint(hx/down%amx%d)

         do imy=down%ay%b, down%ay%n+down%ay%b-1
            jyd = imy - jhy
            jyu = imy + jhy
            if ((jyd.lt.1) .or. (jyd.gt.size(wfld_d,2)) .or. &
            (jyu.lt.1) .or. (jyu.gt.size(wfld_u,2)))cycle

            do imx=down%ax%b, down%ax%n+down%ax%b-1
               jxd = imx - jhx
               jxu = imx + jhx
               if ((jxd.lt.1) .or. (jxd.gt.size(wfld_d,1)) .or. &
               (jxu.lt.1) .or. (jxu.gt.size(wfld_u,1)))cycle

               if(down%smoothme) then
                  wxb=imx-down%rect1
                  sxb=1
                  if(wxb.lt.1) sxb=cx-imx+1
                  wxb=max(1,wxb)

                  wxe=imx+down%rect1
                  sxe=cx+down%rect1
                  if(wxe.gt.down%amx%n) sxe=size(smoother,1)-(wxe-down%amx%n+1)
                  wxe=min(wxe,down%amx%n)

                  if(down%amy%n.gt.1)then
                     wyb=imy-down%rect2
                     syb=1
                     if(wyb.lt.1) syb=cy-imy+1
                     wyb=max(1,wyb)

                     wye=imy+down%rect2
                     sye=cy+down%rect2
                     if(wye.gt.down%amy%n) sye=size(smoother,2)-(wye-down%amy%n+1)
                     wye=min(wxe,down%amx%n)
                  end if
                  denom=sum(cabs(wfld_d(wxb:wxe,wyb:wye,iws))*smoother(sxb:sxe,syb:sye))
               else
                  denom=cabs(wfld_d(imx,imy,iws))
               end if
               if(denom.le.down%eps)denom=1.

               dsliceR(imx-down%ax%b+1, imy-down%ay%b+1, ihx, ihy,ith) = &
               dsliceR(imx-down%ax%b+1, imy-down%ay%b+1, ihx, ihy,ith) &
               + wfld_d(jxd, jyd, iws) &
               * wfld_u(jxu, jyu, iws)/denom
            end do
          end do
       end do
    end do



  end subroutine imaging_decon
!----------------------------------------------------------------
!!Surface Related Multiple Prediction: sum_w(UU)
  integer function multIGop(wfld_u,wfld_d,iws,ith) result(st)
    complex, dimension(:,:,:), allocatable:: wfld_u,wfld_d
    integer,intent(in)::iws,ith




    st=1
    if(down%dcnv)then
       call srmp_decon(wfld_u,iws,ith)
    else
       call srmp(wfld_u,iws,ith)
    end if



    st=0
  end function multIGop
!----------------------------------------------------------------
!! SRMP offset generating imaging: sum_w(UD*)
!! wfld_d conjugated when initialized
  subroutine srmp(wfld,iws,ith)
    complex, dimension(:,:,:), allocatable:: wfld
    integer, intent(in)::iws,ith
    integer :: ihx, ihy, imx, imy
    integer :: jhx, jhy, jxd, jxu, jyd, jyu
    real :: hx, hy




    do ihy=1, down%ahy%n
       hy = down%ahy%o + (ihy-1)*down%ahy%d
       jhy = nint(hy/down%amy%d)

       do ihx=1, down%ahx%n
          hx = down%ahx%o + (ihx-1)*down%ahx%d
          jhx = nint(hx / down%amx%d)

          do imy=down%ay%b, down%ay%n+down%ay%b-1
             jyd = imy - jhy
             jyu = imy + jhy
             if ( (jyd.lt.1) .or. (jyd.gt.size(wfld,2)) .or. &
             (jyu.lt.1) .or. (jyu.gt.size(wfld,2)) )cycle

             do imx=down%ax%b, down%ax%n+down%ax%b-1
                jxd = imx - jhx
                jxu = imx + jhx
                if ( (jxd.lt.1) .or. (jxd.gt.size(wfld,1)) .or. &
                (jxu.lt.1) .or. (jxu.gt.size(wfld,1)) )cycle

                dsliceM (imx-down%ax%b+1, imy-down%ay%b+1, ihx, ihy,ith) &
                = dsliceM(imx-down%ax%b+1, imy-down%ay%b+1, ihx, ihy,ith) &
                   + wfld(jxd, jyd, iws) &
                   * wfld(jxu, jyu, iws)
            end do
          end do
       end do
    end do



  end subroutine srmp
!---------------------------------------------------------------
 subroutine srmp_decon(wfld_u,iws,ith)
   complex, dimension(:,:,:), allocatable :: wfld_u
   integer, intent(in)::iws,ith
   integer :: ihx, ihy, imx, imy
   integer :: jhx, jhy, jxd, jxu, jyd, jyu
   real :: hx, hy,denom
   integer :: wxe=1,wxb=1,wye=1,wyb=1,sxb=1,sxe=1,sye=1,syb=1




   do ihy=1, down%ahy%n
      hy = down%ahy%o + (ihy-1)*down%ahy%d
      jhy = nint(hy/down%amy%d)

      do ihx=1, down%ahx%n
         hx = down%ahx%o + (ihx-1)*down%ahx%d
         jhx = nint(hx/down%amx%d)

         do imy=down%ay%b, down%ay%n+down%ay%b-1
            jyd = imy - jhy
            jyu = imy + jhy
            if ((jyd.lt.1) .or. (jyd.gt.size(wfld_u,2)) .or. &
            (jyu.lt.1) .or. (jyu.gt.size(wfld_u,2)))cycle

            do imx=down%ax%b, down%ax%n+down%ax%b-1
               jxd = imx - jhx
               jxu = imx + jhx
               if ((jxd.lt.1) .or. (jxd.gt.size(wfld_u,1)) .or. &
               (jxu.lt.1) .or. (jxu.gt.size(wfld_u,1)))cycle

               if(down%smoothme) then
                  wxb=imx-down%rect1
                  sxb=1
                  if(wxb.lt.1) sxb=cx-imx+1
                  wxb=max(1,wxb)

                  wxe=imx+down%rect1
                  sxe=cx+down%rect1
                  if(wxe.gt.down%amx%n) sxe=size(smoother,1)-(wxe-down%amx%n+1)
                  wxe=min(wxe,down%amx%n)

                  if(down%amy%n.gt.1)then
                     wyb=imy-down%rect2
                     syb=1
                     if(wyb.lt.1) syb=cy-imy+1
                     wyb=max(1,wyb)

                     wye=imy+down%rect2
                     sye=cy+down%rect2
                     if(wye.gt.down%amy%n) sye=size(smoother,2)-(wye-down%amy%n+1)
                     wye=min(wxe,down%amx%n)
                  end if
                  denom=sum(cabs(wfld_u(wxb:wxe,wyb:wye,iws))*smoother(sxb:sxe,syb:sye))
               else
                  denom=cabs(wfld_u(imx,imy,iws))
               end if
               if(denom.le.down%eps)denom=1.

               dsliceM(imx-down%ax%b+1, imy-down%ay%b+1, ihx, ihy,ith) = &
               dsliceM(imx-down%ax%b+1, imy-down%ay%b+1, ihx, ihy,ith) &
               + wfld_u(jxd, jyd, iws) &
               * wfld_u(jxu, jyu, iws)/denom
            end do
          end do
       end do
    end do



  end subroutine srmp_decon
!----------------------------------------------------------------
!! image normalization field: sum_w(DD*)
  integer function normIGop(wfld_u,wfld_d,iws,ith) result(st)
    complex, dimension(:,:,:), allocatable :: wfld_u,wfld_d
    integer,intent(in)::iws,ith




    st=1
    call normfield(wfld_d,iws,ith)



    st=0
  end function normIGop
!----------------------------------------------------------------
subroutine normfield(wfld,iws,ith)
    complex, dimension(:,:,:), allocatable:: wfld
    integer, intent(in)::iws,ith
    integer :: ihx, ihy, imx, imy
    integer :: jhx, jhy, jxd, jxu, jyd, jyu, i, j
    real :: hx, hy




    do ihy=1, down%ahy%n
       hy = down%ahy%o + (ihy-1)*down%ahy%d
       jhy = nint(hy/down%amy%d)

       do ihx=1, down%ahx%n
          hx = down%ahx%o + (ihx-1)*down%ahx%d
          jhx = nint(hx/down%amx%d)

          do imy=down%ay%b, down%ay%n+down%ay%b-1
             jyd = imy - jhy
             jyu = imy + jhy
             if ((jyd.lt.1) .or. (jyd.gt.size(wfld,2)) .or. &
             (jyu.lt.1) .or. (jyu.gt.size(wfld,2)))cycle

             do imx=down%ax%b, down%ax%n+down%ax%b-1
                jxd = imx - jhx
                jxu = imx + jhx
                if (jxd .lt.1 .or. jxd.gt.size(wfld,1) .or. &
                jxu .lt.1 .or. jxu.gt.size(wfld,1))cycle

                dsliceN(imx-down%ax%b+1,imy-down%ay%b+1 , ihx, ihy,ith) = &
                dsliceN(imx-down%ax%b+1, imy-down%ay%b+1, ihx, ihy,ith) &
                + wfld(jxd, jyd, iws) &
                * conjg(wfld(jxu, jyu, iws))
            end do
          end do
       end do
    end do



  end subroutine normfield
!----------------------------------------------------------------
  subroutine imagespace()
    real::maxi,maxw,sx,sy
    integer::nx,ny
    down%ax%o=down%amx%o;down%ay%o=down%amy%o
    nx=down%ax%n
    ny=down%ay%n
    if(down%amx%n.gt.1) then
       maxi=down%aix%o +(down%aix%n-1)*down%aix%d
       maxw=down%amx%o +(down%amx%n-1)*down%amx%d
       if(maxi.lt.maxw) then
          down%ax%n=1.5+(maxi-down%amx%o)/down%ax%d
       end if
       if(down%amx%o.lt.down%aix%o ) then
          down%ax%o=down%aix%o
          down%ax%b=1.5+(down%ax%o-down%amx%o)/down%ax%d
          down%ax%n=1+down%ax%n-down%ax%b
       end if
    end if
    if(down%amy%n.gt.1) then
       maxi=down%aiy%o +(down%aiy%n-1)*down%aiy%d
       maxw=down%amy%o +(down%amy%n-1)*down%amy%d
       if(maxi.lt.maxw) then
          down%ay%n=1.5+(maxi-down%amy%o)/down%ay%d
       end if
       if(down%amy%o.lt.down%aiy%o ) then
          down%ay%o=down%aiy%o
          down%ay%b=1.5+(down%ay%o-down%amy%o)/down%ay%d
          down%ay%n=1+down%ay%n-down%ay%b
       end if
    end if

! if((down%ax%n.ne.nx).or.(down%ay%n.ne.ny)) call imagearrays()
    if(down%ax%n.ne.nx) call imagearrays()!!if 2axis changes, no big deal
  end subroutine imagespace
!----------------------------------------------------------------
subroutine zeroimage()
  if(down%image) call imagenull(dsliceR)
  if(down%mult) call imagenull(dsliceM)
  if(down%normd) call imagenull(dsliceN)

end subroutine zeroimage
!----------------------------------------------------------------
subroutine imagenull(slice)
complex,dimension(:,:,:,:,:), allocatable::slice
integer:: i1,i2,i3,i4,i5




do i5=1,run%o_nth
   do i4=1,down%ahy%n
      do i3=1,down%ahx%n
         do i2=1,down%ay%n
            do i1=1,down%ax%n
               slice(i1,i2,i3,i4,i5)=cmplx(0.,0.)
            end do
         end do
      end do
   end do
end do



end subroutine imagenull
!----------------------------------------------------------------
subroutine addimagethrds(izs,X)
  type(imag)::X
  integer,intent(in)::izs
  if(trim(X%tag).eq.'R') then
!!$OMP CRITICAL
     call imagesum(dsliceR,X%bin,izs)
!!$OMP END CRITICAL
  end if
  if(trim(X%tag).eq.'M') then
!!$OMP CRITICAL
     call imagesum(dsliceM,X%bin,izs)
!!$OMP END CRITICAL
  end if
  if(trim(X%tag).eq.'N_d')then
!!$OMP CRITICAL
     call imagesum(dsliceN,X%bin,izs)
!!$OMP END CRITICAL
  end if
end subroutine addimagethrds
!----------------------------------------------------------------
subroutine imagesum(dslice,final,izs)
  complex, dimension(:,:,:,:,:), allocatable :: final,dslice
  integer,intent(in)::izs
  integer:: i1,i2,i3,i4,i5
  !! sum up frequency blocks from temp model into actual model
  !! multiply by 2 since only migrating causal frequencies?

  do i5=1,run%o_nth
     do i4=1,down%ahy%n
        do i3=1,down%ahx%n
           do i2=1,down%ay%n
              do i1=1,down%ax%n
                 final(i1,i2,i3,i4,izs)=final(i1,i2,i3,i4,izs)+dslice(i1,i2,i3,i4,i5)!*2.
            end do
         end do
      end do
   end do
end do
end subroutine imagesum
!----------------------------------------------------------------


end module image_mod
