動作環境
GeForce GTX 1070 (8GB)
ASRock Z170M Pro4S [Intel Z170chipset]
Ubuntu 16.04 LTS desktop amd64
TensorFlow v1.2.1
cuDNN v5.1 for Linux
CUDA v8.0
Python 3.5.2
IPython 6.0.0 -- An enhanced Interactive Python.
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609
GNU bash, version 4.3.48(1)-release (x86_64-pc-linux-gnu)
scipy v0.19.1
geopandas v0.3.0
MATLAB R2017b (Home Edition)
x,y,z要素が与えられたベクトルから天頂角と方位角を計算する。
天頂角は[0, pi)の範囲で、方位角は[0, 2pi)の範囲とする。
関連
- ADDA > Euler angle rotation > beta, gammaに対して計算する
- geometry > Euler angle rotation > beta, gamma に関する式 > 天頂角と方位角の式
式の整理
x = cos\gamma sin\beta ... (1) \\
y = sin\gamma sin\beta ... (2) \\
z = cos\beta ...(3)
方向の定義
方向について整理しておく。
- $\beta$ = 0のとき
- x = 0, y = 0, z = 1
-
z axisに対応する
- $\gamma$ = 0のとき
- x = $sin\beta$
- y = 0
- (z = $cos\beta$)
-
x,y平面上ではx axisに対応する
betaの式
式(3)より
\beta = arccos(z) ... (4)
numpy.arccos()を使った場合、betaの範囲は[0, pi)となる。
gammaの式
(追記 2017/12/10)
後述のNumpy実装 v0.2にて下記の式ではなく、arctan2()を使うように変更した
式(2)/式(1)より
\frac{y}{x} = \frac{\sin\gamma}{\cos\gamma} = tan\gamma ... (5.1)
\gamma = arctan(\frac{y}{x}) ... (5.2)
numpy.arctan()を使った場合、gammaの範囲は[-pi/2, pi/2)となる。
象限の扱いが半分だけとなり、希望している[0, 2pi)にはならない。
条件(x<0)を加える。
\gamma = \pi - arctan(\frac{y}{x}) \\
for (x<0) ... (5.3)
x=0の場合はy軸上になるため、下記とする。
\gamma = arctan(y) \\
for (x=0) ... (5.4)
Numpy実装 v0.1
polarAzimuthCalc_171209.py
import numpy as np
import sys
'''
v0.1 Dec.09, 2017
- add Test_group_run()
- add calc_gamma()
- add calc_beta()
'''
def calc_beta(pvec):
'''
polar angle [0, pi]
'''
return np.arccos(pvec[2]) # arccos:[0, pi]
def calc_gamma(pvec):
'''
azimuth angle [0, 2pi]
'''
if abs(pvec[0]) < sys.float_info.epsilon:
gamma = np.arctan(pvec[1]) # arctan: [-pi/2, pi/2]
else:
gamma = np.arctan(pvec[1]/pvec[0]) # arctan: [-pi/2, pi/2]
if pvec[0] < 0.0:
gamma = (np.pi - gamma)
if gamma < 0.0:
gamma += 2 * np.pi
return gamma
def Test_group_run():
pvecs = [
[0, 0, 1], # z direction
[0, 0, -1], # z direction
[0, 1, 0], # y direction
[0, -1, 0], # y direction
[1, 0, 0], # x direction
[-1, 0, 0], # x direction
]
for elem in pvecs:
beta_rad = calc_beta(elem)
gamme_rad = calc_gamma(elem)
fmt = "{0} \tbeta:{1:.5f} gamma:{2:.5f}"
msg = fmt.format(elem, beta_rad, gamme_rad)
print(msg)
if __name__ == '__main__':
Test_group_run()
$ python3 polarAzimuthCalc_171209.py
[0, 0, 1] beta:0.00000 gamma:0.00000
[0, 0, -1] beta:3.14159 gamma:0.00000
[0, 1, 0] beta:1.57080 gamma:0.78540
[0, -1, 0] beta:1.57080 gamma:5.49779
[1, 0, 0] beta:1.57080 gamma:0.00000
[-1, 0, 0] beta:1.57080 gamma:3.14159
Numpy実装v0.2
- gamma, betaの値からx,y,zを計算し、それらからgamma, betaを計算する
- 計算の確認のため
- gamma計算はarctan2()を使うように変更した
関連: https://vvvv.org/blog/polar-spherical-and-geographic-coordinates
polarAzimuthCalc_171209.py
import numpy as np
import sys
'''
v0.2 Dec. 10, 2017
- fix bug: calc_gamma() was too complicated
_+ use arctan2()
- add Test_group_run_random()
v0.1 Dec. 09, 2017
- add Test_group_run_axes()
- add calc_gamma()
- add calc_beta()
'''
def calc_beta(pvec):
'''
polar angle [0, pi]
'''
return np.arccos(pvec[2]) # arccos:[0, pi]
def calc_gamma(pvec):
gamma = np.arctan2(pvec[1], pvec[0])
if gamma < 0.0:
gamma += 2 * np.pi
return gamma
def Test_group_run_axes():
pvecs = [
[0, 0, 1], # z direction
[0, 0, -1], # z direction
[0, 1, 0], # y direction
[0, -1, 0], # y direction
[1, 0, 0], # x direction
[-1, 0, 0], # x direction
]
for elem in pvecs:
beta_rad = calc_beta(elem)
gamma_rad = calc_gamma(elem)
fmt = "{0} \tbeta:{1:.5f} gamma:{2:.5f}"
msg = fmt.format(elem, beta_rad, gamma_rad)
print(msg)
def calc_xyz(beta_rad, gamma_rad):
rad = 1.0 # radius
resx = rad * np.cos(gamma_rad) * np.sin(beta_rad)
resy = rad * np.sin(gamma_rad) * np.sin(beta_rad)
resz = rad * np.cos(beta_rad)
return resx, resy, resz
def Test_group_run_random():
bt_stp = np.pi / 4.0 # beta step
gm_stp = np.pi / 4.0 # gamma step
for bt_in in np.arange(0.0, np.pi, bt_stp):
for gm_in in np.arange(0.0, 2 * np.pi, gm_stp):
# print(bt_in, gm_in)
wrk = calc_xyz(bt_in, gm_in)
bt_out = calc_beta(wrk)
gm_out = calc_gamma(wrk)
if abs(gm_in - gm_out) > sys.float_info.epsilon:
print("%+.5e" % (gm_in - gm_out), end=' ')
print("NG:", end=" ")
else:
continue
msg = ','.join("%+.5f" % elem for elem in wrk)
print(msg, end=' -- ')
fmt = "b{0:.5f} b{1:.5f} g{2:.5f} g{3:.5f}"
msg = fmt.format(bt_in, bt_out, gm_in, gm_out)
print(msg)
if __name__ == '__main__':
Test_group_run_random()
#Test_group_run_axes()
$ python3 polarAzimuthCalc_171209.py | less -SX
+7.85398e-01 NG: +0.00000,+0.00000,+1.00000 -- b0.00000 b0.00000 g0.78540 g0.00000
+1.57080e+00 NG: +0.00000,+0.00000,+1.00000 -- b0.00000 b0.00000 g1.57080 g0.00000
-7.85398e-01 NG: -0.00000,+0.00000,+1.00000 -- b0.00000 b0.00000 g2.35619 g3.14159
+7.85398e-01 NG: -0.00000,-0.00000,+1.00000 -- b0.00000 b0.00000 g3.92699 g3.14159
+1.57080e+00 NG: -0.00000,-0.00000,+1.00000 -- b0.00000 b0.00000 g4.71239 g3.14159
+5.49779e+00 NG: +0.00000,-0.00000,+1.00000 -- b0.00000 b0.00000 g5.49779 g-0.00000
+4.44089e-16 NG: -0.50000,-0.50000,+0.70711 -- b0.78540 b0.78540 g3.92699 g3.92699
+4.44089e-16 NG: -0.70711,-0.70711,+0.00000 -- b1.57080 b1.57080 g3.92699 g3.92699
上記のエラーが出た。
- (0, 0, 1)の時のgammaでエラー
- これは縮退するので一つのgammaしか求まらない
- 下の２つで誤差4.4e-16
- sys.float_info_epsilonは2.220e-16
- その2倍程度の誤差が生じている
answered Mar 2 '12 at 5:51
Ergwun
But don't forget the pitfalls of using it as an absolute error margin for floating point comparisons. E.g. for large numbers, rounding error could exceed epsilon.
pitfallsなのかもしれない。