! 
! Copyright 2008 Shin'ya Murakami <murashin@epa.scitec.kobe-u.ac.jp>
! $Id:$
! 
! 2次元座標のインデックスを保持
! 
module mod_pointindexlist
  
  implicit none
  
  ! 2次元座標のインデックスを表す型
  type PointIndex
     integer ix
     integer iy
  end type PointIndex
  
  ! 2次元座標インデックスのリンクリストのノード
  type PointIndexListNode
     type(PointIndexListNode), pointer :: n
     type(PointIndexListNode), pointer :: p
     type(PointIndex) data
  end type PointIndexListNode

  ! 2次元座標インデックスのリンクリスト
  type PointIndexList
     type(PointIndexListNode), pointer :: n
  end type PointIndexList

  ! 
  ! 演算子定義
  ! 
  interface operator ( == )
     module procedure pi_is_equal
  end interface

contains

  ! 
  ! 初期化
  ! 
  subroutine pil_init(pil)
    
    type(PointIndexList), pointer :: pil
    type(PointIndexListNode), pointer :: node => null()
    
    nullify(pil)
    allocate(pil,node)
    pil%n => node
    nullify(node%n, node%p)
  end subroutine pil_init
  
  ! 
  ! 最初のノードを得る
  ! 
  function pil_get_first(pil) result(retval)
    
    type(PointIndexList), pointer :: pil
    type(PointIndexListNode), pointer :: retval
    
    retval => pil%n%n
    
  end function pil_get_first

  ! 
  ! 最後のノードを得る.
  ! 
  function pil_get_last(pil) result(retval)
    
    type(PointIndexList), pointer :: pil
    type(PointIndexListNode), pointer :: retval, node => null()

    node => pil_each_init(pil)
    retval => node
    do while(pil_each(node))
       retval => node
    end do
  end function pil_get_last
  
  ! 
  ! 指定したノードの次のノードを得る
  ! 
  function piln_get_next(node) result(retval)
    
    type(PointIndexListNode), pointer :: node, retval

    if(associated(node%n))then
       retval => node%n
    else ! last node.
       nullify(retval)
    end if
  end function piln_get_next
  
  ! 
  ! 指定したノードの前のノードを得る
  ! 
  function piln_get_prev(node) result(retval)
    
    type(PointIndexListNode), pointer :: node, retval
    
    if(associated(node%p))then
       retval => node%p
    else
       nullify(retval)
    end if
  end function piln_get_prev
  
  !
  ! 指定したノードの次のノードを削除
  !
  function piln_remove_next(node) result(retval)
    
    type(PointIndexListNode), pointer :: node
    type(PointIndex) retval
    type(PointIndexListNode), pointer :: nnode => null()
    
    if(.not. associated(node%n))then
       ! next node is null
       retval = PointIndex(-999,-999)
    else
       nnode => node%n
       node%n => node%n%n
       if(associated(node%n))then ! 次の次のノードがある
          node%n%p => node
       end if
       retval = nnode%data
       deallocate(nnode)
    end if
  end function piln_remove_next
  
  !
  ! 指定したノードの次にdataを持つノードを挿入
  !
  subroutine piln_add(node, data)
    
    type(PointIndexListNode), pointer  :: node
    type(PointIndex), intent(in) :: data
    
    type(PointIndexListNode), pointer :: node_new => null()
    
    allocate(node_new)
    node_new%data = data
    node_new%n => node%n
    node_new%p => node
    if(associated(node%n))then
       node%n%p => node_new
    end if
    node%n => node_new
  end subroutine piln_add
  
  !
  ! ノードの末尾にdataを持つノードを追加
  !
  subroutine pil_add_last(pil, data)
    
    type(PointIndexList), pointer  :: pil
    type(PointIndex), intent(in) :: data
    
    call piln_add(pil_get_last(pil), data)
  end subroutine pil_add_last
  
  !
  ! リストの先頭にdataを持つノードを追加
  !
  subroutine pil_add_first(pil, data)
    
    type(PointIndexList), pointer  :: pil
    type(PointIndex), intent(in) :: data
    
    call piln_add(pil%n, data)
  end subroutine pil_add_first
  
  !
  ! 指定したノードのデータを取得
  ! 
  function piln_get_data(node) result(retval)
    
    type(PointIndexListNode), pointer   :: node
    type(PointIndex) retval
    
    retval = node%data
  end function piln_get_data
  
  ! 
  ! 指定したノードにデータをセット
  ! 
  subroutine piln_set_data(node, data)
    
    type(PointIndexListNode), pointer   :: node
    type(PointIndex), intent(in) :: data
    
    node%data = data
  end subroutine piln_set_data
  
  ! 
  ! node の次にnode_newを挿入
  ! 
  subroutine piln_add_node(node, node_new)
    
    type(PointIndexListNode), pointer :: node, node_new
    
    node_new%n => node%n
    node_new%p => node
    node%n => node_new
    if(associated(node_new%n))then
       node_new%n%p => node_new
    end if
  end subroutine piln_add_node
  
  !
  ! pilの先頭にnodeを挿入
  !
  subroutine pil_add_node_at_first(pil, node)
    
    type(PointIndexList), pointer :: pil
    type(PointIndexListNode), pointer :: node
    
    call piln_add_node(pil%n, node)
  end subroutine pil_add_node_at_first
  
  !
  ! 指定したノードを削除
  !
  function piln_remove(node) result(retval)
    
    type(PointIndexListNode), pointer :: node, retval
    type(PointIndexListNode), pointer :: nnode => null(), pnode => null()
   
    if(.not. associated(node%n))then
       pnode => node%p
       nullify(pnode%n, node%p)
       retval => node
    else ! others
       pnode => node%p
       nnode => node%n
       pnode%n => nnode
       nnode%p => pnode
       nullify(node%p, node%n)
       retval => node
    end if
  end function piln_remove
  
  ! 
  ! リンクリストを解放する.
  ! 
  subroutine pil_destroy(pil)
    
    type(PointIndexList), pointer :: pil
    type(PointIndexListNode), pointer :: node => null(), nnode => null()
    
    node => pil%n
    deallocate(pil)
    nullify(pil)
    do while(associated(node%n))
       nnode => node%n
       deallocate(node)
       nullify(node)
       node => nnode
    end do
    deallocate(node)
    nullify(node)
  end subroutine pil_destroy
  
  !
  ! while loop用のノード初期化
  !
  function pil_each_init(pil) result(retval)
    
    type(PointIndexList), pointer :: pil
    type(PointIndexListNode), pointer :: retval
    
    retval => pil%n
  end function pil_each_init
  
  !
  ! while loop用
  !
  logical function pil_each(node) result(retval)
    
    type(PointIndexListNode), pointer :: node
    
    node => node%n
    retval = associated(node)
  end function pil_each
  
  !
  ! 要素がいくつあるか数える.
  ! 
  function pil_size(pil) result(retval)
    
    type(PointIndexList), pointer :: pil
    
    type(PointIndexListNode), pointer :: node => null()
    integer n, retval
    
    n = 0
    node => pil_each_init(pil)
    do while(pil_each(node))
       n = n+1
    end do
    retval = n
  end function pil_size
  
  !
  ! 要素がゼロかどうかを返す
  ! 
  function pil_is_empty(pil) result(retval)
    
    type(PointIndexList), pointer :: pil
    logical retval
    
    retval = associated(pil%n%n)
  end function pil_is_empty
  
  !
  ! pil1にpil2をくっつける.
  !
  subroutine pil_attach_list(pil1, pil2)
    
    type(PointIndexList), pointer :: pil1, pil2
    type(PointIndexListNode), pointer :: n1=> null(), n2=> null()
    
    n1 => pil_get_last(pil1)
    n2 => pil2%n%n
    if(associated(n1))then
       if(associated(n2))then
          ! n1, n2 ともにある
          n1%n => n2
          n2%p => n1
          nullify(pil2%n%n)
       else
          ! n1 はある, n2は null
          ! そのまま
       end if
    else if(associated(n2))then
       ! n1 は null, n2はある
       pil1%n%n => n2
       n2%p => pil1%n
       nullify(pil2%n%n)
    else
       ! n1, n2 はともに null
       ! そのまま
    end if
    call pil_destroy(pil2)
  end subroutine pil_attach_list

  ! 
  ! ソート
  ! 
  ! 
  ! TupleListの要素を昇順にソート
  ! 
  subroutine pil_sort(pil)
    
    type(PointIndexList), pointer :: pil
    type(PointIndexListNode), pointer :: n1 => null(), n2 => null()
    
    n1 => pil_each_init(pil)
    do while(pil_each(n1))
       n2 => n1
       do while(pil_each(n2))
          call piln_swap(n1, n2)
       end do
    end do
  end subroutine pil_sort

  ! 
  ! 指定したノードのデータ間を比較して, 昇順になるようにスワップ
  ! 
  subroutine piln_swap(n1, n2)
    
    type(PointIndexListNode), pointer :: n1, n2
    type(PointIndex) tmp
    
    if(n1%data%ix == n2%data%ix)then
       if(n1%data%iy > n2%data%iy)then
          tmp = n2%data
          n2%data = n1%data
          n1%data = tmp
       end if
    else if(n1%data%ix > n2%data%ix)then
       tmp = n2%data
       n2%data = n1%data
       n1%data = tmp
    end if
  end subroutine piln_swap

  ! 
  ! 各隣接要素でデータが同じ場合に削除する.
  ! 
  subroutine pil_uniq(pil)
    
    type(PointIndexList), pointer :: pil
    type(PointIndexListNode), pointer :: n1 =>null(), n2 =>null()
    type(PointIndex) tmp
    
    n1 => pil%n%n
    if(associated(n1))then
       n2 => n1%n
       do while(associated(n2))
          if(n1%data == n2%data)then
             tmp = piln_remove_next(n1)
             n2 => n1%n
          else
             n1 => n1%n
             n2 => n2%n
          end if
       end do
    end if
  end subroutine pil_uniq
  
  
!!!
!!! 配列からリストへの書きだし, リストから配列への書きだし.
!!!
  
  ! 
  ! 配列からリストへの書きだし
  ! 
  function pil_array2list(lary, nx, ny) result(retval)
    
    logical, intent(in) :: lary(0:nx-1, 0:ny-1)
    integer, intent(in) :: nx, ny 
    type(PointIndexList), pointer :: retval
    
    retval => pil_parray2list(lary, nx, ny, 0, nx-1, 0, ny-1)
  end function pil_array2list
  
  ! 
  ! 配列の一部からリストへの書きだし
  ! 
  function pil_parray2list(lary, nx, ny, ix0, ix1, iy0, iy1) result(retval)
    
    logical, intent(in) :: lary(0:nx-1, 0:ny-1)
    integer, intent(in) :: nx, ny, ix0, ix1, iy0, iy1
    type(PointIndexList), pointer :: pil => null(), retval
    integer i,j
    type(PointIndex) tmp
    
    call pil_init(pil)
    do j=iy0,iy1
       do i=ix0,ix1
          if(lary(i,j))then
             tmp%ix = i
             tmp%iy = j
             call pil_add_first(pil, tmp)
          end if
       end do
    end do
    retval => pil
  end function pil_parray2list
  
  ! 
  ! リストから配列への書き戻し. 
  ! 配列は適切に初期化されているとする.
  ! 
  subroutine pil_list2array(lary, nx, ny, pil)
    
    logical, intent(inout) :: lary(0:nx-1, 0:ny-1)
    integer, intent(in) :: nx, ny
    type(PointIndexList), pointer :: pil
    type(PointIndexListNode), pointer :: node
    type(PointIndex) tmp
    
    node => pil_each_init(pil)
    do while(pil_each(node))
       tmp = piln_get_data(node)
       lary(tmp%ix, tmp%iy) = .true.
    end do
  end subroutine pil_list2array
  
  ! 
  ! pi%i >= 0, pi%j >=0 なpi1, pi2が隣接しているか比較.
  ! 
  ! !差をとって, 絶対値の和を調べて, 1なら隣接.
  ! 
  logical function pi_is_adjacent(pi1, pi2) result(retval)
    
    type(PointIndex), intent(in) :: pi1, pi2
    
    if((abs(pi1%ix-pi2%ix) + abs(pi1%iy-pi2%iy)) == 1)then
       retval = .true.
    else
       retval = .false.
    end if
  end function pi_is_adjacent
  
  !
  ! PointIndexListのダンプ
  !
  subroutine pil_dump(pil)
    
    type(PointIndexList), pointer :: pil
    type(PointIndexListNode), pointer :: node
    integer i 
    
    i = 0
    print *, "dumping point index list...."
    node => pil_each_init(pil)
    do while(pil_each(node))
       i = i + 1
       print *, "pi: ", node%data
    end do
    print *, "end of dumping"
  end subroutine pil_dump

  !
  ! 重心位置を計算する.
  !
  subroutine pil_centroid(pil, xc, yc)
    
    type(PointIndexList), pointer :: pil
    real(8), intent(out) :: xc, yc

    type(PointIndexListNode), pointer :: node
    integer n

    xc = 0.0d0
    yc = 0.0d0

    n = pil_size(pil)

    node => pil_each_init(pil)
    do while(pil_each(node))
       xc = xc + node%data%ix
       yc = yc + node%data%iy
    end do
    xc = xc/dble(n)
    yc = yc/dble(n)
  end subroutine pil_centroid
  
  ! 
  ! PointIndexの比較.
  ! 
  function pi_is_equal(a,b) result(retval)
    
    type(PointIndex), intent(in) :: a,b
    logical retval
    
    if(a%ix == b%ix .and. a%iy == b%iy)then
       retval = .true.
    else
       retval = .false.
    end if
  end function pi_is_equal

end module mod_pointindexlist
