動作環境

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)の範囲とする。

関連

式の整理

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なのかもしれない。