module OW_dotst

  use OW_type
  use OW_exec_time
  use OW_parms_areal
  use OW_taper_areal
  use OW_wei_areal

  implicit none

  contains

  subroutine owdotst(taper,recv,source,vref,vel)

    integer				:: i,j,k
    real					:: taper(:)
    real					:: vel(:,:),vref(:,:,:)
    complex				:: recv(:,:),source(:,:)
    real, allocatable			:: deltaI1(:,:,:),deltaS1(:,:)
    real, allocatable			:: deltaI2(:,:,:),deltaS2(:,:)
    real, dimension (2)		 	:: dot1, dot2
    real, allocatable, dimension(:) 		:: mod1, mod2
    real, allocatable, dimension(:) 		:: dat1, dat2

    allocate(mod1(img%z%n*img%xm%n))
    allocate(mod2(img%z%n*img%xm%n))
    allocate(dat1(img%z%n*img%sh%n*img%xm%n))
    allocate(dat2(img%z%n*img%sh%n*img%xm%n))
    allocate(deltaS1(img%z%n,img%xm%n))
    allocate(deltaS2(img%z%n,img%xm%n))
    allocate(deltaI1(img%z%n,img%sh%n,img%xm%n))
    allocate(deltaI2(img%z%n,img%sh%n,img%xm%n))

    call random_number(deltaS1)
    call random_number(deltaI2)

    deltaS2=0.;deltaI1=0.
    deltaI2=deltaI2*1000
    deltaS1=deltaS1*1000


    do i=1,img%xm%n
       do j=1,img%sh%n
          do k=1,img%z%n
             dat1(((i-1)*img%sh%n+(j-1))*img%z%n+k)=deltaI1(k,j,i)
             dat2(((i-1)*img%sh%n+(j-1))*img%z%n+k)=deltaI2(k,j,i)
          end do
       end do
    end do

    do j=1,img%xm%n
       do k=1,img%z%n
          mod1((j-1)*img%z%n+k)=deltaS1(k,j)
          mod2((j-1)*img%z%n+k)=deltaS2(k,j)
       end do
    end do

    write(0,*) sum(abs(dat1)),sum(abs(dat2)),sum(abs(mod1)),sum(abs(mod2))

    call owwemva(.false.,taper,recv,source,vref,vel,deltaS1,deltaI1)
    call owwemva(.true., taper,recv,source,vref,vel,deltaS2,deltaI2)

    do i=1,img%xm%n
       do j=1,img%sh%n
          do k=1,img%z%n
             dat1(((i-1)*img%sh%n+(j-1))*img%z%n+k)=deltaI1(k,j,i)
             dat2(((i-1)*img%sh%n+(j-1))*img%z%n+k)=deltaI2(k,j,i)
          end do
       end do
    end do

    do j=1,img%xm%n
       do k=1,img%z%n
          mod1((j-1)*img%z%n+k)=deltaS1(k,j)
          mod2((j-1)*img%z%n+k)=deltaS2(k,j)
       end do
    end do

    write(0,*) sum(abs(dat1)),sum(abs(dat2)),sum(abs(mod1)),sum(abs(mod2))

    dot1( 1) = dot_product( dat1, dat2)
    dot1( 2) = dot_product( mod1, mod2)
    write (0,*) dot1

    dot2( 1) = dot_product( dat1, dat2)
    dot2( 2) = dot_product( mod1, mod2)
    write (0,*) dot2

  end subroutine
end module
