module iqtransf_

  implicit none

  integer, private :: nmx, nax, nz
  integer, private :: xwin,     zwin
  logical, private :: verb

  real, dimension(:,:),   allocatable :: tmp
  real, dimension(:,:),   allocatable :: tmp1,tmp2
contains
  !----------------------------------------------------------------
  subroutine iqinit(nmxin,naxin,nzin,xwinin,zwinin,verbin)
    integer, intent(in) ::  nmxin, naxin,  nzin
    integer, intent(in) :: xwinin,       zwinin
    logical, intent(in) :: verbin

    nmx=nmxin
    nax=naxin
    nz = nzin

    xwin = xwinin
    zwin = zwinin

    verb=verbin

    allocate(tmp(nz,nmx))
    allocate(tmp1(nz,nmx),tmp2(nz,nmx))
  end subroutine iqinit
  !----------------------------------------------------------------
  !! Semblance = SUM(a)^2/SUM(a^2)
  subroutine se(adcig,iqmap)
    real :: adcig(:,:,:)
    real :: iqmap(:,:)

    integer :: iz,imx,iax
    integer :: jz,jmx,jax
    real    :: sum1, sum2, stmp, srat
    integer :: kont

    if(verb) write(0,*) "compute temp"
    do imx  =1,nmx
       do iz=1, nz
          
          sum1 = 0.
          sum2 = 0.
          kont = 0
          do jax=1,nax
             sum1 = sum1 + adcig(iz,jax,imx)
             sum2 = sum2 + adcig(iz,jax,imx)**2
             kont = kont + 1
          end do
          tmp1(iz,imx) = sum1**2 / kont
          tmp2(iz,imx) = sum2

       end do
    end do

    if(verb) write(0,*) "compute final"
    do imx  =1+xwin,nmx-xwin
       do iz=1+zwin, nz-zwin
          iqmap(iz,imx) = 0.

          sum1=0
          sum2=0
          do jmx=imx-xwin,imx+xwin
             do jz=iz-zwin,iz+zwin
                sum1=sum1+tmp1(jz,jmx)
                sum2=sum2+tmp2(jz,jmx)
             end do
          end do
          iqmap(iz,imx) = sum1/sum2 / (2*zwin+1) / (2*xwin+1)
       end do
    end do
  end subroutine se
  !----------------------------------------------------------------
  !! Differential Semblance = SQRT(SUM(da^2)/SUM(a^2))
  subroutine ds(adcig,iqmap)
    real :: adcig(:,:,:)
    real :: iqmap(:,:)

    integer :: iz,imx,iax
    integer :: jz,jmx,jax
    real    :: sum1, sum2, stmp, srat

    do imx=1+xwin,nmx-xwin
       do iz=1+zwin,nz-zwin
          iqmap(iz,imx) = 0.

          do jmx=imx-xwin,imx+xwin
             do jz=iz-zwin,iz+zwin

                sum1 = 0
                sum2 = 0
                srat = 0
                do jax=2,nax-1
                   sum1 = sum1 + ( adcig(jz,jax+1,jmx) - adcig(jz,jax-1,jmx))**2
                   sum2 = sum2 +   adcig(jz,jax  ,jmx)**2
                end do
                if(sum2 /=0) srat = sum1 / sum2

                iqmap(iz,imx) = & 
                iqmap(iz,imx) + sqrt(srat) / (2*zwin+1) / (2*xwin+1)
             end do
          end do

       end do
    end do
  end subroutine ds
  !----------------------------------------------------------------
  !! Stack Power = SUM(a^2)
  subroutine sp(adcig,iqmap)
    real :: adcig(:,:,:)
    real :: iqmap(:,:)

    integer :: iz,imx,iax
    integer :: jz,jmx,jax
    real    :: sum1, sum2, stmp, srat

    do imx=1+xwin,nmx-xwin
       do iz=1+zwin,nz-zwin

          sum1=0
          do jmx=imx-xwin,imx+xwin
             do jz=iz-zwin,iz+zwin
                do jax=1,nax
                   sum1=sum1+adcig(jz,jax,jmx)**2
                end do
             end do
          end do
          iqmap(iz,imx) = sum1
       end do
    end do
  end subroutine sp
  !----------------------------------------------------------------
  !! Stack intensity = SUM(|a|)
  subroutine si(adcig,iqmap)
    real :: adcig(:,:,:)
    real :: iqmap(:,:)

    integer :: iz,imx,iax
    integer :: jz,jmx,jax
    real    :: sum1, sum2, stmp, srat

    do imx=1+xwin,nmx-xwin
       do iz=1+zwin,nz-zwin

          sum1=0



          do jmx=imx-xwin,imx+xwin
             do jz=iz-zwin,iz+zwin
                do jax=1,nax
                   sum1=sum1+abs(adcig(jz,jax,jmx))

                end do
             end do
          end do
          iqmap(iz,imx) = sum1**2



       end do
    end do
  end subroutine si
  !----------------------------------------------------------------
end module iqtransf_









