module wei_kmap
  !! K-domain mapping

  use wei_pro
  use wei_type
  !---------------
  implicit none

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  type(wekmap), public :: kmap
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!

contains
  !----------------------------------------------------------------
  subroutine weikmap_init()
    integer :: imx,ihx,imy,ihy
    integer :: er
    character(len=128), parameter :: unit='weikmap_init'
#ifdef DEBUG
    if(wpro%dbg>0) call in(unit)
#endif

    call weimem_add(4,wpro%amx%n,unit)
    allocate(kmap%kmx(wpro%amx%n),stat=er)
    if(er/=0) call seperr('cant allocate kmx')
    call klinemap(kmap%kmx,wpro%amx%n,wpro%amx%d)

    call weimem_add(4,wpro%amy%n,unit)
    allocate(kmap%kmy(wpro%amy%n),stat=er)
    if(er/=0) call seperr('cant allocate kmy')
    call klinemap(kmap%kmy,wpro%amy%n,wpro%amy%d)

    call weimem_add(4,wpro%ahx%n,unit)
    allocate(kmap%khx(wpro%ahx%n),stat=er)
    if(er/=0) call seperr('cant allocate khx')
    call klinemap(kmap%khx,wpro%ahx%n,wpro%ahx%d)

    call weimem_add(4,wpro%ahy%n,unit)
    allocate(kmap%khy(wpro%ahy%n),stat=er)
    if(er/=0) call seperr('cant allocate khy')
    call klinemap(kmap%khy,wpro%ahy%n,wpro%ahy%d)

    !! ksx(amx,ahx), krx(amx,ahx)
    call weimem_add(4,wpro%amx%n*wpro%ahx%n,unit)
    allocate(kmap%ksx(wpro%amx%n,wpro%ahx%n),stat=er)
    if(er/=0) call seperr('cant allocate ksx')

    call weimem_add(4,wpro%amx%n*wpro%ahx%n,unit)
    allocate(kmap%krx(wpro%amx%n,wpro%ahx%n),stat=er)
    if(er/=0) call seperr('cant allocate krx')
    !! ksx = kmx - khx
    !! krx = kmx + khx
    do ihx=1,wpro%ahx%n
       do imx=1,wpro%amx%n
          kmap%ksx(imx,ihx) = kmap%kmx(imx)-kmap%khx(ihx)
          kmap%krx(imx,ihx) = kmap%kmx(imx)+kmap%khx(ihx)
       end do
    end do

    !! ksy(amy,ahy), kry(amy,ahy)
    call weimem_add(4,wpro%amy%n*wpro%ahy%n,unit)
    allocate(kmap%ksy(wpro%amy%n,wpro%ahy%n),stat=er)
    if(er/=0) call seperr('cant allocate ksy')

    call weimem_add(4,wpro%amy%n*wpro%ahy%n,unit)
    allocate(kmap%kry(wpro%amy%n,wpro%ahy%n),stat=er)
    if(er/=0) call seperr('cant allocate kry')
    !! ksy = kmy - khy
    !! kry = kmy + khy
    do ihy=1,wpro%ahy%n
       do imy=1,wpro%amy%n
          kmap%ksy(imy,ihy) = kmap%kmy(imy)-kmap%khy(ihy)
          kmap%kry(imy,ihy) = kmap%kmy(imy)+kmap%khy(ihy)
       end do
    end do

#ifdef DEBUG
    if(wpro%dbg>0) call out(unit)
#endif
  end subroutine weikmap_init
  !----------------------------------------------------------------
  subroutine weikmap_report()

    write(0,*)'        ************************************************'
    if(associated(kmap%kmx)) then
       write(0,*)"        *    size(kmx)",size(kmap%kmx),"/",wpro%amx%n
    end if
    if(associated(kmap%kmy)) then
       write(0,*)"        *    size(kmy)",size(kmap%kmy),"/",wpro%amy%n
    end if
    if(associated(kmap%khx)) then
       write(0,*)"        *    size(khx)",size(kmap%khx),"/",wpro%ahx%n
    end if
    if(associated(kmap%khy)) then
       write(0,*)"        *    size(khy)",size(kmap%khy),"/",wpro%ahy%n
    end if
    if(associated(kmap%ksx)) then
       write(0,*)"        *    size(ksx)",size(kmap%ksx),"/",wpro%amx%n*wpro%ahx%n
    end if
    if(associated(kmap%krx)) then
       write(0,*)"        *    size(krx)",size(kmap%krx),"/",wpro%amx%n*wpro%ahx%n
    end if
    if(associated(kmap%ksy)) then
       write(0,*)"        *    size(ksy)",size(kmap%ksy),"/",wpro%amy%n*wpro%ahy%n
    end if
    if(associated(kmap%kry)) then
       write(0,*)"        *    size(kry)",size(kmap%kry),"/",wpro%amy%n*wpro%ahy%n
    end if
    write(0,*)'        ************************************************'

  end subroutine weikmap_report
  !----------------------------------------------------------------
  subroutine klinemap(k,n,d)
    real, dimension(:), pointer :: k
    integer, intent(in) :: n
    real,    intent(in) :: d
    integer             :: i,m

    k=0
    if(n/=1) then
       m = n/2

       do i=1,m
          k(i  ) =    i-1
          k(i+m) = -m+i-1
       end do

    end if

    k = k / (n*d)       !! 1/length
    k = 2 * wpro%pi * k !! radians
  end subroutine klinemap
  !----------------------------------------------------------------
end module wei_kmap







