module wei_kmap_mod
  use sep3d_struct_mod
  use wei_wavespace_type_mod
  !! K-domain mapping
  !---------------
  implicit none

contains
  !----------------------------------------------------------------
  subroutine weikmap_init(wave)
    integer :: imx,ihx, imy,ihy
    integer :: er
    type(wspace) :: wave

#ifdef DEBUG
   write(0,*) "IN WEIKMAP_INIT"
#endif

    allocate(wave%kmx(wave%npad(1)),stat=er)
    if(er/=0) call seperr('cant allocate kmx')
    call klinemap(wave%kmx,wave%npad(1),wave%d(1),0)

    allocate(wave%kmy(wave%npad(2)),stat=er)
    if(er/=0) call seperr('cant allocate kmy')
    call klinemap(wave%kmy,wave%npad(2),wave%d(2),0)

    allocate(wave%khx(wave%npad(3)),stat=er)
    if(er/=0) call seperr('cant allocate khx')
    call klinemap(wave%khx,wave%npad(3),wave%d(3),0)

    allocate(wave%khy(wave%npad(4)),stat=er)
    if(er/=0) call seperr('cant allocate khy')
    call klinemap(wave%khy,wave%npad(4),wave%d(4),0)


    !! ksx(amx,ahx), krx(amx,ahx)
    allocate(wave%ksx(wave%npad(1), wave%npad(3)),stat=er)
    if(er/=0) call seperr('cant allocate ksx')
    allocate(wave%krx(wave%npad(1), wave%npad(3)),stat=er)
    if(er/=0) call seperr('cant allocate krx')

    !! ksy(amy,ahy), kry(amy,ahy)
    allocate(wave%ksy(wave%npad(2), wave%npad(4)),stat=er)
    if(er/=0) call seperr('cant allocate ksy')
    allocate(wave%kry(wave%npad(2), wave%npad(4)),stat=er)
    if(er/=0) call seperr('cant allocate kry')

    call form_ksx_kry(wave)

#ifdef DEBUG
    write(0,*) "OUT KMAP INIT"
#endif
  end subroutine weikmap_init

  subroutine form_ksx_kry(wave)
  type(wspace) :: wave
  integer :: ihx,ihy,imx,imy
#ifdef DEBUG
   write(0,*) "IN FORM KSX_KRY"
#endif
    !! ksx = kmx - khx
    !! krx = kmx + khx
    do ihx=1,wave%npad(3)
       do imx=1,wave%npad(1)
          wave%ksx(imx,ihx) = wave%kmx(imx)-wave%khx(ihx)
          wave%krx(imx,ihx) = wave%kmx(imx)+wave%khx(ihx)
       end do
    end do

    !! ksy = kmy - khy
    !! kry = kmy + khy
    do ihy=1,wave%npad(4)
       do imy=1,wave%npad(2)
          wave%ksy(imy,ihy) = wave%kmy(imy)-wave%khy(ihy)
          wave%kry(imy,ihy) = wave%kmy(imy)+wave%khy(ihy)
       end do
    end do
#ifdef DEBUG
   write(0,*) "OUT FORM KSX_KRY"
#endif

  end subroutine


  !----------------------------------------------------------------
  subroutine klinemap(k,n,d,shift)
    real, dimension(:), pointer :: k
    integer, intent(in) :: n,shift
    real,    intent(in) :: d
    integer             :: i,m
#ifdef DEBUG
    write(0,*) "IN KLINEMAP"
#endif
    k=0.
    if(n/=1) then
       m = n/2
       do i=1,m
          k(i+shift  ) =    i-1
          k(i+m+shift) = -m+i-1
       end do

!       do i=1,n
!          if(i <= m) then
!            k(i+shift) = i -1
!          else
!            k(i+shift) = i - n -1
!          end if
!       end do


    else
       k=0. !! make sure k=0 if n=1
    end if

    k = k / (n*d)      !! 1/length
    k = 2 * atan(1.)*4 * k !! radians
#ifdef DEBUG
    write(0,*) "OUT KLINEMAP"
#endif
  end subroutine klinemap

end module
