!> \~japanese
!! \brief 球面三角法に関連した演算を行う手続きを提供するモジュール.
!! \~english
!! \brief This module provides some subroutines for performing the operations associated with spherical trigonometry.
!! \~
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module igmcore_spherical_trigonometry
 
  ! モジュール引用 ; Use statements
  !

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

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

  ! 宣言文 ; Declaration statements
  !

  implicit none
  private

  ! 公開手続き
  ! Public statements
  !
  public :: geodesic_arc_length
  public :: spherical_tri_area, spherical_penta_area, spherical_hex_area

  interface spherical_tri_area
    module procedure spherical_tri_area1
    module procedure spherical_tri_area2
  end interface spherical_tri_area
  
contains

!
!> \~japanese
!! \brief 球面上に存在する点 A と点 B の間の測地距離を計算する.
!!
!! @param[in] radius  球の半径.
!! @param[in] p1      直交座標系における点 A の位置ベクトル.
!! @param[in] p2      直交座標系における点 B の位置ベクトル.
!! @return            球面上に存在する 2 点間の測地距離.
!!
!! \~english
!! \brief Calculates the geodesic distance between point A and point B on a sphere.
!! 
!! @param[in] radius  The radius of a sphere.
!! @param[in] p1      A position vector of point A in cartesian coordinate system.
!! @param[in] p2      A position vector of point B in cartesian coordinate system.
!! @return            The geodesic distance between two points on a sphere.
!!
function geodesic_arc_length( &
  & radius, p1, p2            &  ! (in)
  & ) result(length)

  
  ! 宣言文 ; Declaration statements
  !
  
  real(DP), intent(in) :: radius ! 2 点が存在する球面をもつ球の半径. 
  real(DP), intent(in) :: p1(3)  ! 点 A の位置ベクトル. 
  real(DP), intent(in) :: p2(3)  ! 点 B の位置ベクトル. 
  real(DP) length                ! 球面上に存在する 2 点間の測地距離. 

  ! 作業変数
  ! Work variable
  !
  real(DP) tmp(3)

  ! 実行文 ; Executable statements
  !
  
  !length = radius * dacos( 1.0d0 - vec_length(p1(1:3) - p2(1:3))**2 / ( 2.0d0 * radius**2 ) )
  tmp(1:3) = p1(1:3) - p2(1:3)
  length = radius * acos( 1.0d0 - dot_product(tmp(1:3), tmp(1:3) ) / ( 2.0d0 * radius**2 ) )
end function geodesic_arc_length


!
!> \~japanese
!! \brief 球面上の 3 点 P1, P2, P3 によって作られる球面三角形の面積を求める.
!!
!! @param[in] radius   球の半径.
!! @param[in] orth_p1  直交座標系における点 P1 の位置ベクトル.
!! @param[in] orth_p2  直交座標系における点 P2 の位置ベクトル.
!! @param[in] orth_p3  直交座標系における点 P3 の位置ベクトル.
!! @return             球面三角形の面積.
!!
!! \~english
!! \brief Calculates the area of a spherical triangle which consist of three points(point P1, P2 and P3) on a sphere.
!!
!! @param[in] radius   The radius of a sphere.
!! @param[in] orth_p1  A position vector of point P1 in cartesian coordinate system.
!! @param[in] orth_p2  A position vector of point P2 in cartesian coordinate system.
!! @param[in] orth_p3  A position vector of point P3 in cartesian coordinate system.
!! @return             The area of spherical triangle.
!! 
function spherical_tri_area1( &
  & radius, orth_p1, orth_p2, orth_p3 &  ! (in)
  & ) result(area)

  ! 宣言文 ; Declaration statements
  !

  real(DP), intent(in) :: radius         ! 三点 P1, P2, P3 が球面に存在する球の半径.
  real(DP), intent(in) :: orth_p1(3)     ! 点 P1 の位置ベクトル(ベクトルの成分は直交座標系における成分で表す).
  real(DP), intent(in) :: orth_p2(3)     ! 点 P2 の位置ベクトル(ベクトルの成分は直交座標系における成分で表す).
  real(DP), intent(in) :: orth_p3(3)     ! 点 P3 の位置ベクトル(ベクトルの成分は直交座標系における成分で表す).
  real(DP) :: area                       ! 球面三角形の面積.

  ! 作業変数
  ! Work variables
  !
  real(DP) arc1, arc2, arc3

  ! 実行文 ; Executable statements
  !

  ! 球面三角形の各辺(各弧)の長さ(単位球の中心に対する角度)を計算する
  arc1 = geodesic_arc_length_angle(radius, orth_p2, orth_p3)
  arc2 = geodesic_arc_length_angle(radius, orth_p3, orth_p1)
  arc3 = geodesic_arc_length_angle(radius, orth_p1, orth_p2)

  ! 球面三角法の公式に基づき, 球面三角形の面積を求める
  area = ( calc_angle(arc1, arc2, arc3)    &
    &      + calc_angle(arc2, arc3, arc1)  &
    &      + calc_angle(arc3, arc1, arc2)  &
    &      - PI                            &
    &     ) * radius**2

end function spherical_tri_area1

!
!> \~japanese
!! \brief 球面上の 3 点 P1, P2, P3 によって作られる球面三角形の面積を求める.
!!
!! @param[in] radius   球の半径.
!! @param[in] orth_p1  直交座標系における点 P1 の位置ベクトル.
!! @param[in] orth_p2  直交座標系における点 P2 の位置ベクトル.
!! @param[in] orth_p3  直交座標系における点 P3 の位置ベクトル.
!! @return             球面三角形の面積.
!!
!! \~english
!! \brief Calculates the area of a spherical triangle which consist of three points(point P1, P2 and P3) on a sphere.
!!
!! @param[in] radius   The radius of a sphere.
!! @param[in] orth_p1  A position vector of point P1 in cartesian coordinate system.
!! @param[in] orth_p2  A position vector of point P2 in cartesian coordinate system.
!! @param[in] orth_p3  A position vector of point P3 in cartesian coordinate system.
!! @return             The area of spherical triangle.
!!
function spherical_tri_area2( &
  & radius, orth_pts           &  ! (in)
  & ) result(area)

  ! 宣言文 ; Declaration statements
  !

  real(DP), intent(in) :: radius           ! 三点 P1, P2, P3 が球面に存在する球の半径.
  real(DP), intent(in) :: orth_pts(3,3)    !  点 P1, P2, P3 の位置ベクトル(ベクトルの成分は直交座標系における成分で表す).
  real(DP) :: area                         ! 球面三角形の面積.

  ! 実行文 ; Executable statements
  !

  area = spherical_tri_area1(radius, orth_pts(:,1), orth_pts(:,2), orth_pts(:,3))

end function spherical_tri_area2

!
!> \~japanese
!! \brief 球面上の 5 点 P1, P2, P3, P4, P4, P5によって作られる球面五角形の面積を求める.
!!
!! @param[in] radius   球の半径.
!! @param[in] orth_p1  直交座標系における点 P1 の位置ベクトル.
!! @param[in] orth_p2  直交座標系における点 P2 の位置ベクトル.
!! @param[in] orth_p3  直交座標系における点 P3 の位置ベクトル.
!! @param[in] orth_p4  直交座標系における点 P4 の位置ベクトル.
!! @param[in] orth_p5  直交座標系における点 P5 の位置ベクトル.
!! @return             球面五角形の面積.
!!
!! \~english
!! \brief Calculates the area of a spherical pentagon which consist of five points(point P1, P2, P3, P4 and P5) on a sphere.
!!
!! @param[in] radius   The radius of a sphere.
!! @param[in] orth_p1  A position vector of point P1 in cartesian coordinate system.
!! @param[in] orth_p2  A position vector of point P2 in cartesian coordinate system.
!! @param[in] orth_p3  A position vector of point P3 in cartesian coordinate system.
!! @param[in] orth_p4  A position vector of point P4 in cartesian coordinate system.
!! @param[in] orth_p5  A position vector of point P5 in cartesian coordinate system.
!! @return             The area of spherical pentagon.
!! 
function spherical_penta_area( &
  & radius,                    &                   ! (in)
  & orth_p1, orth_p2, orth_p3, orth_p4, orth_p5 &  ! (in)
  ) result(area)

  ! 宣言文 ; Declaration statements
  !

  real(DP), intent(in) :: radius                              ! 5 点が球面に存在する球の半径
  real(DP), intent(in) :: orth_p1(3), orth_p2(3), orth_p3(3), orth_p4(3), orth_p5(3)   ! 5 点それぞれの位置ベクトル
  real(DP) area                                               ! 球面五角形の面積

  area = spherical_tri_area(radius, orth_p1, orth_p2, orth_p3) + spherical_tri_area(radius, orth_p1, orth_p3, orth_p4) &
    &    + spherical_tri_area(radius, orth_p1, orth_p4, orth_p5)
end function spherical_penta_area

!
!> \~japanese
!! \brief 球面上の 6 点 P1, P2, P3, P4, P4, P5, P6 によって作られる球面六角形の面積を求める.
!!
!! @param[in] radius   球の半径.
!! @param[in] orth_p1  直交座標系における点 P1 の位置ベクトル.
!! @param[in] orth_p2  直交座標系における点 P2 の位置ベクトル.
!! @param[in] orth_p3  直交座標系における点 P3 の位置ベクトル.
!! @param[in] orth_p4  直交座標系における点 P4 の位置ベクトル.
!! @param[in] orth_p5  直交座標系における点 P5 の位置ベクトル.
!! @param[in] orth_p5  直交座標系における点 P6 の位置ベクトル.
!! @return             球面五角形の面積.
!!
!! \~english
!! \brief Calculates the area of a spherical hexagon which consist of five points(point P1, P2, P3, P4, P5 and P6) on a sphere.
!!
!! @param[in] radius   The radius of a sphere.
!! @param[in] orth_p1  A position vector of point P1 in cartesian coordinate system.
!! @param[in] orth_p2  A position vector of point P2 in cartesian coordinate system.
!! @param[in] orth_p3  A position vector of point P3 in cartesian coordinate system.
!! @param[in] orth_p4  A position vector of point P4 in cartesian coordinate system.
!! @param[in] orth_p5  A position vector of point P5 in cartesian coordinate system.
!! @param[in] orth_p5  A position vector of point P6 in cartesian coordinate system.
!! @return             The area of spherical hexagon.
!! 
function spherical_hex_area( &
  & radius, &                                               ! (in)
  & orth_p1, orth_p2, orth_p3, orth_p4, orth_p5, orth_p6 &  ! (in)
  & ) result(area)

  ! 宣言文 ; Declaration statements
  !

  real(DP), intent(in) :: radius                                   ! 6 点が球面に存在する球の半径
  real(DP), intent(in) :: orth_p1(3), orth_p2(3), orth_p3(3), &    !
    &                     orth_p4(3), orth_p5(3), orth_p6(3)       ! 6 点の位置ベクトル

  real(DP) area                                                    ! 球面六角形の面積

  ! 実行文 ; Executable statements
  !

  ! 球面六角形を 3 個の球面三角形として, 面積を求める. 
  area =   spherical_tri_area(radius, orth_p1, orth_p2, orth_p3) &
    &    + spherical_tri_area(radius, orth_p1, orth_p3, orth_p4) &
    &    + spherical_tri_area(radius, orth_p1, orth_p4, orth_p5) &
    &    + spherical_tri_area(radius, orth_p1, orth_p5, orth_p6)

end function spherical_hex_area

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

!
!> \~japanese
!! \brief 球面上における 2 点(A, B)間の距離を単位球の中心に対する角度として計算する.
!!
!! @param[in] radius 2 点が存在する球面をもつ球の半径.
!! @param[in] p1 点 A の位置ベクトル.
!! @param[in] p2 点 B の位置ベクトル.
!! @return 球面上の二点間の距離(単位球の中心に対する角度).
!!
!! \~english
!! \brief 球面上における 2 点間の距離を単位球の中心に対する角度として計算する.
!!
!! @param[in] radius 2 点が存在する球面をもつ球の半径.
!! @param[in] p1 点 A の位置ベクトル.
!! @param[in] p2 点 B の位置ベクトル.
!! @return 球面上の二点間の距離(単位球の中心に対する角度).
!!
function geodesic_arc_length_angle( &
  & radius, p1, p2 &   ! (in)
  & ) result(length)

  ! 宣言文 ; Declaration statements
  !

  real(DP), intent(in) :: radius ! 2 点が存在する球面をもつ球の半径
  real(DP), intent(in) :: p1(3)  ! 点の位置ベクトル
  real(DP), intent(in) :: p2(3)  ! 点の位置ベクトル
  real(DP) length                ! 球面上の二点間の距離(単位球の中心に対する角度)

  ! 作業変数
  ! Work variables
  !
  real(DP) tmp(3)

  !length = dacos( 1.0d0 - vec_length(p1(1:3) - p2(1:3))**2 / ( 2.0d0 * radius**2 ) )
  tmp(1:3) = p1(1:3) - p2(1:3)
  length = dacos( 1.0d0 - dot_product(tmp(1:3), tmp(1:3) ) / ( 2.0d0 * radius**2 ) )
end function geodesic_arc_length_angle

!
!> \~japanese
!! \brief 球面三角形 ABC の一角である角 BAC を 3 辺の弧の長さ(この弧の長さは単位球の中心に対する角度を使って表す)から求める.
!!
!! @param[in] coresp_arc 対象となる角と向かい合う弧 BC の長さ.
!! @param[in] arc1 対象となる角を挟む弧 AB の長さ.
!! @param[in] arc2 対象となる角を挟む弧 AC の長さ.
!! @return 角 BAC の角度.
!!
!! \~english
!! \brief 球面三角形 ABC の一角である角 BAC を 3 辺の弧の長さ(この弧の長さは単位球の中心に対する角度を使って表す)から求める.
!!
!! @param[in] coresp_arc 対象となる角と向かい合う弧 BC の長さ.
!! @param[in] arc1 対象となる角を挟む弧 AB の長さ.
!! @param[in] arc2 対象となる角を挟む弧 AC の長さ.
!! @return 角 BAC の角度.
!!
function calc_angle( &
  & coresp_arc, arc1, arc2 &  ! (in)
  & ) result(angle)


  ! 宣言文 ; Declaration statements
  !

  real(DP), intent(in) :: coresp_arc  ! 対象となる角と向かい合う弧の長さ
  real(DP), intent(in) :: arc1        ! 対象となる角を挟む弧の長さ
  real(DP), intent(in) :: arc2        ! 対象となる角を挟む弧の長さ
  real(DP) angle                      ! 求めた角度

  ! 実行文 ; Executable statements
  !
  
  angle = dacos( &
      ( dcos(coresp_arc) - dcos(arc1) * dcos(arc2) ) / ( dsin(arc1) * dsin(arc2) ) &
    )

end function calc_angle

end module igmcore_spherical_trigonometry
