module oc_dso_mod
  !!=======================
  !! Paul Sava    (C) 2000
  !! paul@sep.stanford.edu
  !! oclib v2.0
  !!=======================
  use sep
  use oc_file_mod
  use oc_adjnull_mod
  use oc_weight_mod
  implicit none

  integer,                       private :: ocmaxmem
  integer,                       private :: nzs,nzb, nmx,nmy,nhx,nhy,nz
  complex, dimension(:,:,:,:,:), private :: image
  real,    dimension(:),         private :: hx,hy
  real,    dimension(:,:),       private :: hh

  character(len=128), private :: w_
  logical,            private :: dowght
contains
  !----------------------------------------------------------------
  subroutine oc_dso_init(nzs_in,nzb_in,nmx_in,nmy_in &
  ,                      nhx_in,ohx,dhx &
  ,                      nhy_in,ohy,dhy, maxmem_in)
    integer,            intent(in) :: nzs_in,nzb_in
    integer,            intent(in) :: nmx_in,nmy_in
    integer,            intent(in) :: nhx_in,nhy_in
    real,               intent(in) :: ohx,dhx,ohy,dhy
    integer,            optional   :: maxmem_in
    integer                        :: ihx,ihy

    ocmaxmem=0; if(present(maxmem_in)) ocmaxmem=maxmem_in

    nzs=nzs_in
    nzb=nzb_in
    nmx=nmx_in
    nmy=nmy_in
    nhx=nhx_in
    nhy=nhy_in

    allocate(image(nmx,nmy,nhx,nhy,nzs))
    allocate(   hh(        nhx,nhy    ))
    allocate(   hx(        nhx        ))
    allocate(   hy(            nhy    ))

    do ihy=1,nhy
       hy(ihy)=ohy + (ihy-1)*dhy
    end do

    do ihx=1,nhx
       hx(ihx)=ohx + (ihx-1)*dhx
    end do

    do ihy=1,nhy
       do ihx=1,nhx
          hh(ihx,ihy) = sqrt(hy(ihy)**2 + hx(ihx)**2)
       end do
    end do
    hh = hh / maxval(hh)

    deallocate(hx,hy)

  end subroutine oc_dso_init
  !----------------------------------------------------------------
  integer function oc_dso(adj,add, x_,yy_, op1,op2,op3,op4,op5,op6,op7,op8,op9) result(stat)
    logical,               intent(in)  :: adj,add
    character(len=128),    intent(in)  :: x_,yy_
    integer, external,     optional    :: op1,op2,op3,op4,op5,op6,op7,op8,op9
    type(fileinfo)                     :: filex,filey
    integer                            :: izb,izs,imx,imy,ihx,ihy

    call oc_adjnull(adj,add,x_,yy_)

    call oc_allocatefile(filex, x_,  ocmaxmem); call sseek(x_, 0,0)
    call oc_allocatefile(filey, yy_, ocmaxmem); call sseek(yy_,0,0)
    call oc_checksimilarity(filex,filey,"oc_dso")

    do izb=1,nzb
       if(adj) then
          call sep_read (image, yy_)
       else
          call sep_read (image, x_ )
       end if

       do izs=1,nzs
          do imy=1,nmy
             do imx=1,nmx
                do ihy=1,nhy
                   do ihx=1,nhx
                      image(imx,imy,ihx,ihy,izs) = &
                      image(imx,imy,ihx,ihy,izs) * hh(ihx,ihy)
                   end do
                end do
             end do
          end do
       end do

       if(adj) then
          call sep_write(image,  x_)
       else
          call sep_write(image, yy_)
       end if
    end do

  end function oc_dso
  !----------------------------------------------------------------
  subroutine oc_dsowgh_init(nzs,nzb,nmx,nmy &
  ,                      nhx,ohx,dhx &
  ,                      nhy,ohy,dhy, maxmem, w_)
    integer,            intent(in) :: nzs,nzb
    integer,            intent(in) :: nmx,nmy
    integer,            intent(in) :: nhx,nhy
    real,               intent(in) :: ohx,dhx,ohy,dhy
    integer,            optional   :: maxmem
    character(len=128), optional   :: w_

    call oc_dso_init(nzs,nzb,nmx,nmy,nhx,ohx,dhx,nhy,ohy,dhy, maxmem)

    if(present(w_)) then
       call oc_weight_init(w_,maxmem)
    end if

  end subroutine oc_dsowgh_init
  !----------------------------------------------------------------
  integer function oc_dsowgh(adj,add, x_,yy_, op1,op2,op3,op4,op5,op6,op7,op8,op9) result(stat)
    logical,               intent(in)  :: adj,add
    character(len=128),    intent(in)  :: x_,yy_
    integer, external,     optional    :: op1,op2,op3,op4,op5,op6,op7,op8,op9
    character(len=128)                 :: t_,name
    type(fileinfo)                     :: filex,filey
    
    call oc_adjnull(adj,add,x_,yy_)
    
    call oc_allocatefile(filex, x_,  ocmaxmem); call sseek(x_, 0,0)
    call oc_allocatefile(filey, yy_, ocmaxmem); call sseek(yy_,0,0)
    call oc_checksimilarity(filex,filey, "oc_dsowgh")

    name='_t.dsowgh'; t_ = oc_clone(T, x_,name,ocmaxmem)
    call oc_assign(t_,0.,ocmaxmem)

    if(adj) then
       stat = oc_dso   (T,  F,t_,yy_)
       stat = oc_weight(T,add,x_, t_)
    else
       stat = oc_weight(F,  F,x_, t_)
       stat = oc_dso   (F,add,t_,yy_)
    end if

    call oc_deallocatefile(filex); call sseek(x_, 0,0)
    call oc_deallocatefile(filey); call sseek(yy_,0,0)
    stat = 0
  end function oc_dsowgh
  !----------------------------------------------------------------
end module oc_dso_mod
