geometry > x,y,zベクトルから天頂角と方位角を計算する > 式の整理 | Numpy実装 polarAzimuthCalc_171209.py > v0.1, v0.2 | sys.float_info_epsilonのpitfalls

```txt:動作環境
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に対して計算する](https://qiita.com/7of9/items/717c08d9c161e2d66f23)
- [geometry > Euler angle rotation > beta, gamma に関する式 > 天頂角と方位角の式](https://qiita.com/7of9/items/8dc230ecec4e51aed18e)

## 式の整理

```math
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に対応する

![DSC_0270.JPG](https://qiita-image-store.s3.amazonaws.com/0/32870/313901c0-8f63-7a45-4947-5c47e229bba4.jpeg)


### betaの式

式(3)より

```math
\beta = arccos(z) ... (4)
```

[numpy.arccos()](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.arccos.html)を使った場合、betaの範囲は[0, pi)となる。

### gammaの式
**(追記 2017/12/10)**
**後述のNumpy実装 v0.2にて下記の式ではなく、arctan2()を使うように変更した**


式(2)/式(1)より

```math
\frac{y}{x} = \frac{\sin\gamma}{\cos\gamma} = tan\gamma ... (5.1)
```

```math
\gamma = arctan(\frac{y}{x}) ... (5.2)
```

[numpy.arctan()](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.arctan.html)を使った場合、gammaの範囲は[-pi/2, pi/2)となる。
象限の扱いが半分だけとなり、希望している[0, 2pi)にはならない。

条件(x<0)を加える。

```math
\gamma = \pi - arctan(\frac{y}{x}) \\
 for (x<0) ... (5.3)
```

x=0の場合はy軸上になるため、下記とする。

```math
\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()

```

```bash
$ 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://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.arctan2.html)を使うように変更した

関連: 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()

```

```bash
$ 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倍程度の誤差が生じている

https://stackoverflow.com/questions/9528421/value-for-epsilon-in-python
> 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なのかもしれない。