module wei_kmap
  !! K-domain mapping

  use wei_process
  use wei_type
  !---------------
  implicit none
  integer, private, save :: nhx_use,nhy_use

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

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

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

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

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

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

    nhx_use=pro%ahx%n
    nhy_use=pro%ahy%n

    !! ksx(amx,ahx), krx(amx,ahx)
    allocate(kmap%ksx(pro%amx%n, pro%ahx%n),stat=er)
    if(er/=0) call seperr('cant allocate ksx')
    allocate(kmap%krx(pro%amx%n, pro%ahx%n),stat=er)
    if(er/=0) call seperr('cant allocate krx')

    !! ksy(amy,ahy), kry(amy,ahy)
    allocate(kmap%ksy(pro%amy%n, pro%ahy%n),stat=er)
    if(er/=0) call seperr('cant allocate ksy')
    allocate(kmap%kry(pro%amy%n, pro%ahy%n),stat=er)
    if(er/=0) call seperr('cant allocate kry')

    call form_ksx_kry()

#ifdef DEBUG
    if(pro%debug) call out(unit)
#endif
  end subroutine weikmap_init

  subroutine redo_kmaps(x_beg,x_end,y_beg,y_end)
   integer,intent(in) :: x_beg,x_end,y_beg,y_end
   integer :: nx,ny
   logical :: redo
    redo=.false.
    nx=x_end-x_beg+1
    ny=y_end-y_beg+1
    if(nx .ne. nhx_use) then
      nhx_use=nx; redo=.true.
       call klinemap(kmap%khx,nx,pro%ahx%d,x_beg-1)
    end if
    if(ny .ne. nhy_use) then
       call klinemap(kmap%khy,ny,pro%ahy%d,y_beg-1)
      nhy_use=ny; redo=.true.
    end if
    if(redo) then
      call form_ksx_kry()
    end if
   

  end subroutine 

  subroutine form_ksx_kry()
  integer :: ihx,ihy,imx,imy
    !! ksx = kmx - khx
    !! krx = kmx + khx
    do ihx=1,pro%ahx%n
       do imx=1,pro%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 = kmy - khy
    !! kry = kmy + khy
    do ihy=1,pro%ahy%n
       do imy=1,pro%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

  end subroutine


  !----------------------------------------------------------------
  subroutine weikmap_report()

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

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

k=0.
    if(n/=1) then
       m = n/2
       
       do i=1,n
          if(i <= m) then
          k(i+shift) = i -1
          else
          k(i+shift) = i - n -1
          end if
       end do
         
       
!       do i=1,m
!          k(i+shift  ) =    i-1
!          k(i+m+shift) = -m+i-1
!       end do
       
    else
       k=0. !! make sure k=0 if n=1
    end if

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







