src/core/math/igmcore_spherical_trigonometry.f90

Go to the documentation of this file.

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

contains

!
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

!
function spherical_tri_area( &
  & 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_area

!
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

!
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

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

!
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

!
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