module mva_cpsfou_mod
  use  base_slofft_mod
  implicit none

  real, dimension(:), pointer         , private :: kx, ky
  real                                , private :: dx, dy, dz
  integer                             , private :: nx, ny
contains

  integer function mva_cpsfou_init(kx_in, ky_in, dx_in, dy_in, dz_in, nx_in, ny_in) result(stat)
    real, dimension(:), target :: kx_in, ky_in
    real                       :: dx_in, dy_in, dz_in
    integer                    :: nx_in, ny_in
 
    kx => kx_in
    ky => ky_in
    dx =  dx_in
    dy =  dy_in
    dz =  dz_in
    nx =  nx_in
    ny =  ny_in

    stat = 0
  end function 


  integer function mva_cpsfou_clos() result(stat)
    dx = 0.; dy = 0.; dz = 0.
    nx = 0 ; ny = 0
    stat = 0
  end function

  !!.. G = G(P,S)
  integer function mva_cpsfou_oper(G, P, S, ws) result(stat)
    complex, dimension(:,:) :: G, P
    real   , dimension(:,:) :: S
    real                    :: ws

    real                    :: dx2, w2, S2, mk2, kp2
    complex                 :: iwdz, lapxo
    integer                 :: ikx, iky

    iwdz = cmplx(0., 1.)*ws*dz
    dx2  = dx*dx
    w2   = ws*ws
    S2   = (sum(S)/nx/ny)**2

    call FFT2DX(P, +1)

    do iky=1, ny
    do ikx=1, nx
       mk2 = kx(ikx)**2 + ky(iky)**2
       kp2 = 1-mk2/w2/S2 
       if (kp2>0.0001) then
          G(ikx,iky) = iwdz*P(ikx,iky)/sqrt(kp2)
       else
          G(ikx,iky) = 0.
       end if
    end do
    end do

    call FFT2DX(P, -1)
    call FFT2DX(G, -1)
    stat = 0
  end function


end module
