!!.. performs local windowing and its adjoint
module hes_locwin_mod
  implicit none


  integer, private :: dim1, dim2, dim3

 

contains

  integer function hes_locwin_init() result(stat)
  end function



  integer function hes_locwin_oper() result(stat)

  end function




  subroutine fdecon_kernel(modl, data, hess, imag)
    real, dimension(:,:,:) :: modl
    real, dimension(:,:,:) :: data
    real, dimension(:,:,:) :: hess
    real, dimension(:,:,:) :: imag

    integer                :: ix_modl, iy_modl, iz_modl
    integer                :: ix_offd, iy_offd, iz_offd
    integer                :: jx_offd, jy_offd, jz_offd
    integer                :: index_x_modl, index_y_modl, index_z_modl
    real                   :: x_offd, y_offd, z_offd

    call read_data

    do iz_modl=1, nz_modl 
    do iy_modl=1, ny_modl
    do ix_modl=1, nx_modl

       do iz_offd=1, hgrid%az_offd; 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 = iz_modl+jz_offd                            ; if(index_z_offd<1.or.index_z_offd>dim3) cycle

       do iy_offd=1, hgrid%ay_offd; 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 = iy_modl+jy_offd                            ; if(index_y_offd<1.or.index_y_offd>dim2) cycle

       do ix_offd=1, hgrid%ax_offd; 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 = ix_modl+jx_offd                            ; if(index_x_offd<1.or.index_x_offd>dim1) cycle

          data(ix_offd,iy_offd,iz_offd) = imag(index_x_offd,index_y_offd,index_z_offd)

       end do
       end do
       end do

       call deconv_fileio(hess, hgrid, ix_modl, iy_modl, iz_modl, filetag, "reed")

       call fft3d(data, +1)
       call fft3d(hess, +1)

       modl = data / (hess+eps)

       call fft3d(modl, -1) 

       call assign_image_back

    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



end module
