module wem_hespos_mod
  implicit none

contains

  integer function wem_hespos_init() result(stat)
  end function

  integer function wem_hespos_oper() result(stat)
    do iz_offd=1, hgrid%az_offd%n; hz_offd = hgrid%az_offd%o+(iz_offd-1)*hgrid%az_offd%d; jz_offd = nint(hz_offd/pgrid%az%d)
    do iy_offd=1, hgrid%ay_offd%n; hy_offd = hgrid%ay_offd%o+(iy_offd-1)*hgrid%ay_offd%d; jy_offd = nint(hy_offd/pgrid%ay%d)
    do ix_offd=1, hgrid%ax_offd%n; hx_offd = hgrid%ax_offd%o+(ix_offd-1)*hgrid%ax_offd%d; jx_offd = nint(hx_offd/pgrid%ax%d)

       call hespos_fileio(hess, hgrid, ogrid, ix_offd, iy_offd, iz_offd, filetag, "reed")

       do iz_modl=1, ogrid%az%n; index_z_modl = max(1, min(indz+iz_modl-1,pgrid%az%n))      ; z_modl = pgrid%az%o + (index_z_modl-1)*pgrid%az%d
                                 index_z_offd = max(1, min(index_z_modl+jz_offd,pgrid%az%n)); z_offd = pgrid%az%o + (index_z_offd-1)*pgrid%az%d
       do iy_modl=1, ogrid%ay%n; index_y_modl = max(1, min(indy+iy_modl-1,pgrid%ay%n))      ; y_modl = pgrid%ay%o + (index_y_modl-1)*pgrid%ay%d
                                 index_y_offd = max(1, min(index_y_modl+jy_offd,pgrid%ay%n)); y_offd = pgrid%ay%o + (index_y_offd-1)*pgrid%ay%d
       do ix_modl=1, ogrid%ax%n; index_x_modl = max(1, min(indx+ix_modl-1,pgrid%ax%n))      ; x_modl = pgrid%ax%o + (index_x_modl-1)*pgrid%ax%d
                                 index_x_offd = max(1, min(index_x_modl+jx_offd,pgrid%ax%n)); x_offd = pgrid%ax%o + (index_x_offd-1)*pgrid%ax%d

          greens_sou =       wfds(index_x_modl,index_y_modl,index_z_modl)  &
                     * conjg(wfds(index_x_offd,index_y_offd,index_z_offd))

          select case (hgrid%hestyp)
            case (1) !!.. source intensity, assuming receiver side Green's function to be 1
                 greens_rec = cmplx(1., 0.)
            case (2) !!.. approximate the receiver side Green's function as a aperture dependent constant (not ready yet)
                 if (pgrid%ay%n>1 .and. pgrid%ax%n>1) then
                    wt_modl = rec_green3d(x_modl, y_modl, pgrid%gxmin, pgrid%gxmax, pgrid%gymin, pgrid%gymax, z_modl)
                    wt_offd = rec_green3d(x_offd, y_offd, pgrid%gxmin, pgrid%gxmax, pgrid%gymin, pgrid%gymax, z_offd)
                 else
                    wt_modl = rec_green2d(x_modl, pgrid%gxmin, pgrid%gxmax, z_modl)
                    wt_offd = rec_green2d(x_offd, pgrid%gxmin, pgrid%gxmax, z_offd)
                 end if
                 greens_rec = wt_modl
            case (3) !!.. Hessian via phase-encoding 
                 greens_rec = conjg(wfdr(index_x_modl,index_y_modl,index_z_modl)) &
                            *       wfdr(index_x_offd,index_y_offd,index_z_offd)
          end select

          hess(ix_modl,iy_modl,iz_modl) = hess(ix_modl,iy_modl,iz_modl) + real(greens_sou*greens_rec)

       end do
       end do
       end do

       call hesapp_fileio(hess, hgrid, ogrid, ix_offd, iy_offd, iz_offd, filetag, "rite")

    end do
    end do
    end do

  end function

end module
