!!.. Fourier domain non-stationary deconvolution
module hes_fdecon_mod2
  use  base_ioutil_mod
  use  base_struct_mod
  use  base_fftwv3_mod
  use  base_slofft_mod
  use  wem_wnbmap_mod
!  use  hes_locprd_mod
  use  csolver_tiny_mod
  use  ccgstep_mod
!  use  solver_tiny_mod
!  use  cgstep_mod
  implicit none

  type(axis)                            , private :: akx, aky, akz
  type(hes_grid)                        , private :: hgrid
  character(len=256)                    , private :: filetag

  integer*8                             , private :: plan_ldat_3dxF
  integer*8                             , private :: plan_lhes_3dxF
  integer*8                             , private :: plan_lmod_3dxI

  complex, dimension(:,:,:), allocatable, private :: ldat, lmod, lhes
  real   , dimension(:,:,:), allocatable, private :: hess
  real   , dimension(:)    , allocatable, private :: kx, ky, kz

  integer                               , private :: nkx, nky, nkz
  integer                               , private :: nx_wind, ny_wind, nz_wind
  integer, external                     , private :: sreed_window, srite_window
 

contains

  integer function hes_fdecon_init(hgrid_in, filetag_in) result(stat)
    type(hes_grid)   :: hgrid_in
    character(len=*) :: filetag_in
 
    stat = 1
    
    hgrid   = hgrid_in

    filetag = filetag_in

    nkx = hgrid%ax_offd%n + hgrid%nx_offd_pad_beg + hgrid%nx_offd_pad_end
    nky = hgrid%ay_offd%n + hgrid%ny_offd_pad_beg + hgrid%ny_offd_pad_end
    nkz = hgrid%az_offd%n + hgrid%nz_offd_pad_beg + hgrid%nz_offd_pad_end

write(0,*)"nkx,nky,nkz", nkx, nky, nkz

    allocate(ldat(nkx,nky,nkz))
    allocate(lmod(nkx,nky,nkz))
    allocate(lhes(nkx,nky,nkz))
    allocate(hess(hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n))
    allocate(kx(nkx),ky(nky),kz(nkz))


    nx_wind = hgrid%ax_modl%n / hgrid%ax_offd%n
    ny_wind = hgrid%ay_modl%n / hgrid%ay_offd%n
    nz_wind = hgrid%az_modl%n / hgrid%az_offd%n

write(0,*)"nx_wind, ny_wind, nz_wind", nx_wind, ny_wind, nz_wind

akx%n = nkx; aky%n = nky; akz%n = nkz
akx%o = 0. ; aky%o = 0. ; akz%o = 0.
akx%d = 1. ; aky%d = 1. ; akz%d = 1.

    ldat = 0.; lmod = 0.; lhes = 0.
    kx   = 0.; ky   = 0.; kz   = 0.

!    call k_map(kx, nkx, hgrid%ax_offd%d)
!    call k_map(ky, nky, hgrid%ay_offd%d)
!    call k_map(kz, nkz, hgrid%az_offd%d)

    !!.. create the plan for the forward Fourier transform of the local data
    call sfftw_dft_3dx_init(plan_ldat_3dxF, ldat, nkx, nky, nkz, +1)

    !!.. create the plan for the forward Fourier transform of the local Hessian
    call sfftw_dft_3dx_init(plan_lhes_3dxF, lhes, nkx, nky, nkz, +1)

    !!.. create the plan for the inverse Fourier transform of the local image
    call sfftw_dft_3dx_init(plan_lmod_3dxI, lmod, nkx, nky, nkz, -1)

    stat = 0
  end function

  integer function hes_fdecon_clos() result(stat)
    deallocate(ldat)
    deallocate(lmod)
    deallocate(lhes)
    deallocate(hess)
    stat = 0
  end function

  integer function hes_fdecon_oper(modl, data) result(stat)
    real, dimension(:,:,:) :: modl  !!.. the totoal image size, output
    real, dimension(:,:,:) :: data  !!.. the totoal image size, input

    stat = 1

    call fdecon_kernel(modl, data)

    stat = 0
  end function




  subroutine fdecon_kernel(modl, data)
    real, dimension(:,:,:) :: modl  !!.. the totoal image size, output
    real, dimension(:,:,:) :: data  !!.. the totoal image size, input

    integer                :: ix_modl, iy_modl, iz_modl
    integer                :: ix_offd, iy_offd, iz_offd
    integer                :: jx_offd, jy_offd, jz_offd
    integer                :: jx_modl, jy_modl, jz_modl
    integer                :: ikx, iky, ikz
    integer                :: ix_wind, iy_wind, iz_wind
    integer                :: index_x_modl, index_y_modl, index_z_modl
    integer                :: index_x_offd, index_y_offd, index_z_offd

    integer                :: stat
    real                   :: x_offd, y_offd, z_offd
  call auxout("test_hes")
  call houaxis("test_hes", akx, 1)
  call houaxis("test_hes", aky, 2)
  call houaxis("test_hes", akz, 3)
  call houaxis("test_hes", hgrid%ax_modl, 4)
  call houaxis("test_hes", hgrid%ay_modl, 5)
  call houaxis("test_hes", hgrid%az_modl, 6)
  call auxputch("esize", "i", 8, "test_hes")

  call auxout("test_dat")
  call houaxis("test_dat", akx, 1)
  call houaxis("test_dat", aky, 2)
  call houaxis("test_dat", akz, 3)
  call houaxis("test_dat", hgrid%ax_modl, 4)
  call houaxis("test_dat", hgrid%ay_modl, 5)
  call houaxis("test_dat", hgrid%az_modl, 6)
  call auxputch("esize", "i", 8, "test_dat")

  call auxout("test_mod")
  call houaxis("test_mod", akx, 1)
  call houaxis("test_mod", aky, 2)
  call houaxis("test_mod", akz, 3)
  call houaxis("test_mod", hgrid%ax_modl, 4)
  call houaxis("test_mod", hgrid%ay_modl, 5)
  call houaxis("test_mod", hgrid%az_modl, 6)
  call auxputch("esize", "i", 8, "test_mod")


    modl = 0.

    do iz_wind=1, nz_wind; write(0,*)"iz_modl=",iz_wind 
    do iy_wind=1, ny_wind 
    do ix_wind=1, nx_wind

       ldat = cmplx(0.,0.)
       lhes = cmplx(0.,0.)
       lmod = cmplx(0.,0.)

       index_x_modl = (ix_wind-1)*hgrid%ax_offd%n+hgrid%ax_offd%n/2+1
       index_y_modl = (iy_wind-1)*hgrid%ay_offd%n+hgrid%ay_offd%n/2+1
       index_z_modl = (iz_wind-1)*hgrid%az_offd%n+hgrid%az_offd%n/2+1

       !!.. extract the locally windowd images
       do iz_offd=1, hgrid%az_offd%n; z_offd = hgrid%az_offd%o+(iz_offd-1)*hgrid%az_offd%d; jz_offd = nint(z_offd/hgrid%az_modl%d);
                                      index_z_offd = index_z_modl+jz_offd; if(index_z_offd<1.or.index_z_offd>hgrid%az_modl%n) cycle

       do iy_offd=1, hgrid%ay_offd%n; y_offd = hgrid%ay_offd%o+(iy_offd-1)*hgrid%ay_offd%d; jy_offd = nint(y_offd/hgrid%ay_modl%d);
                                      index_y_offd = index_y_modl+jy_offd; if(index_y_offd<1.or.index_y_offd>hgrid%ay_modl%n) cycle

       do ix_offd=1, hgrid%ax_offd%n; x_offd = hgrid%ax_offd%o+(ix_offd-1)*hgrid%ax_offd%d; jx_offd = nint(x_offd/hgrid%ax_modl%d);
                                      index_x_offd = index_x_modl+jx_offd; if(index_x_offd<1.or.index_x_offd>hgrid%ax_modl%n) cycle

          jx_modl = ix_offd + hgrid%nx_offd_pad_beg
          jy_modl = iy_offd + hgrid%ny_offd_pad_beg
          jz_modl = iz_offd + hgrid%nz_offd_pad_beg 

          ldat(jx_modl,jy_modl,jz_modl) = cmplx(data(index_x_offd,index_y_offd,index_z_offd),0.)

       end do
       end do
       end do




       !!.. read the off-diagonals of the Hessian
       call deconv_fileio(hess, hgrid, index_x_modl, index_y_modl, index_z_modl, filetag, "reed")

       lhes = 0.

       lhes(hgrid%nx_offd_pad_beg+1:hgrid%nx_offd_pad_beg+hgrid%ax_offd%n &
           ,hgrid%ny_offd_pad_beg+1:hgrid%ny_offd_pad_beg+hgrid%ay_offd%n &
           ,hgrid%nz_offd_pad_beg+1:hgrid%nz_offd_pad_beg+hgrid%az_offd%n ) = hess

       !!.. Fourier transform of the Hessian
       call sfftw_dft_3dx_oper(plan_lhes_3dxF, lhes, nkx, nky, nkz)
       
       !!.. Fourier transform of the local image
       call sfftw_dft_3dx_oper(plan_ldat_3dxF, ldat, nkx, nky, nkz)

       !!.. normalize the Hessian spectrum
       lhes = lhes/maxval(cabs(lhes))
       

       !!.. inversion in the Fourier domain
       lmod = ldat*conjg(lhes)!*lhes/(lhes*conjg(lhes)+0.000001)


       !!.. inverse Fourier transform
       call sfftw_dft_3dx_oper(plan_lmod_3dxI, lmod, nkx, nky, nkz) 


       do iz_offd=1, hgrid%az_offd%n; z_offd = hgrid%az_offd%o+(iz_offd-1)*hgrid%az_offd%d; jz_offd = nint(z_offd/hgrid%az_modl%d);
                                      index_z_offd = index_z_modl+jz_offd; if(index_z_offd<1.or.index_z_offd>hgrid%az_modl%n) cycle

       do iy_offd=1, hgrid%ay_offd%n; y_offd = hgrid%ay_offd%o+(iy_offd-1)*hgrid%ay_offd%d; jy_offd = nint(y_offd/hgrid%ay_modl%d);
                                      index_y_offd = index_y_modl+jy_offd; if(index_y_offd<1.or.index_y_offd>hgrid%ay_modl%n) cycle

       do ix_offd=1, hgrid%ax_offd%n; x_offd = hgrid%ax_offd%o+(ix_offd-1)*hgrid%ax_offd%d; jx_offd = nint(x_offd/hgrid%ax_modl%d);
                                      index_x_offd = index_x_modl+jx_offd; if(index_x_offd<1.or.index_x_offd>hgrid%ax_modl%n) cycle

          jx_modl = ix_offd + hgrid%nx_offd_pad_beg
          jy_modl = iy_offd + hgrid%ny_offd_pad_beg
          jz_modl = iz_offd + hgrid%nz_offd_pad_beg

          modl(index_x_offd,index_y_offd,index_z_offd) = modl(index_x_offd,index_y_offd,index_z_offd)+real(lmod(jx_modl,jy_modl,jz_modl))

       end do
       end do
       end do


    end do
    end do
    end do
    return
  end subroutine


  subroutine deconv_fileio(hess, hgrid, index_x_modl, index_y_modl, index_z_modl, filetag, flag)
    real, dimension(:,:,:) :: hess
    type(hes_grid)         :: hgrid
    integer                :: index_x_modl, index_y_modl, index_z_modl
    character(len=*)       :: filetag
    character(len=4)       :: flag

    integer, dimension(6)  :: n_cube, j_wind, n_wind, f_wind
    integer                :: n_dim_cube, err

    n_dim_cube = 6
    n_cube = (/hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n,hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n/)
    j_wind = (/1,1,1,1,1,1/)
    n_wind = (/hgrid%ax_offd%n,hgrid%ay_offd%n,hgrid%az_offd%n,1,1,1/)
    f_wind = (/0,0,0,index_x_modl-1,index_y_modl-1,index_z_modl-1/)

    if (flag == "rite") then
       err = srite_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,hess)
    else if (flag == "reed") then
       hess = 0.
       err = sreed_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,4,hess)
    end if

    return
  end subroutine


  subroutine deconv_fileio_test(hess, hgrid, index_x_modl, index_y_modl, index_z_modl, filetag, flag)
    complex, dimension(:,:,:) :: hess
    type(hes_grid)            :: hgrid
    integer                   :: index_x_modl, index_y_modl, index_z_modl
    character(len=*)          :: filetag
    character(len=4)          :: flag

    integer, dimension(6)     :: n_cube, j_wind, n_wind, f_wind
    integer                   :: n_dim_cube, err

    n_dim_cube = 6
    n_cube = (/nkx,nky,nkz,hgrid%ax_modl%n,hgrid%ay_modl%n,hgrid%az_modl%n/)
    j_wind = (/1,1,1,1,1,1/)
    n_wind = (/nkx,nky,nkz,1,1,1/)
    f_wind = (/0,0,0,index_x_modl-1,index_y_modl-1,index_z_modl-1/)

    if (flag == "rite") then
       err = srite_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,hess)
    else if (flag == "reed") then
       hess = 0.
       err = sreed_window(filetag,n_dim_cube,n_cube,n_wind,f_wind,j_wind,8,hess)
    end if

    return
  end subroutine




end module
