!> \~japanese
!! \brief 正二十面体格子点およびコントロールボリュームの座標データ等の計算を行う手続きを提供するモジュール.
!!
!! - STDGrid_Builder モジュールについて
!!   - 実際の格子点座標やコントロールボリュームの頂点座標などを求める計算は, 本モジュールが提供する手続きを通して行われ,
!!     その結果は構造体 IcGrid2D_FVM の対応する成分にそれぞれ格納される.
!!   - このモジュールの手続きを通して, 正二十面格子(STD-grid)の格子点座標, コントロールボリュームの頂点座標,
!!     コントロールボリュームの面積を計算することができる.
!!   - 上で挙げた処理を順番に全て行うことを「STD-grid を構築する」と呼び, construct_STDGrid を呼ぶことによって行える.
!!
!! - 正二十面格子 STD-grid  の生成法について
!!   - 単位球に内包されている正二十面体の各辺を球面に投影し, 測地弧の中点をつなぎ合わせて,
!!     各球面三角形から 4 個のサブ球面三角形が生成される.
!!   - この過程を繰り返すことによって, 目標の解像度まで格子細分化を行う.
!!
!! \~english
!! \brief This module provides some subroutines which calculate coordinate data of icosahedral gridpoints and control volumes.
!!
!! - About STDGrid_Builder module
!!   - 実際の格子点座標やコントロールボリュームの頂点座標などを求める計算は, 本モジュールが提供する手続きを通して行われ,
!!     その結果は構造体 IcGrid2D_FVM の対応する成分にそれぞれ格納される.
!!   - このモジュールの手続きを通して, 正二十面格子(STD-grid)の格子点座標, コントロールボリュームの頂点座標,
!!     コントロールボリュームの面積を計算することができる.
!!   - 上で挙げた処理を順番に全て行うことを「STD-grid を構築する」と呼び, construct_STDGrid を呼ぶことによって行える.
!!
!! -  STD-Grid construction method
!!  - Each side of the icosahedron which is embedded into unit sphere are projected onto the sphere.
!!  - By connecting the midpoints of the geodesic arcs, four spherical subtriangles are generated from each triangle.
!!  - By iterating such a process until the target resolution, grid refinement proceeds.
!!
!! \~
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module STDGrid_Builder

  ! モジュール引用 ; Use statements
  !

  ! 種類型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP           ! 倍精度実数型. Double precision

  ! 数学
  ! Mathematic
  !
  use igmcore_math, only: PI     ! 円周率. PI.

  ! 線形代数
  ! Linear algebra
  !
  use igmcore_linear_algebra, only: &
    & vec_length, rotateZ

  ! 座標変換
  ! Coordinate conversion
  !
  use igmcore_coordinate_conversion, only: &
    & geo_to_orth_pos

  ! 球面三角関数
  ! Spherical trigonometry
  !
  use igmcore_spherical_trigonometry, only: &
    & spherical_penta_area, &
    & spherical_hex_area

  ! 正二十面体格子データの管理
  ! Manager of Icosahedral grid data
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, IcGrid2D_FVM_CV, &
    & paste_margin_width, calc_icosahedron_vertex, set_rcregion_bounds, &
    & get_IdMin, get_EffSize_Min, get_EffSize_Max, get_IcRadius, &
    & generate_CV5_GPindex, generate_CV6_GPindex, &
    & ICOSAHEDRON_VERTEX_NUM, RC_REGIONS_NUM

  ! 宣言文 ; Declaration statements
  !
  implicit none 
  private

  ! 公開手続き
  ! Public procedures
  !
  public :: construct_STDGrid
  public :: calc_STDGrid_GP_vertex, calc_CV_vertex, calc_CV_Area
  public :: split_rectangle_region, project_on_sphere

contains 

!
!> \~japanese
!! \brief STD-grid を構築する.
!!
!! @param[in,out] icgrid 構造型 IcGrid2D_FVM の変数.
!!
!! \~english
!! \brief Constructs the STD-grid.
!!
!! @param[in,out] icgrid The variable of derived type IcGrid2D_FVM.
!!
subroutine construct_STDGrid(icgrid)

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid

  ! 実行文 ; Executable statements
  !

  call calc_STDGrid_GP_vertex(icgrid)
  call calc_CV_vertex(icgrid)
  call calc_CV_Area(icgrid)

end subroutine construct_STDGrid

!
!> \~japanese
!! \brief STD-Grid の格子点座標を計算する.
!!
!! 求めた格子点座標は, 引数で指定された構造体 IcGrid2D_FVM の成分 rcs_AGrid に保存される.
!!
!! \param[in,out] icgrid 構造型 IcGrid2D_FVM の変数.
!!
!! \~english
!! \brief Calculates the coordinate data of STD-Grid grid points.
!!
!! 求めた格子点座標は, 引数で指定された構造体 IcGrid2D_FVM の成分 rcs_AGrid に保存される.
!!
!! \param[in,out] icgrid The variable of derived type IcGrid2D_FVM.
!!
subroutine calc_STDGrid_GP_vertex( &
  & icgrid &  ! (inout)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid

  ! 作業変数
  ! Work variables
  !
  real(DP) :: orth_icvertex(ICOSAHEDRON_VERTEX_NUM, 3)
  integer :: ringID
  integer :: rcID

  integer :: EMin, EMax
  integer :: upper_IDs(6) = (/ 2, 3, 4, 5, 6, 2 /)
  integer :: lower_IDs(6) = (/ 7, 8, 9, 10, 11, 7 /)

  ! 実行文 ; Executable statements
  !

  !write(*,*) 'calc STDGP_vertex'
  ! (単位)正二十面体の 12 個の頂点の座標を計算する.
  call calc_icosahedron_vertex(orth_icvertex)

  ! 球の半径分だけ拡大する
  orth_icvertex(:,:) = get_IcRadius(icgrid) * orth_icvertex(:,:)

  ! 各矩形領域の端にあたる格子点が, 正二十面体の頂点に重なるように, 座標を設定する.
  do ringID=1,5
    ! 上側の矩形領域に対して
    call set_rcregion_bounds( &
      & self=icgrid, &                                    ! (inout)
      & rcID=ringID, &                                      ! (in)
      & top=orth_icvertex(1,:), &                           ! (in)
      & left=orth_icvertex(upper_IDs(ringID),:), &          ! (in)
      & bottom=orth_icvertex(lower_IDs(ringID), :), &       ! (in)
      & right=orth_icvertex(upper_IDs(ringID+1),:) )        ! (in)

    ! 下側の矩形領域に対して
    call set_rcregion_bounds( &
      & self=icgrid, &                                    ! (inout)
      & rcID=ringID+5, &                                    ! (in)
      & top=orth_icvertex(upper_IDs(ringID+1),:), &         ! (in)
      & left=orth_icvertex(lower_IDs(ringID),:), &          ! (in)
      & bottom=orth_icvertex(12,:), &                       ! (in)
      & right=orth_icvertex(lower_IDs(ringID+1),:) )        ! (in)

  end do

  ! 配列インデックスの境界に関する情報を取得
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)

  ! 各矩形領域の端の格子点座標が決まったので, それをもとに再帰的に各矩形領域を分割していく.
  do rcID=1, RC_REGIONS_NUM
    call split_rectangle_region( &
      & rc_AGrid=icgrid%rcs_AGrid(:, EMin:EMax,EMin:EMax,rcID), & ! (inout)
      & ic_radius=get_IcRadius(icgrid), &                         ! (in)
      & row1=EMin, row2=EMax, &                                        ! (in)
      & col1=EMin, col2=EMax )                                         ! (in)
  end do

  ! のりしろを貼り合わせる
  call paste_margin_width(icgrid)

end subroutine calc_STDGrid_GP_vertex

!
!> \~japanese
!! \brief 正二十面格子の各格子点に付随するコントロールボリュームの頂点座標を計算する.
!!
!! 求めたコントロールボリュームの頂点座標は, 引数で指定された構造体 IcGrid2D_FVM の成分 rcs_CV に保存される.
!!
!! \param[in,out] icgrid 構造体 IcGrid2D_FVM の変数.
!!
!! \~english
!! \brief Calculates the vertex coordinate of control volumes surrounding eqch grid point of icosahedral grid.
!!
!! 求めたコントロールボリュームの頂点座標は, 引数で指定された構造体 IcGrid2D_FVM の成分 rcs_CV に保存される.
!!
!! \param[in,out] icgrid The variable of derived type IcGrid2D_FVM.
!!
subroutine calc_CV_vertex(icgrid)

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid

  ! 作業変数
  ! Work variable
  !
  type(IcGrid2D_FVM_CV), pointer :: cvInfo
  integer :: CV6_GPindx(6, 3, 2)
  integer :: CV5_GPindx(5, 3, 2)
  integer :: rcID, i, j, cvID
  integer :: EMin, EMax, idMin
  real(DP) :: ic_radius

  ! 実行文 ; Executable statemnets
  !

!  write(*,*) 'calc STDCV_vertex'

  ! 配列インデックスの境界に関する情報を取得
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)
  idMin = get_IdMin(icgrid)
  ! 半径の取得
  ic_radius = get_IcRadius(icgrid)

  ! 各格子点について, そのコントロールボリュームの頂点座標を計算する.
  do rcID=1, RC_REGIONS_NUM
    do j=EMin, EMax
      do i=EMin, EMax
         
        cvInfo => icgrid%rcs_CVInfo(i, j, rcID)
        
        if ( ( rcID > 5 .and. i==EMin .and. j==EMin  ) .or. ( i==EMin .and. j== EMax ) &
          & .or. ( i==EMax .and. j== EMin ) .or. ( rcID <= 5 .and. i==EMax .and. j== EMax ) ) then
          ! 矩形領域の左右端の特異格子点のコントロールボリュームの頂点(5個)座標を求める.

          CV5_GPindx(:,:,:) = generate_CV5_GPindex(icgrid, GP_i=i, GP_j=j, rcID=rcID )
          do cvID=1,5
            cvInfo%cv_pts(1:3, cvID) = &
              & calc_each_CV_vertex( icgrid%rcs_AGrid(:,:,:, rcID), CV5_GPindx(cvID,:,:), ic_radius, idMin )
          end do

        else if( rcID <= 5 .and. i==EMin .and. j==EMin  ) then
          ! 北極にある格子点が持つコントロールボリュームの頂点(5個)座標を求める.
          do cvID=1,5
             cvInfo%cv_pts(1:3, cvID) = &
               & (  icgrid%rcs_AGrid(:,EMin,EMin, cvID)    &
               &  + icgrid%rcs_AGrid(:,EMin+1,EMin, cvID)  &
               &  + icgrid%rcs_AGrid(:,EMin,EMin+1, cvID)  &
               &  ) / 3.0d0

             call project_on_sphere( ic_radius, cvInfo%cv_pts(1:3, cvID) )
           end do

        else if( rcID > 5 .and. i==EMax .and. j== EMax ) then
          ! 南極にある格子点が持つコントロールボリュームの頂点(5個)座標を求める.
          do cvID=1,5
             cvInfo%cv_pts(1:3, cvID) = &
               & (   icgrid%rcs_AGrid(:, EMax,EMax, 11-cvID) &
               &   + icgrid%rcs_AGrid(:, EMax-1,EMax, 11-cvID) &
               &   + icgrid%rcs_AGrid(:, EMax,EMax-1, 11-cvID) &
               &  ) / 3.0d0

             call project_on_sphere(ic_radius, cvInfo%cv_pts(1:3, cvID))
           end do

        else
          ! 特異点以外の格子点のコントロールボリュームの頂点(6個)座標を求める.
          CV6_GPindx(:,:,:) =  generate_CV6_GPindex(icgrid, GP_i=i, GP_j=j, rcID=rcID )

          do cvID=1,6
            cvInfo%cv_pts(1:3, cvID) = &
              & calc_each_CV_vertex(icgrid%rcs_AGrid(:,:,:,rcID), CV6_GPindx(cvID,:,:), ic_radius, idMin)
          end do

        end if
      end do
    end do
  end do

end subroutine calc_CV_vertex

!
!> \~japanese
!! \brief 正二十面格子の各格子点に付随するコントロールボリュームの面積を計算する.
!!
!! 求めたコントロールボリュームの面積は, 引数で指定された構造体 IcGrid2D_FVM の成分 rcs_CV に保存される.
!!
!! \param[in,out] icgrid 構造体 IcGrid2D_FVM の変数.
!!
!! \~english
!! \brief Calculates the area of the control volumes surrounding each grid point of icosahedral grid.
!!
!! 求めたコントロールボリュームの面積は, 引数で指定された構造体 IcGrid2D_FVM の成分 rcs_CVArea に保存される.
!!
!! \param[in,out] icgrid The variable of derived type IcGrid2D_FVM.
!!
subroutine calc_CV_Area( &
  & icgrid &
  & )

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid


  ! 作業変数
  ! Work variables
  !
  type(IcGrid2D_FVM_CV), pointer :: cv_info
  integer :: rcID, i, j
  integer :: EMin, EMax, idMin, k
  real(DP) :: ic_radius
  real(DP) :: CV5(3,5), CV6(3,6)

  ! 実行文 ; Executable statements
  !

!  write(*,*) 'calc CV_Area'

  ! 配列インデックスの境界に関する情報を取得
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)
  ! 半径の取得
  ic_radius = get_IcRadius(icgrid)

  ! それぞれの格子点について, それに付随するコントールボリュームの面積を求める.
  do rcID=1,RC_REGIONS_NUM
    do j=EMin, EMax
      do i=EMin, EMax

        cv_info => icgrid%rcs_CVInfo(i, j, rcID)

        if ( ( i==EMin .and. j==EMin ) .or. ( i==EMax .and. j== EMin ) .or. &
          &  ( i==EMin .and. j== EMax ) .or. ( i==EMax .and. j== EMax ) ) then
          ! 特異格子点のコントールボリューム(頂点数 5 個)に対して

          CV5(:,:) = cv_info%cv_pts(:, 1:5)

          cv_info%cv_area = &
            & spherical_penta_area( ic_radius, CV5(:,1), CV5(:,2), CV5(:,3),CV5(:,4),CV5(:,5) )

        else
          ! 通常の格子点のコントールボリューム(頂点数 6 個)に対して
         
         CV6(:,:) = cv_info%cv_pts(:, 1:6)

         cv_info%cv_area = &
            & spherical_hex_area( ic_radius, CV6(:,1), CV6(:,2), CV6(:,3),CV6(:,4),CV6(:,5), CV6(:,6) )

        end if
        
       !      if(isNAN(icgrid%rcs_CVInfo(i,j,rcID)%cv_area)) write(*,*) rcID,i,j, CV6 
      end do
    end do
  end do

end subroutine calc_CV_Area

!
!> \~japanese
!! \brief ある矩形領域を 4 分割し, 新しく生成された 5 個の格子点の座標を計算する.
!!
!! \param[in,out] rc_AGrid 対象となる矩形領域内の格子点の位置ベクトルを保持する配列.
!! \param[in] ic_radius 正二十面体格子を内包する球の半径.
!! \param[in] row1
!! \param[in] row2
!! \param[in] col1
!! \param[in] col2
!!
!! \~english
!! \brief ある矩形領域を 4 分割し, 新しく生成された 5 個の格子点の座標を計算する.
!!
!! \param[in,out] rc_AGrid The array containing the position vectors of grid-points in the target rectangle region.
!! \param[in] ic_radius The radius of a sphere which an icosahedron is embedded into.
!! \param[in] row1
!! \param[in] row2
!! \param[in] col1
!! \param[in] col2
!!
recursive subroutine split_rectangle_region( &
  & rc_AGrid,                         &  ! (inout)
  & ic_radius, row1, row2, col1, col2 &  ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(inout) :: rc_AGrid(:,:,:)
  real(DP), intent(in) :: ic_radius
  integer, intent(in) :: row1
  integer, intent(in) :: row2
  integer, intent(in) :: col1
  integer, intent(in) :: col2

  ! 作業変数
  ! Work variable
  !
  integer mid_row, mid_col

  ! 実行文 ; Executable statements
  !

  if( mod(row1+row2, 2) == 1 ) then
    ! 分割予定の格子点が存在しないので, 再帰処理を終了する.
    return
  end if

  ! 対象の矩形領域の各辺の中点と矩形領域内の中心点の座標を計算する.
  mid_row = (row1 + row2) / 2
  mid_col = (col1 + col2) / 2

  rc_AGrid(:, row1, mid_col) = 0.5d0 * ( rc_AGrid(:, row1, col1) + rc_AGrid(:, row1, col2) )
  rc_AGrid(:, row2, mid_col) = 0.5d0 * ( rc_AGrid(:, row2, col1) + rc_AGrid(:, row2, col2) )
  rc_AGrid(:, mid_row, col1) = 0.5d0 * ( rc_AGrid(:, row1, col1) + rc_AGrid(:, row2, col1) )
  rc_AGrid(:, mid_row, col2) = 0.5d0 * ( rc_AGrid(:, row1, col2) + rc_AGrid(:, row2, col2) )
  rc_AGrid(:, mid_row, mid_col) = 0.5d0 * ( rc_AGrid(:, row2, col1) + rc_AGrid(:, row1, col2) )

  ! 平面図形として矩形平面上に存在する 5 点を, 原点を中心として球面上に射影する.
  call project_on_sphere(ic_radius, rc_AGrid(:, row1, mid_col))
  call project_on_sphere(ic_radius, rc_AGrid(:, row2, mid_col))
  call project_on_sphere(ic_radius, rc_AGrid(:, mid_row, col1))
  call project_on_sphere(ic_radius, rc_AGrid(:, mid_row, col2))
  call project_on_sphere(ic_radius, rc_AGrid(:, mid_row, mid_col))

  ! 頂点の座標が決まった４つの矩形を, それぞれさらに分割する.
  call split_rectangle_region(rc_AGrid, ic_radius, row1, mid_row, col1, mid_col)
  call split_rectangle_region(rc_AGrid, ic_radius, row1, mid_row, mid_col, col2)
  call split_rectangle_region(rc_AGrid, ic_radius, mid_row, row2, col1, mid_col)
  call split_rectangle_region(rc_AGrid, ic_radius, mid_row, row2, mid_col, col2)

end subroutine split_rectangle_region

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! モジュール外非公開手続き
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!> \~japanese
!! \brief 矩形領域の平面上にある一点を原点に対して球面上へ射影する.
!!
!! @param[in]     radius 正二十面体格子を内包する球の半径.
!! @param[in,out] vec 球面に射影する点の位置ベクトル. その結果は vec に代入される.
!!
!! \~english
!! \brief 矩形領域の平面上にある一点を原点に対して球面上へ射影する.
!!
!! @param[in]     radius The radius of a sphere which an icosahedron is embedded into.
!! @param[in,out] vec 球面に射影する点の位置ベクトル. その結果は vec に代入される.
!!
subroutine project_on_sphere( &
  & radius, &  ! (in)
  & vec     &  ! (inout)
  & )

  ! 宣言文 ; Declaration statements
  !
  real(8), intent(in) :: radius    ! 正二十面体格子を内包する球の半径
  real(8), intent(inout) :: vec(3) ! 射影したい点の位置ベクトル

  ! 実行文 ; Executable statements
  !

  vec = radius / vec_length(vec) * vec
end subroutine project_on_sphere

!
!> \~japanese
!! \brief コントールボリュームの複数ある頂点のうちの一個の頂点座標を計算する.
!!
!! @param[in] idMin 配列 rc_AGrid の第 1 or 2 添字の下限インデックス.
!! @param[in] rc_AGrid 対象となるコントロールボリュームを含む矩形領域内の格子点の位置ベクトルを保持する配列.
!! @param[in] GP3_index コントロールボリュームの頂点座標の計算に必要な格子点と対応する配列添字の情報を保持する配列.
!! @param[in] ic_radius 正二十面体格子を内包する球の半径.
!! @return コントロールボリュームの頂点の位置ベクトルを保持する配列.
!!
!! \~english
!! \brief Calculates the vertex coordinate of one of control volume vertecies.
!!
!! @param[in] idMin The lower bound first (or second) index of the array rc_AGrid.
!! @param[in] rc_AGrid The array containg the position vectors data of grid-point in a rectangle region where there is the target control volume.
!! @param[in] GP3_index The array containg the information about the indecies of the array which correspond to the grid-points
!!  necessary for calculating a vertex of the control volume.
!! @param[in] ic_radius The radius of a sphere which an icosahedron is embedded into.
!! @return The array containg the information about a position vector of a vertex of the control volume.
!!
function calc_each_CV_vertex( &
  & rc_AGrid, GP3_index, ic_radius, idMin &  ! (in)
  & ) result(cv_pt)

  ! 宣言文 ; Declaration statements
  !
  integer, intent(in) :: idMin
  real(DP), intent(in) :: rc_AGrid(:,idMin:,idMin:)
  integer, intent(in) :: GP3_index(:,:)
  real(DP), intent(in) :: ic_radius
  real(DP) cv_pt(3)

  ! 作業変数
  ! Work variables
  !
  integer  k

  ! 実行文 ; Executable statements
  !

  cv_pt(:) = 0.0d0

  do k=1,3
    cv_pt(:) = cv_pt(:) + rc_AGrid(:, GP3_index(k, 1), GP3_index(k, 2))
  end do

  cv_pt(:) = cv_pt(:) / 3.0d0
  call project_on_sphere(ic_radius, cv_pt(:))

end function calc_each_CV_vertex

end module STDGrid_Builder
