module wei_lint

  use wei_util
  use wei_process

  !---------------
  implicit none

  integer, private :: nx,ny
  real,    private :: ox,oy
  real,    private :: dx,dy

  real,    private :: xmin, xmax
  real,    private :: ymin, ymax
contains
  !----------------------------------------------------------------  
  subroutine weilint_init()

    ox = pro%amx%o
    oy = pro%amy%o

    dx = pro%amx%d
    dy = pro%amy%d

    nx = pro%amx%n
    ny = pro%amy%n

    xmin = ox
    ymin = oy

    xmax = ox + (nx-1)*dx
    ymax = oy + (ny-1)*dy

  end subroutine weilint_init
  !----------------------------------------------------------------
  subroutine weilint(adj,xc,yc,vals,cube)
    logical, intent(in)                :: adj
    real,    intent(in)                :: xc,yc
    complex, dimension(:)              :: vals
    complex, dimension(:,:,:), pointer :: cube

    real    :: f,fx,fy,gx,gy
    integer :: i,ix0,ix1,iy0,iy1
    logical :: inmapx, inmapy

    if(nx==1) then
       ix0=1
       ix1=1
       fx=0.
       gx=1.
       inmapx=.true.
    else
       if(xc<xmin .or. xc>=xmax) then
          inmapx=.false.
       else
          inmapx=.true.
          f  = (xc-ox)/dx
          i  = f
          ix0= 1+i
          ix1= 1+i+1
          fx = f-i
          gx = 1.-fx
       end if
    end if

    if(ny==1) then
       iy0=1
       iy1=1
       fy=0.
       gy=1.
       inmapy=.true.
    else
       if(yc<ymin .or. yc>=ymax) then
          inmapy=.false.
       else
          inmapy=.true.
          f  = (yc-oy)/dy
          i  = f
          iy0= 1+i
          iy1= 1+i+1
          fy = f-i
          gy = 1.-fy
       end if
    end if

    if(inmapx .and. inmapy) then
       if(adj) then !! 4 pull
          cube(ix0,iy0,:) = cube(ix0,iy0,:) + gx * gy * vals(:)
          cube(ix1,iy0,:) = cube(ix1,iy0,:) + fx * gy * vals(:)
          cube(ix0,iy1,:) = cube(ix0,iy1,:) + gx * fy * vals(:)
          cube(ix1,iy1,:) = cube(ix1,iy1,:) + fx * fy * vals(:)
       else         !! 4 push
          vals(:)         =                   gx * gy * cube(ix0,iy0,:) 
          vals(:)         = vals(:)         + fx * gy * cube(ix1,iy0,:)
          vals(:)         = vals(:)         + gx * fy * cube(ix0,iy1,:)
          vals(:)         = vals(:)         + fx * fy * cube(ix1,iy1,:)
       end if

    else
       if(adj) then
          ! do nothing
       else
          vals(:)=0
       end if
    end if
  end subroutine weilint
  !----------------------------------------------------------------
end module wei_lint
