module mva_scater_mod
  implicit none

  complex, dimension(:,:),     pointer, private :: P
  real   , dimension(:,:),     pointer, private :: S
  real                                , private :: dx, dy, dz
  integer                             , private :: nx, ny
contains

  integer function mva_scater_init(P_in, S_in, dx_in, dy_in, dz_in, nx_in, ny_in) result(stat)
    complex, dimension(:,:), target :: P_in
    real   , dimension(:,:), target :: S_in
    real                            :: dx_in, dy_in, dz_in
    integer                         :: nx_in, ny_in
 
    P => P_in
    S => S_in
    
    dx = dx_in
    dy = dy_in
    dz = dz_in
    nx = nx_in
    ny = ny_in

    stat = 0
  end function 


  integer function mva_scater_clos() result(stat)
    nullify(P)
    nullify(S)
    dx = 0.; dy = 0.; dz = 0.
    nx = 0 ; ny = 0
    stat = 0
  end function

  !!.. forward: dP = G (P,s) dS
  !!.. adjoint: dS = G'(P,s) dP
  !!.. G(P,s) = idz(dkz/ds)P
  !!------------------------------------------------------------------------------------
  !!.. dS    : perturbed slowness
  !!.. dP    : perturbed wavefield
  !!.. G(P,s): the scattering wavefield, it is a function of the background wavefield P
  integer function mva_scater_oper(adj, add, dS, dP, ws) result(stat)
    logical                 :: adj, add 
    real   , dimension(:,:) :: dS
    complex, dimension(:,:) :: dP
    real                    :: ws

    real                    :: dx2, w2, S2
    complex                 :: iwdz, lapxo
    integer                 :: ix

    iwdz = cmplx(0., 1.)*ws*dz
    dx2  = dx*dx
    w2   = ws*ws

    if (.not. add) then
       if (adj) then
           dS = 0.
       else
           dP = 0.
       end if
    end if

    do ix=1, nx
       S2    = S(ix,1)*S(ix,1) 
       if (ix == 1) then
          lapxo = (P(ix+1,1)-P(ix,1))/dx
       else if (ix == nx) then
          lapxo = (P(ix,1)-P(ix-1,1))/dx
       else
          lapxo = (P(ix-1,1)-2.*P(ix,1)+P(ix+1,1))/dx2
       end if

       if (adj) then
          dS(ix,1) = dS(ix,1) + conjg(iwdz*(P(ix,1)+0.5/w2/S2*lapxo))*dP(ix,1)
       else
          dP(ix,1) = dP(ix,1) +       iwdz*(P(ix,1)+0.5/w2/S2*lapxo) *dS(ix,1)
       end if 
    end do
    stat = 0
  end function


end module
