module wei_lint
  use wei_util
  use wei_pro
  use wei_par
  use wei_S
  !---------------
  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 = wpro%amx%o
    oy = wpro%amy%o

    dx = wpro%amx%d
    dy = wpro%amy%d

    nx = wpro%amx%n
    ny = wpro%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
    real    :: w00,w01,w10,w11,sumw


    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
          if( abs(f-nint(f))<1000*epsilon(f) ) f=nint(f)

          i  = floor(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
          if( abs(f-nint(f))<1000*epsilon(f) ) f=nint(f)

          i  = floor(f)
          iy0= 1+i
          iy1= 1+i+1
          fy = f-i
          gy = 1.-fy
       end if
    end if

    !! 00   10       --> x
    !!    .         |
    !! 01   11      v
    !!               y
    w00=(1.-fx)*(1.-fy)
    w10=(   fx)*(1.-fy)
    w01=(1.-fx)*(   fy)
    w11=(   fy)*(   fy)

    !! normalize the sum of the weights
    sumw = w00+w10+w01+w11
    w00 = w00 / sumw
    w10 = w10 / sumw
    w01 = w01 / sumw
    w11 = w11 / sumw

    if(inmapx .and. inmapy) then
       if(adj) then !! 4 pull
          cube(ix0,iy0,:) = cube(ix0,iy0,:) + w00 * vals(:)
          cube(ix1,iy0,:) = cube(ix1,iy0,:) + w10 * vals(:)
          cube(ix0,iy1,:) = cube(ix0,iy1,:) + w01 * vals(:)
          cube(ix1,iy1,:) = cube(ix1,iy1,:) + w11 * vals(:)
       else         !! 4 push
          vals(:)         =                   w00 * cube(ix0,iy0,:) 
          vals(:)         = vals(:)         + w10 * cube(ix1,iy0,:)
          vals(:)         = vals(:)         + w01 * cube(ix0,iy1,:)
          vals(:)         = vals(:)         + w11 * cube(ix1,iy1,:)
       end if
    else
       if(adj) then
          ! do nothing
       else
          vals(:)=0
       end if
    end if

  end subroutine weilint
  !----------------------------------------------------------------
  subroutine weilintc(adj,X,line,xc,yc)
    type(weslow), intent(inout) :: X
    logical,                   intent(in) :: adj
    complex, dimension(:,:,:), pointer    :: cube
    complex, dimension(:),     pointer    :: line
    real,                      intent(in) :: xc,yc
    real    :: xx,yy
    integer :: ix,iy, jx,jy
    real    :: fx,fy
    real    :: w00,w01,w10,w11, sumw
    real    :: kx,ky

    cube=>X%bin
    xx = xc
    yy = yc

    !! don't allow coordinates outside the slowness grid
    xx = min(xx,ox + (nx-1)*dx)
    xx = max(xx,ox)
    yy = min(yy,oy + (ny-1)*dy)
    yy = max(yy,oy)

    kx=(xx-ox) / dx
    ky=(yy-oy) / dy

!!$    !! DANGER: 1e-3 due to seplib float truncations
!!$    !! lower index
!!$    if( abs((1+kx)-nint(1+kx))<1e-3) then
!!$       ix= nint(1+kx)
!!$    else
!!$       ix=floor(1+kx)
!!$    end if
!!$    if( abs((1+ky)-nint(1+ky))<1e-3) then
!!$       iy= nint(1+ky)
!!$    else
!!$       iy=floor(1+ky)
!!$    end if

    ix=floor(1+kx)
    iy=floor(1+ky)

    !! keep indices in bounds
    ix = min(ix,nx)
    iy = min(iy,ny)
    ix = max(ix, 1); jx=min(ix+1,nx)
    iy = max(iy, 1); jy=min(iy+1,ny)

    !! linear weights
    fx = kx - (ix-1)
    fy = ky - (iy-1)
    fx = max(0., fx)
    fy = max(0., fy)

    !! 00   10       --> x
    !!    .         |
    !! 01   11      v
    !!               y
    w00=(1.-fx)*(1.-fy)
    w10=(   fx)*(1.-fy)
    w01=(1.-fx)*(   fy)
    w11=(   fy)*(   fy)

    !! normalize the sum of the weights
    sumw = w00+w10+w01+w11
    w00 = w00 / sumw
    w10 = w10 / sumw
    w01 = w01 / sumw
    w11 = w11 / sumw

    if(adj) then !! push
       cube(ix,iy,:) = cube(ix,iy,:) + w00 * line
       cube(jx,iy,:) = cube(jx,iy,:) + w10 * line
       cube(ix,jy,:) = cube(ix,jy,:) + w01 * line
       cube(jx,jy,:) = cube(jx,jy,:) + w11 * line
    else         !! pull
       line    = &
       w00*real(cube(ix,iy,:)) + &
       w10*real(cube(jx,iy,:)) + &
       w01*real(cube(ix,jy,:)) + &
       w11*real(cube(jx,jy,:))
    end if
  end subroutine weilintc

end module wei_lint



