module OW_imp_scatter

  use OW_type
  use OW_parms
  use OW_fftw

  implicit none
  
contains
  !---------------------------------------------
  ! implicit scattering operator derived from 45 degree extrapolator
  !---------------------------------------------

subroutine owimpscatter(iwa,ddws,bwr,slow,dslow,sign,adj)

  complex                                  :: ddws(:), bwr(:)
  real                                     :: slow(:), dslow(:)
  complex, dimension(kxm%n)                :: A, B, C, U, aa, bb, cc, cdslow
  integer                                  :: ix, iwa
  real                                     :: dx2
  integer                                  :: sign                ! +1: source ; -1: receiver
  logical                                  :: adj                 !  0: forward;  1: adjoint


  if (rec%w%all(iwa).ne.0.) then
        A = cmplx(1., -rec%w%all(iwa)*slow*img%z%d)
        A = 2 * A*A * rec%w%all(iwa)**2 * slow**2
        B = cmplx(1., -3. * rec%w%all(iwa)*slow*img%z%d) * \
        cmplx(1., rec%w%all(iwa)*slow*img%z%d)
        

        dx2 = img%xm%d**2

        if (adj) then          ! Get ds from dp
        
           if (sign==1) then   
              A = conjg(A)
              B = conjg(B)
           end if
           C = cmplx(0., 4.) * img%z%d * ddws
           aa = -B / dx2
           bb = A + 2*B/dx2
           cc = aa
           cdslow = cmplx(dslow,0)
           U = bwr * C
           call ctridag(aa,bb,cc,U,cdslow,kxm%n)
           
        else                  ! Get dp from ds

           if (sign==-1) then
              A = conjg(A)
              B = conjg(B)
           end if
           C = cmplx(0., 4.) * img%z%d * dslow
           
           aa = -B / dx2
           bb = A + 2*B/dx2
           cc = a
           
           U = bwr * C
           call ctridag(a,b,c,U,ddws,kxm%n)
           
        end if
  else
     ddws=0.
  end if

end subroutine owimpscatter


  !!---------------------------------------------------------------------------
  !!.. complex tridiagonal solver
  !!---------------------------------------------------------------------------
  !!.. a(n), sub-diagonal elemnets, only elements from 2:n are used
  !!.. b(n), diagonal elemetns
  !!.. c(n), super-diagonal elements, only elements from 1:n-1 are used
  !!
  !!   /  b(1), c(1)                                     \ / u(  1) \     / r(  1) \ 
  !!   |  a(2), b(2), c(2)                               | | u(  2) |     | r(  2) |
  !!   |        a(3), b(3), c(3)                         | | u(  3) |  =  | r(  3) |
  !!   |              ... ...                            | |  ...   |     |  ...   |
  !!   |                          a(n-1), b(n-1), c(n-1) | | u(n-1) |     | r(n-1) |
  !!   \                                  a(  n), b(  n) / \ u(  n) /     \ r(  n) /
  
  subroutine ctridag(a,b,c,r,u,n)
    integer            :: n
    complex            :: a(n), b(n), c(n), r(n), u(n)
    integer            :: j
    integer, parameter :: NMAX = 5000
    complex            :: bet, gam(NMAX)
    bet  = b(1)
    u(1) = r(1)/bet
    do j=2, n
       gam(j) = c(j-1)/bet
       bet    = b(j)-a(j)*gam(j)
       u(j) = (r(j)-a(j)*u(j-1))/bet
    end do
    do j=n-1,1,-1
       u(j) = u(j) - gam(j+1)*u(j+1)
    end do
    return
  end subroutine ctridag

end module OW_imp_scatter
