! 
! Copyright 2008-2010 Shin'ya Murakami <murashin@epa.scitec.kobe-u.ac.jp>
! $Id:$
! 
! 渦の統計を取るときにつかう.
! 

module VortexCensus
  
  use d2df
  use mod_realnumberlist
  use mod_pointlist
  use mod_integerlist
  use misc
  
  implicit none
  
contains
  
  ! 
  ! lcondが .true. の格子点列をPatchedArea型にして返す.
  ! 
  function vc_array2pa(lcond, nx, ny) result(retval)
    
    logical, intent(in) :: lcond(0:nx-1,0:ny-1)
    integer, intent(in) :: nx, ny
    type(PatchedArea), pointer :: retval
    
    retval => d2df_array2patchedarea(lcond, nx, ny, .true., .true.)
  end function vc_array2pa
  
  ! 
  ! thr 以下の個数の格子点から構成される塊領域を取り除く
  !
  subroutine vc_pa_small_structure_filter(pa, thr)
    
    type(PatchedArea), pointer :: pa
    integer, intent(in) :: thr
    
    type(PatchedAreaNode), pointer :: pa_node => null(), pa_node_tmp => null()
    
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       if(pil_size(pa_node%i) <= thr)then
          pa_node_tmp => pa_node%p
          call pil_destroy(pan_remove_next(pa_node%p))
          pa_node => pa_node_tmp
       end if
    end do
  end subroutine vc_pa_small_structure_filter

  ! 
  ! 渦領域の数を計算する.
  ! 
  integer function vc_count_vortex(pa) result(retval)
    
    type(PatchedArea), pointer :: pa
    
    retval = pa_size(pa)
  end function vc_count_vortex
  
  !
  ! 渦領域の格子点数を計算する.
  !
  integer function vc_count_gridpoints(pa) result(retval)
    
    type(PatchedArea), pointer :: pa
    type(PatchedAreaNode), pointer :: pa_node => null()
    integer vgp
    
    vgp = 0
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       vgp = vgp + pil_size(pa_node%i)
    end do
    retval = vgp
  end function vc_count_gridpoints
  
  !
  ! 渦の平均半径 R を求める.
  ! 真円を仮定して, R = sqrt(渦の面積/pi) とする.
  !
  real(8) function vc_mean_radii(pa, nx, ny, dkx, dky) result(retval)
    
    type(PatchedArea), pointer :: pa
    integer, intent(in) :: nx, ny
    real(8), intent(in) :: dkx, dky
    type(PatchedAreaNode), pointer :: pa_node => null()
    real(8) r, dpi, pi, a
    
    pi = acos(-1.0d0)
    dpi = 2.0d0*pi
    a = dpi**2/(dble(nx*ny)*dkx*dky) ! 格子点数 * a = 格子点が占める面積
    
    r = 0.0d0
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       r = r + dsqrt(a * pil_size(pa_node%i)/pi)
    end do
    retval = r/dble(pa_size(pa))
  end function vc_mean_radii
  
  ! 
  ! それぞれの渦の循環の大きさを渦ごとに平均した値を返す
  ! 
  real(8) function vc_mean_circulation(pa, vor, nx, ny, dkx, dky) result(retval)
    
    type(PatchedArea), pointer :: pa
    real(8), intent(in) :: vor(0:nx-1, 0:ny-1), dkx, dky
    integer, intent(in) :: nx, ny
    
    type(PatchedAreaNode), pointer :: pa_node => null()
    type(PointIndexListNode), pointer :: pil_node => null()
    real(8) circ, circ_each
    
    circ = 0.0d0
    circ_each = 0.0d0
    
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       circ_each = 0.0d0
       pil_node => pil_each_init(pa_node%i)
       do while(pil_each(pil_node))
          circ_each = circ_each + vor(pil_node%data%ix, pil_node%data%iy)
       end do
       circ = circ + abs(circ_each)
    end do
    circ = circ/dble(pa_size(pa))
    retval = circ/dble(nx*ny)*area(dkx, dky)
  end function vc_mean_circulation
  
  ! 
  ! 渦の保有エンストロフィーの平均値を返す
  ! 
  real(8) function vc_mean_enstrophy(pa, vor, nx, ny, dkx, dky) result(retval)
    
    type(PatchedArea), pointer :: pa
    real(8), intent(in) :: vor(0:nx-1, 0:ny-1), dkx, dky
    integer, intent(in) :: nx, ny
    
    type(PatchedAreaNode), pointer :: pa_node => null()
    type(PointIndexListNode), pointer :: pil_node => null()
    real(8) enst
    
    enst = 0.0d0
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       pil_node => pil_each_init(pa_node%i)
       do while(pil_each(pil_node))
          enst = enst + vor(pil_node%data%ix, pil_node%data%iy)**2
       end do
    end do
    enst = 0.5d0*enst !渦領域中の全格子点におけるエンストロフィーの和
    enst = enst/dble(nx*ny)*area(dkx, dky) ! 全渦領域中のエンストロフィーの和
    retval = enst/dble(pa_size(pa)) ! 渦領域１個あたりのエンストロフィー
  end function vc_mean_enstrophy
  
  ! 
  ! 渦の保有エネルギーの平均値を返す.
  ! ただし, エネルギーは -\zeta \psiの空間積分とする.
  ! 
  real(8) function vc_mean_energy(pa, vor, str, nx, ny, dkx, dky) result(retval)
    
    type(PatchedArea), pointer :: pa
    real(8), intent(in) :: vor(0:nx-1, 0:ny-1), str(0:nx-1, 0:ny-1), dkx, dky
    integer, intent(in) :: nx, ny
    
    type(PatchedAreaNode), pointer :: pa_node => null()
    type(PointIndexListNode), pointer :: pil_node => null()
    real(8) energy
    
    energy = 0.0d0
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       pil_node => pil_each_init(pa_node%i)
       do while(pil_each(pil_node))
          energy = energy  - vor(pil_node%data%ix, pil_node%data%iy) &
               &            *str(pil_node%data%ix, pil_node%data%iy)
       end do
    end do
    energy = energy/dble(pa_size(pa))
    retval = 0.5d0*energy/dble(nx*ny)*area(dkx, dky)
  end function vc_mean_energy
  
  ! 
  ! 渦の保有パリンストロフィーの平均値を返す
  ! 
  real(8) function vc_mean_palinstrophy(pa, palin, nx, ny, dkx, dky) result(retval)
    
    type(PatchedArea), pointer :: pa
    real(8), intent(in) :: palin(0:nx-1, 0:ny-1), dkx, dky
    integer, intent(in) :: nx, ny
    
    type(PatchedAreaNode), pointer :: pa_node => null()
    type(PointIndexListNode), pointer :: pil_node => null()
    real(8) palin_avg
    
    palin_avg = 0.0d0
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       pil_node => pil_each_init(pa_node%i)
       do while(pil_each(pil_node))
          palin_avg = palin_avg + palin(pil_node%data%ix, pil_node%data%iy)
       end do
    end do
    palin_avg = 0.5d0*palin_avg !渦領域中の全格子点におけるパリンストロフィーの和
    palin_avg = palin_avg/dble(nx*ny)*area(dkx, dky) !全渦領域中のパリンストロフィーの和
    retval = palin_avg/dble(pa_size(pa)) !渦領域1個あたりのパリンストロフィー
  end function vc_mean_palinstrophy
  
  ! 
  ! 渦の中心位置のリンクリストを生成.
  ! 単に位置の平均値を求める.
  ! 渦度重み付きの位置平均値にした方が良いかもしれない.
  !
  ! 境界を跨ぐ渦がある場合には, 渦の大きさが領域の4分の1未満の
  ! 大きさの場合のみうまく計算できる.
  ! 
  function vc_pa2positions_of_center(pa, x, y, nx, ny, lx, ly) result(retval)
    
    type(PatchedArea), pointer :: pa
    integer, intent(in) :: nx, ny
    real(8), intent(in) :: x(0:nx-1), y(0:ny-1), lx, ly
    type(PointList), pointer :: retval
    
    type(PointList), pointer :: pl => null()
    type(PatchedAreaNode), pointer :: pa_node => null()
    type(PointIndexListNode), pointer :: pil_node => null()
    type(IntegerList), pointer :: il => null()
    type(IntegerListNode), pointer :: il_node => null()
    
    type(Point) p
    real(8) x0, y0
    logical is_over_x, is_over_y
    
    ! 境界をまたぐ渦領域に対応するノードのid のリストを作成
    il => vc_bpil2il(pa%b)
    
    call pl_init(pl)
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       p = Point(0.0d0, 0.0d0)
       if(il_include_data(il, pa_node%id, il_node))then
          ! 境界をまたぐノードの場合は特別あつかい.
          call bpil_is_over_direction(pa%b, pa_node%id, is_over_x, is_over_y)
          pil_node => pil_each_init(pa_node%i)
          do while(pil_each(pil_node))
             x0 = x(pil_node%data%ix)
             y0 = y(pil_node%data%iy)
             if(is_over_x)then
                ! ly/2 < y0 なる位置 y0 は, y0 - Ly で勘定する
                if(0.5d0*ly < y0)then
                   p = p + Point(0.0d0, y0-ly)
                else
                   p = p + Point(0.0d0, y0)
                end if
             else
                p = p + Point(0.0d0, y0)
             end if
             if(is_over_y)then
                ! lx/2 < x0 なる位置 x0 は, x0 - Lx で勘定する
                if(0.5d0*lx < x0)then
                   p = p + Point(x0-lx, 0.0d0)
                else
                   p = p + Point(x0, 0.0d0)
                end if
             else
                p = p + Point(x0, 0.0d0)
             end if
          end do
          p = p/dble(pil_size(pa_node%i)) ! 位置の平均値
          ! [0, lx)\times[0, ly) から飛び出ている場合は中に入れる.
          if(p%x < 0.0d0)then
             p%x = p%x + lx
          end if
          if(p%y < 0.0d0)then
             p%y = p%y + ly
          end if
       else
          pil_node => pil_each_init(pa_node%i)
          do while(pil_each(pil_node))
             p = p + Point(x(pil_node%data%ix), y(pil_node%data%iy))
          end do
          p = p/dble(pil_size(pa_node%i)) ! 位置の平均値
       end if
       call pl_add_first(pl, p)
    end do
    call il_destroy(il)
    retval => pl
  end function vc_pa2positions_of_center
  
  ! 
  ! 境界情報から, 境界格子点を含む渦領域のidのリストを得る
  ! 
  function vc_bpil2il(bpil) result(retval)
    
    type(BoundaryPointIndexList), pointer :: bpil
    type(IntegerList), pointer :: il => null(), retval
    
    type(BoundaryPointIndexListNode), pointer :: bpil_node
    
    call il_init(il)
    bpil_node => bpil_each_init(bpil)
    do while(bpil_each(bpil_node))
       call il_add_first(il, bpil_node%data%area_id)
    end do
    call il_sort(il)
    call il_uniq(il)
    retval => il
  end function vc_bpil2il
  
  !
  ! 渦間の平均距離を求める.
  !
  real(8) function vc_mean_distance_of_vortices(pl, lx, ly) result(retval)
    
    type(PointList), pointer :: pl
    real(8), intent(in) :: lx, ly
    type(PointListNode), pointer :: pl_node1 => null(), pl_node2 => null()
    ! integer(8) i
    integer(8) n
    real(8) r
    
    ! i = 0
    n = pl_size(pl)
    r = 0.0d0
    pl_node1 => pl_each_init(pl)
    do while(pl_each(pl_node1))
       if(associated(pl_node1%n))then
          pl_node2 => pl_node1
          do while(pl_each(pl_node2))
             r = r + pl_distance(pl_node1%data, pl_node2%data, &
                  &                     .true., .true., lx, ly)
             ! i = i + 1
          end do
       end if
    end do
    ! print *, "number of comparing", (n*(n-1))/2, i
    retval = r/dble((n*(n-1))/2)
  end function vc_mean_distance_of_vortices
  
  !
  ! 渦間の平均距離を求める, その2. 桁落ちを意識した平均をとる.
  ! ! 一般に, 結果はほとんど変わらない....
  !
  real(8) function vc_mean_distance_of_vortices2(pl, lx, ly) result(retval)
    
    type(PointList), pointer :: pl
    real(8), intent(in) :: lx, ly
    type(RealNumberList), pointer :: rnl => null()
    type(PointListNode), pointer :: pl_node1 => null(), pl_node2 => null()
    real(8) r
    
    call rnl_init(rnl)
    r = 0.0d0
    pl_node1 => pl_each_init(pl)
    do while(pl_each(pl_node1))
       if(associated(pl_node1%n))then
          pl_node2 => pl_node1
          do while(pl_each(pl_node2))
             r = pl_distance(pl_node1%data, pl_node2%data, &
                  &                 .true., .true., lx, ly)
             call rnl_add_first(rnl, r)
          end do
       end if
    end do
    ! print *, "number of comparing", (n*(n-1))/2, i
    retval = rnl_mean(rnl)
    call rnl_destroy(rnl)
  end function vc_mean_distance_of_vortices2
  
  !
  ! 各渦と最も近い渦との間の平均距離を求める.
  !
  real(8) function vc_mean_distance_of_nearest_vortices(pl, lx, ly) result(retval)
    
    type(PointList), pointer :: pl
    real(8), intent(in) :: lx, ly
    type(PointListNode), pointer :: pl_node1 => null(), pl_node2 => null()
    type(PointListNode), pointer :: pl_node_tmp => null()
    integer(8) n
    type(Point) p
    real(8) r_sum, r, tmp
    
    n = pl_size(pl)
    r_sum = 0.0d0
    pl_node1 => pl_each_init(pl)
    do while(pl_each(pl_node1))
       pl_node_tmp => pl_node1%p
       p = pln_remove_next(pl_node_tmp)
       r = 1.0d10 ! 十分大きなrからスタート
       pl_node2 => pl_each_init(pl)
       do while(pl_each(pl_node2))
          tmp = pl_distance(p, pl_node2%data,.true.,.true.,lx,ly)
          if(tmp < r)then
             r = tmp
          end if
       end do
       call pln_add(pl_node_tmp, p)
       r_sum = r_sum + r
    end do
    retval = r_sum/dble(n)
  end function vc_mean_distance_of_nearest_vortices

!   ! 
!   ! 各々の渦領域での関数の極値を求める.
!   ! 
!   function vc_function_extremum_each_region(pa, f, nx, ny) result(retval)
    
!     type(PatchedArea), pointer :: pa
!     real(8), intent(in) :: f(0:nx-1, 0:ny-1)
!     integer, intent(in) :: nx, ny
!     type(PointList), pointer :: retval
!     type(PointList), pointer :: pl => null()
!     type(PatchedAreaNode), pointer :: pa_node => null()
!     type(PointIndexListNode), pointer :: pil_node => null()
    
!     type(Point) p
    
!     call pl_init(pl)
!     pa_node => pa_each_init(pa)
!     do while(pa_each(pa_node))
!        p = Point(0.0d0, 0.0d0)
!        pil_node => pil_each_init(pa_node%i)
!        do while(pil_each(pil_node))
!           p = p + Point(x(pil_node%data%ix), y(pil_node%data%iy))
!        end do
!        p = p/dble(pil_size(pa_node%i)) ! 位置の平均値
!        call pl_add_first(pl, p)
!     end do
!     retval => pl
!   end function vc_function_extremum_each_region
  
  ! 
  ! pl で指定された位置における関数値 f を求め, 平均する.
  ! 関数値 f はスペクトルから直接求める.
  ! 
  function vc_function_mean_at_points(pl, f, nx, ny, dkx, dky, ky) result(retval)
    
    type(PointList), pointer :: pl
    complex(8), intent(in) :: f(0:nx/2, 0:ny-1)
    integer, intent(in) :: nx, ny, ky(0:ny-1)
    real(8), intent(in) :: dkx, dky
    real(8) retval, fsum
    type(PointListNode), pointer :: pl_node
    
    fsum = 0.0d0
    
    pl_node => pl_each_init(pl)
    do while(pl_each(pl_node))
       fsum = fsum + funcval2d(f, pl_node%data%x, pl_node%data%y, nx, ny, dkx, dky, ky)
    end do
    retval = fsum/dble(pl_size(pl))
    
  end function vc_function_mean_at_points
  
  ! 
  ! pl で指定された位置における関数値 f を求め, 平均する.
  ! 関数値 f は実空間場から線形補間する.
  ! 
  real(8) function vc_function_mean_at_points2(pl, f, x, y, nx, ny, dkx, dky) result(retval)
    
    type(PointList), pointer :: pl
    real(8), intent(in) :: f(0:nx/2, 0:ny-1)
    integer, intent(in) :: nx, ny
    real(8), intent(in) :: dkx, dky, x(0:nx-1), y(0:ny-1)
    real(8) fsum, x0, y0
    type(PointListNode), pointer :: pl_node
    
    fsum = 0.0d0
    
    pl_node => pl_each_init(pl)
    do while(pl_each(pl_node))
       x0 = pl_node%data%x
       y0 = pl_node%data%y
       fsum = fsum + InterpolateLinear2d(f, x0, y0, x, y, nx, ny, dkx, dky)
    end do
    retval = fsum/dble(pl_size(pl))
    
  end function vc_function_mean_at_points2
  
  ! 
  ! pl で指定された位置における関数値 f の絶対値を求め, 平均する.
  ! 関数値 f は実空間場から線形補間する.
  ! 
  real(8) function vc_abs_function_mean_at_points2(pl, f, x, y, nx, ny, dkx, dky) result(retval)
    
    type(PointList), pointer :: pl
    real(8), intent(in) :: f(0:nx/2, 0:ny-1), dkx, dky, x(0:nx-1), y(0:ny-1)
    integer, intent(in) :: nx, ny
    real(8) fsum, x0, y0
    type(PointListNode), pointer :: pl_node
    
    fsum = 0.0d0
    
    pl_node => pl_each_init(pl)
    do while(pl_each(pl_node))
       x0 = pl_node%data%x
       y0 = pl_node%data%y
       fsum = fsum + abs(InterpolateLinear2d(f, x0, y0, x, y, nx, ny, dkx, dky))
    end do
    retval = fsum/dble(pl_size(pl))
    
  end function vc_abs_function_mean_at_points2
  
  !
  ! 渦領域が楕円形であると仮定して, 偏平率を計算する.
  ! 
  ! 1. 領域の境界の格子点列を得たとする.
  ! 2. 領域の境界の格子点列から一つずつ選び, 
  !    全ての境界の格子点から, 最も距離が遠いものを選ぶ.
  !    格子点の組を覚えておく.
  ! 3. これを全ての境界格子点に対して行う.
  !    最も距離が遠かった組を長軸とする.
  !    (これは必ずしも中心を通らない.)
  ! 4. 長軸から短軸を求める.
  !    面倒なので, 長軸とx軸の角を求めて, 
  !    長軸がx軸に合うように座標変換する.
  !    スペクトルでやるのは面倒なので, 線形内挿でやる.
  !    y軸の長さが短軸の長さになる.
  ! 5. 偏平率を求める.
  ! 
  ! 問題: 楕円形が仮定されているが, 
  !       そもそも離散的な格子上で理想的な楕円は表現できないわけで, 
  !       理想的な楕円を格子状に離散化したときにどのくらいの誤差が出るか.
  !       また, 楕円じゃない場合でも一応計算できてしまうが...
  !       
  
  ! 
  ! 領域の格子点列から, 領域の境界格子点列を得る.
  ! 
  ! 
  ! まず, 初期値として, y軸,x軸ともに最も小さい点を探す.
  ! ! と思ったけど, 周期境界だからこれではだめだ.
  ! ! 周期境界を跨いでいる場合にだけ, 最大値をとるようにしたりとか細工するか?
  ! 
  ! 境界を跨いでない場合は, 初期の点から, 
  ! 上へ, 左へ, 下へ, 右へ, の順で探索.
  ! 一周したら失敗.
  ! 上行って左がなかったら, そこは境界として登録.
  ! 見付からなかったらその度にその点を境界として登録.
  ! これだとジグザグに連結された領域は得られない.....
  
  
  ! 
  ! 個々の渦の数密度を, 渦面積をインデックスとするリストとして求める.
  ! 
  subroutine vc_pa2number_density_with_area(pa, pdf, n, xmin, dx, dA)
    
    type(PatchedArea), pointer :: pa
    integer, intent(in)  :: n
    real(8), intent(in)  :: xmin, dx, dA
    real(8), intent(out) :: pdf(0:n-1)
    
    type(PatchedAreaNode), pointer :: pa_node => null()
    real(8) logsize
    integer i
    
    pdf = 0.0d0
    
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       logsize = log(dble(pil_size(pa_node%i))*dA)
       i = coord_to_index(logsize, xmin, dx)
       pdf(i) = pdf(i) + 1
    end do
  end subroutine vc_pa2number_density_with_area
  
  ! 
  ! 渦面積 A の最大値と最小値を計算する
  ! 
  subroutine vc_count_maxmin_vortexarea(pa, xmax, xmin)
    
    type(PatchedArea), pointer :: pa
    real(8), intent(out)  :: xmin, xmax
    
    type(PatchedAreaNode), pointer :: pa_node => null()
    real(8) size
    
    ! patchedarea でループ
    ! 個々のpatchedareaについて, 格子点数を数える
    ! 分類して配列に格納. 
    
    xmax = 1.0d-100
    xmin = 1.0d100
    
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       size = pil_size(pa_node%i)
       if(xmax < size)then
          xmax = size
       end if
       if(xmin > size)then
          xmin = size
       end if
    end do
  end subroutine vc_count_maxmin_vortexarea

  ! 
  ! 最大の格子点数で構成される渦のみを取り出す.
  !
  subroutine vc_pa_choose_maximum_area(pa)
    
    type(PatchedArea), pointer :: pa

    type(PatchedAreaNode), pointer :: pa_node => null(), pa_node_tmp => null()
    real(8) xmax, xmin
    
    call vc_count_maxmin_vortexarea(pa, xmax, xmin)
    
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       if(pil_size(pa_node%i) < xmax)then
          pa_node_tmp => pa_node%p
          call pil_destroy(pan_remove_next(pa_node%p))
          pa_node => pa_node_tmp
       end if
    end do
  end subroutine vc_pa_choose_maximum_area

  ! 
  ! 重心が(x0,y0)に最も近い渦のみを取り出す.
  !
  subroutine vc_pa_choose_nearest(pa, x0, y0, nx, ny, lx, ly)
    
    type(PatchedArea), pointer :: pa
    integer, intent(in) :: nx, ny
    real(8), intent(in) :: x0, y0, lx, ly
    
    type(PatchedAreaNode), pointer :: pa_node => null(), pa_node_tmp => null()
    real(8) xc, yc, xmin, d2
    integer i, xmin_index

    xmin = 1.0d3
    xmin_index = -1
    
    ! 重心が(x0,y0)に最も近い要素の存在するインデックスを調べる
    i = 0
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       i = i + 1
       call pil_centroid(pa_node%i, xc, yc) ! インデックスの重心が帰って来る
       xc = xc/dble(nx)*lx
       yc = yc/dble(ny)*ly
       d2 = (xc-x0)**2 +(yc-y0)**2
       if(d2 < xmin)then
          xmin = d2
          xmin_index = i
       end if
    end do
    ! xmin_index番目以外を削除
    i = 0
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       i = i + 1
       if( xmin_index /= i)then
          pa_node_tmp => pa_node%p
          call pil_destroy(pan_remove_next(pa_node%p))
          pa_node => pa_node_tmp
       end if
    end do
  end subroutine vc_pa_choose_nearest

  ! 
  ! 渦領域の物理量の自乗平均値の最大の領域を取り出す.
  ! 境界を跨ぐ場合は考えてない.
  !
  subroutine vc_pa_choose_maximum(pa, phys, nx, ny)
    
    type(PatchedArea), pointer :: pa
    real(8), intent(in) :: phys(0:nx-1, 0:ny-1)
    integer, intent(in) :: nx, ny
    
    type(PatchedAreaNode), pointer :: pa_node => null(), pa_node_tmp => null()
    type(PointIndexListNode), pointer :: pil_node => null()
    real(8) phys_sum_max, phys_sum
    integer i, max_index, n
    
    phys_sum_max = 1.0d-15
    max_index = -1
    
    ! 物理量の自乗平均値が最大の領域になるインデックスを探す
    i = 0
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       i = i + 1
       phys_sum = 0.0d0
       n = pil_size(pa_node%i)
       pil_node => pil_each_init(pa_node%i)
       do while(pil_each(pil_node))
          phys_sum = phys_sum + phys(pil_node%data%ix, pil_node%data%iy)**2
       end do
       phys_sum = phys_sum/dble(n)
       if(phys_sum_max < phys_sum)then
          phys_sum_max = phys_sum
          max_index = i
       end if
    end do
    ! max_index 番目以外を削除
    i = 0
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       i = i + 1
       if( max_index /= i)then
          pa_node_tmp => pa_node%p
          call pil_destroy(pan_remove_next(pa_node%p))
          pa_node => pa_node_tmp
       end if
    end do
  end subroutine vc_pa_choose_maximum

end module VortexCensus
