! 
! Copyright 2008 Shin'ya Murakami <murashin@epa.scitec.kobe-u.ac.jp>
! $Id:$
! 
! 2次元座標を保持
! 
module mod_pointlist
  
  implicit none
  
  ! 2次元座標を表す型
  type Point
     real(8) x
     real(8) y
  end type Point
  
  ! 2次元座標のリンクリストのノード
  type PointListNode
     type(PointListNode), pointer :: n
     type(PointListNode), pointer :: p
     type(Point) data
  end type PointListNode
  
  ! 2次元座標のリンクリスト
  type PointList
     type(PointListNode), pointer :: n
  end type PointList
  
!!!
!!! 演算子定義
!!!

  interface operator ( + )
     module procedure point_plus
  end interface
      
  interface operator ( - )
     module procedure point_minus
  end interface
        
  interface operator ( * )
     module procedure point_multiply_integer
     module procedure point_multiply_real
  end interface
  
  interface operator ( / )
     module procedure point_divide
  end interface

contains
  
  ! 
  ! 初期化
  ! 
  subroutine pl_init(pl)
    
    type(PointList), pointer :: pl
    type(PointListNode), pointer :: node => null()
    
    nullify(pl)
    allocate(pl,node)
    pl%n => node
    nullify(node%n, node%p)
  end subroutine pl_init
  
  ! 
  ! 最初のノードを得る
  ! 
  function pl_get_first(pl) result(retval)
    
    type(PointList), pointer :: pl
    type(PointListNode), pointer :: retval
    
    retval => pl%n%n    
  end function pl_get_first
  
  ! 
  ! 最後のノードを得る.
  ! 
  function pl_get_last(pl) result(retval)
    
    type(PointList), pointer :: pl
    type(PointListNode), pointer :: retval, node => null()
    
    node => pl_each_init(pl)
    retval => node
    do while(pl_each(node))
       retval => node
    end do
  end function pl_get_last
  
  ! 
  ! 指定したノードの次のノードを得る
  ! 
  function pln_get_next(node) result(retval)
    
    type(PointListNode), pointer :: node, retval

    if(associated(node%n))then
       retval => node%n
    else ! last node.
       nullify(retval)
    end if
  end function pln_get_next
  
  ! 
  ! 指定したノードの前のノードを得る
  ! 
  function pln_get_prev(node) result(retval)
    
    type(PointListNode), pointer :: node, retval
    
    if(associated(node%p))then
       retval => node%p
    else
       nullify(retval)
    end if
  end function pln_get_prev
  
  !
  ! 指定したノードの次のノードを削除
  !
  function pln_remove_next(node) result(retval)
    
    type(PointListNode), pointer :: node
    type(Point) retval
    type(PointListNode), pointer :: nnode => null()
    
    if(.not. associated(node%n))then
       ! next node is null
       retval = Point(-999,-999)
    else 
       nnode => node%n
       node%n => node%n%n
       if(associated(nnode%n))then ! 次の次のノードがある
          nnode%n%p => node
       end if
       retval = nnode%data
       deallocate(nnode)
    end if
  end function pln_remove_next
  
  ! 
  ! 指定したノードの次にdataを持つノードを挿入
  ! 
  subroutine pln_add(node, data)
    
    type(PointListNode), pointer :: node
    type(Point), intent(in) :: data
    type(PointListNode), 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 pln_add
  
  !
  ! リストの末尾にdataを持つノードを追加
  !
  subroutine pl_add_last(pl, data)
    
    type(PointList), pointer :: pl
    type(Point), intent(in) :: data
    
    call pln_add(pl_get_last(pl), data)
  end subroutine pl_add_last
  
  !
  ! リストの先頭にdataを持つノードを追加
  !
  subroutine pl_add_first(pl, data)
    
    type(PointList), pointer :: pl
    type(Point), intent(in) :: data
    
    call pln_add(pl%n, data)
  end subroutine pl_add_first
  
  !
  ! 指定したノードのデータを取得
  ! 
  function pln_get_data(node) result(retval)
    
    type(PointListNode), pointer :: node
    type(Point) retval
    
    retval = node%data
  end function pln_get_data
  
  ! 
  ! 指定したノードにデータをセット
  ! 
  subroutine pln_set_data(node, data)
    
    type(PointListNode), pointer :: node
    type(Point), intent(in) :: data
    
    node%data = data
  end subroutine pln_set_data
  
  !
  ! node の次に node_new を挿入
  !
  subroutine pln_add_node(node, node_new)
    
    type(PointListNode), 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 pln_add_node
  
  !
  ! pl の先頭にnodeを挿入
  !
  subroutine pl_add_node_at_first(pl, node)
    
    type(PointList), pointer :: pl
    type(PointListNode), pointer :: node
    
    call pln_add_node(pl%n, node)
  end subroutine pl_add_node_at_first
  
  !
  ! 指定したノードを削除
  !
  function pln_remove(node) result(retval)
    
    type(PointListNode), pointer :: node, retval
    type(PointListNode), 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 pln_remove
  
  ! 
  ! リンクリストを解放する.
  ! 
  subroutine pl_destroy(pl)
    
    type(PointList), pointer :: pl
    type(PointListNode), pointer :: node => null(), nnode => null()
    
    node => pl%n
    deallocate(pl)
    nullify(pl)
    do while(associated(node%n))
       nnode => node%n
       deallocate(node)
       nullify(node)
       node => nnode
    end do
    deallocate(node)
    nullify(node)
  end subroutine pl_destroy
  
  !
  ! while loop用のノード初期化
  !
  function pl_each_init(pl) result(retval)
    
    type(PointList), pointer :: pl
    type(PointListNode), pointer :: retval
    
    retval => pl%n
  end function pl_each_init
  
  !
  ! while loop用
  !
  logical function pl_each(node) result(retval)
    
    type(PointListNode), pointer :: node
    
    node => node%n
    retval = associated(node)
  end function pl_each
  
  ! 
  ! ノードがいくつあるか数える.
  ! 特殊な先頭ノードは数えない.
  ! 
  integer function pl_size(pl) result(retval)
    
    type(PointList), pointer :: pl
    type(PointListNode), pointer :: node => null()
    integer n
    
    n = 0
    node => pl_each_init(pl)
    do while(pl_each(node))
       n = n + 1
    end do
    retval = n
  end function pl_size
  
  ! 
  ! pl1にpl2をくっつける.
  ! 
  subroutine pl_attach_list(pl1, pl2)
    
    type(PointList), pointer :: pl1, pl2
    type(PointListNode), pointer :: n1=> null(), n2=> null()
    
    n1 => pl_get_last(pl1)
    n2 => pl2%n%n
    if(associated(n1))then
       if(associated(n2))then
          ! n1, n2 ともにある
          n1%n => n2
          n2%p => n1
          nullify(pl2%n%n)
       else
          ! n1 はある, n2は null
          ! そのまま
       end if
    else if(associated(n2))then
       ! n1 は null, n2はある
       pl1%n%n => n2
       n2%p => pl1%n
       nullify(pl2%n%n)
    else
       ! n1, n2 はともに null
       ! そのまま
    end if
    call pl_destroy(pl2)
  end subroutine pl_attach_list
  
  ! 
  ! 2点間の距離を求める.
  ! x_is_periodic, y_is_periodic は x, y方向に
  ! 周期境界条件かどうかの論理値
  ! lx, lyはドメインの大きさ
  !
  real(8) function pl_distance(p1, p2, &
       &              x_is_periodic, y_is_periodic, lx, ly) result(retval)
    
    type(Point), intent(in) :: p1, p2
    logical, intent(in) :: x_is_periodic, y_is_periodic
    real(8), intent(in) :: lx, ly
    
    retval = sqrt(distance1d(p1%x, p2%x, x_is_periodic, lx) &
         &       +distance1d(p1%y, p2%y, y_is_periodic, ly))
  end function pl_distance
  
  ! 
  ! 直線上の2点間の二乗距離を返す
  ! 
  real(8) function distance1d(x1, x2, is_periodic, lx) result(retval)
    
    real(8), intent(in) :: x1, x2, lx
    logical, intent(in) :: is_periodic
    
    real(8) tmp1, tmp2, tmp3
    
    tmp1 = (x1 - x2)**2
    tmp2 = (x1 - x2 - lx)**2
    tmp3 = (x1 - x2 + lx)**2
    
    if(is_periodic)then
       if(tmp1 < tmp2)then
          if(tmp1 < tmp3)then
             retval = tmp1
          else
             retval = tmp3
          end if
       else ! tmp1 > tmp2
          if(tmp2 < tmp3)then
             retval = tmp2
          else
             retval = tmp3
          end if
       end if
    else ! 周期境界でない
       retval = tmp1
    end if
  end function distance1d
  
  ! 
  ! PointListのダンプ
  ! 
  subroutine pl_dump(pl)
    
    type(PointList), pointer :: pl
    type(PointListNode), pointer :: node
    integer i 
    
    i = 0
    print *, "dumping point index list...."
    node => pl_each_init(pl)
    do while(pl_each(node))
       i = i + 1
       print *, "pi: ", node%data
    end do
    print *, "end of dumping"
  end subroutine pl_dump
  
  ! 
  ! PointListのdeep copy
  ! 
  function pl_deepcopy(pl) result(retval)
    
    type(PointList), pointer :: pl, retval, pl2 => null()
    type(PointListNode), pointer :: node => null(), node2 => null()
    
    call pl_init(pl2)
    node => pl_each_init(pl)
    node2 => pl2%n
    do while(pl_each(node))
       call pln_add(node2, node%data)
       node2 => node2%n
    end do
    retval => pl2
  end function pl_deepcopy
  
  ! 
  ! plの要素の平均値を計算する.
  ! 桁落ちを避けるため, 二つの要素ごとの和を再帰的にとる.
  ! 
  function pl_mean(pl) result(retval)
    
    type(PointList), pointer :: pl
    type(Point) retval
    
    type(PointList), pointer :: pl2
    type(PointListNode), pointer :: pl_node, pl_node_tmp
    type(Point) p
    
    pl2 => pl_deepcopy(pl)
    
    do while(pl_size(pl2) > 1)
       pl_node => pl_each_init(pl2)
       do while(pl_each(pl_node))
          if(associated(pl_node%n))then
             p = pln_remove_next(pl_node)
             pl_node%data = pl_node%data + p
          else
             pl_node_tmp => pl_node%p
             p = pln_remove_next(pl_node%p)
             pl_node_tmp%data = pl_node_tmp%data + p
             pl_node => pl_node_tmp
          end if
       end do
    end do
    pl_node => pl_get_first(pl2)
    retval = pl_node%data/dble(pl_size(pl))
    call pl_destroy(pl2)
  end function pl_mean

  !
  ! plから2次元配列を作る.
  ! 
  subroutine pl_to_array(pl, array)
    
    type(PointList), pointer :: pl
    real(8), pointer :: array(:,:)

    integer i,n 
    type(PointListNode), pointer :: node
    
    n = pl_size(pl)
    allocate(array(1:2,1:n))
    
    i = 1
    node => pl_each_init(pl)
    do while(pl_each(node))
       array(1,i) = node%data%x
       array(2,i) = node%data%y
       i = i+1
    end do
  end subroutine pl_to_array

!!!
!!! 和差積商
!!!
  
  ! 
  ! Pointどうしの足し算
  ! 
  function point_plus(p1, p2) result(retval)
    
    type(Point), intent(in) :: p1, p2
    type(Point) retval
    
    retval = Point(p1%x + p2%x, p1%y + p2%y)
  end function point_plus
  
  ! 
  ! Pointどうしの引き算
  ! 
  function point_minus(p1, p2) result(retval)
    
    type(Point), intent(in) :: p1, p2
    type(Point) retval
    
    retval = Point(p1%x - p2%x, p1%y - p2%y)
  end function point_minus

  !
  ! Pointとscalarの積 
  !
  function point_multiply_integer(p, a) result(retval)
    
    type(Point), intent(in) :: p
    integer, intent(in) :: a
    type(Point) retval
    
    retval = Point(dble(a)*p%x, dble(a)*p%y)
    
  end function point_multiply_integer
  
  !
  ! Pointとscalarの積 
  !
  function point_multiply_real(p, a) result(retval)
    
    type(Point), intent(in) :: p
    real(8), intent(in) :: a
    type(Point) retval
    
    retval = Point(a*p%x, a*p%y)
  end function point_multiply_real
  
  !
  ! Pointとscalarの商
  !
  function point_divide(p, a) result(retval)
    
    type(Point), intent(in) :: p
    real(8), intent(in) :: a
    type(Point) retval
    
    retval = Point(p%x/a, p%y/a)
  end function point_divide

end module mod_pointlist
