IGMBaseLib 1.0
|
00001 00012 module igmcore_spherical_trigonometry 00013 00014 ! モジュール引用 ; Use statements 00015 ! 00016 00017 ! 種類型パラメタ 00018 ! Kind type parameter 00019 ! 00020 use dc_types, only: DP ! 倍精度実数型. Double precision. 00021 00022 ! 数学 00023 ! Mathematic 00024 ! 00025 use igmcore_math, only: PI ! 円周率 00026 00027 ! 宣言文 ; Declaration statements 00028 ! 00029 00030 implicit none 00031 private 00032 00033 ! 公開手続き 00034 ! Public statements 00035 ! 00036 public :: geodesic_arc_length 00037 public :: spherical_tri_area, spherical_penta_area, spherical_hex_area 00038 00039 contains 00040 00041 ! 00058 function geodesic_arc_length( & 00059 & radius, p1, p2 & ! (in) 00060 & ) result(length) 00061 00062 00063 ! 宣言文 ; Declaration statements 00064 ! 00065 00066 real(DP), intent(in) :: radius ! 2 点が存在する球面をもつ球の半径. 00067 real(DP), intent(in) :: p1(3) ! 点 A の位置ベクトル. 00068 real(DP), intent(in) :: p2(3) ! 点 B の位置ベクトル. 00069 real(DP) length ! 球面上に存在する 2 点間の測地距離. 00070 00071 ! 作業変数 00072 ! Work variable 00073 ! 00074 real(DP) tmp(3) 00075 00076 ! 実行文 ; Executable statements 00077 ! 00078 00079 !length = radius * dacos( 1.0d0 - vec_length(p1(1:3) - p2(1:3))**2 / ( 2.0d0 * radius**2 ) ) 00080 tmp(1:3) = p1(1:3) - p2(1:3) 00081 length = radius * acos( 1.0d0 - dot_product(tmp(1:3), tmp(1:3) ) / ( 2.0d0 * radius**2 ) ) 00082 end function geodesic_arc_length 00083 00084 ! 00103 function spherical_tri_area( & 00104 & radius, orth_p1, orth_p2, orth_p3 & ! (in) 00105 & ) result(area) 00106 00107 ! 宣言文 ; Declaration statements 00108 ! 00109 00110 real(DP), intent(in) :: radius ! 三点 P1, P2, P3 が球面に存在する球の半径. 00111 real(DP), intent(in) :: orth_p1(3) ! 点 P1 の位置ベクトル(ベクトルの成分は直交座標系における成分で表す). 00112 real(DP), intent(in) :: orth_p2(3) ! 点 P2 の位置ベクトル(ベクトルの成分は直交座標系における成分で表す). 00113 real(DP), intent(in) :: orth_p3(3) ! 点 P3 の位置ベクトル(ベクトルの成分は直交座標系における成分で表す). 00114 real(DP) area ! 球面三角形の面積. 00115 00116 ! 作業変数 00117 ! Work variables 00118 ! 00119 real(DP) arc1, arc2, arc3 00120 00121 ! 実行文 ; Executable statements 00122 ! 00123 00124 ! 球面三角形の各辺(各弧)の長さ(単位球の中心に対する角度)を計算する 00125 arc1 = geodesic_arc_length_angle(radius, orth_p2, orth_p3) 00126 arc2 = geodesic_arc_length_angle(radius, orth_p3, orth_p1) 00127 arc3 = geodesic_arc_length_angle(radius, orth_p1, orth_p2) 00128 00129 ! 球面三角法の公式に基づき, 球面三角形の面積を求める 00130 area = ( calc_angle(arc1, arc2, arc3) & 00131 & + calc_angle(arc2, arc3, arc1) & 00132 & + calc_angle(arc3, arc1, arc2) & 00133 & - PI & 00134 & ) * radius**2 00135 00136 end function spherical_tri_area 00137 00138 ! 00161 function spherical_penta_area( & 00162 & radius, & ! (in) 00163 & orth_p1, orth_p2, orth_p3, orth_p4, orth_p5 & ! (in) 00164 ) result(area) 00165 00166 ! 宣言文 ; Declaration statements 00167 ! 00168 00169 real(DP), intent(in) :: radius ! 5 点が球面に存在する球の半径 00170 real(DP), intent(in) :: orth_p1(3), orth_p2(3), orth_p3(3), orth_p4(3), orth_p5(3) ! 5 点それぞれの位置ベクトル 00171 real(DP) area ! 球面五角形の面積 00172 00173 area = spherical_tri_area(radius, orth_p1, orth_p2, orth_p3) + spherical_tri_area(radius, orth_p1, orth_p3, orth_p4) & 00174 & + spherical_tri_area(radius, orth_p1, orth_p4, orth_p5) 00175 end function spherical_penta_area 00176 00177 ! 00202 function spherical_hex_area( & 00203 & radius, & ! (in) 00204 & orth_p1, orth_p2, orth_p3, orth_p4, orth_p5, orth_p6 & ! (in) 00205 & ) result(area) 00206 00207 ! 宣言文 ; Declaration statements 00208 ! 00209 00210 real(DP), intent(in) :: radius ! 6 点が球面に存在する球の半径 00211 real(DP), intent(in) :: orth_p1(3), orth_p2(3), orth_p3(3), ! 00212 orth_p4(3), orth_p5(3), orth_p6(3) ! 6 点の位置ベクトル 00213 00214 real(DP) area ! 球面六角形の面積 00215 00216 ! 実行文 ; Executable statements 00217 ! 00218 00219 ! 球面六角形を 3 個の球面三角形として, 面積を求める. 00220 area = spherical_tri_area(radius, orth_p1, orth_p2, orth_p3) & 00221 & + spherical_tri_area(radius, orth_p1, orth_p3, orth_p4) & 00222 & + spherical_tri_area(radius, orth_p1, orth_p4, orth_p5) & 00223 & + spherical_tri_area(radius, orth_p1, orth_p5, orth_p6) 00224 00225 end function spherical_hex_area 00226 00227 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00228 ! 00229 ! モジュール外非公開手続き 00230 ! 00231 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00232 00233 ! 00250 function geodesic_arc_length_angle( & 00251 & radius, p1, p2 & ! (in) 00252 & ) result(length) 00253 00254 ! 宣言文 ; Declaration statements 00255 ! 00256 00257 real(DP), intent(in) :: radius ! 2 点が存在する球面をもつ球の半径 00258 real(DP), intent(in) :: p1(3) ! 点の位置ベクトル 00259 real(DP), intent(in) :: p2(3) ! 点の位置ベクトル 00260 real(DP) length ! 球面上の二点間の距離(単位球の中心に対する角度) 00261 00262 ! 作業変数 00263 ! Work variables 00264 ! 00265 real(DP) tmp(3) 00266 00267 !length = dacos( 1.0d0 - vec_length(p1(1:3) - p2(1:3))**2 / ( 2.0d0 * radius**2 ) ) 00268 tmp(1:3) = p1(1:3) - p2(1:3) 00269 length = dacos( 1.0d0 - dot_product(tmp(1:3), tmp(1:3) ) / ( 2.0d0 * radius**2 ) ) 00270 end function geodesic_arc_length_angle 00271 00272 ! 00289 function calc_angle( & 00290 & coresp_arc, arc1, arc2 & ! (in) 00291 & ) result(angle) 00292 00293 00294 ! 宣言文 ; Declaration statements 00295 ! 00296 00297 real(DP), intent(in) :: coresp_arc ! 対象となる角と向かい合う弧の長さ 00298 real(DP), intent(in) :: arc1 ! 対象となる角を挟む弧の長さ 00299 real(DP), intent(in) :: arc2 ! 対象となる角を挟む弧の長さ 00300 real(DP) angle ! 求めた角度 00301 00302 ! 実行文 ; Executable statements 00303 ! 00304 00305 angle = dacos( & 00306 ( dcos(coresp_arc) - dcos(arc1) * dcos(arc2) ) / ( dsin(arc1) * dsin(arc2) ) & 00307 ) 00308 00309 end function calc_angle 00310 00311 end module igmcore_spherical_trigonometry