geometry > calculate polar and azimuth angles from (x,y,z) vector > Equations and Numpy Implementation > polarAzimuthCalc_171209.py > v0.1..v0.3

16911 ワード
ADDA borgWarp geometry pySpherepts geometry テキストリンク
Environment
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)
ADDA v.1.3b6

This is an English version of this article.

The equation of gamma is modified to use arctan2().
In addition, the code is modified (refactored) a little.

About

We introduce the calculation of the zenith and azimuth angles from a (x,y,z) vector.
The zenith angle is in [0, pi) and the azimuth is in [0, 2pi), where [ and ( denote open and closed, respectively.

Equations

Definition of the zenith and azimuth angles

We define the zenith and azimuth angles as $\beta$ and $\gamma$.

The (x,y,z) and the ($\beta$ and $\gamma$) is related as follows:

x = cos\gamma sin\beta ... (1) \\
y = sin\gamma sin\beta ... (2) \\
z = cos\beta ...(3)

Definition of the directions

We define $\beta$ and $\gamma$ as follows in terms of the directions.

  • for $\beta$ = 0
    • x = 0, y = 0, z = 1
    • > corresponds to z axis
  • for $\gamma$ = 0
    • x = $sin\beta$
    • y = 0
    • (z = $cos\beta$)
    • > corresponds to x axis in (x,y) plane

Equation of beta

From Eq.(3)

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

By using numpy.arccos(), the range of beta is [0, pi).

Equation of gamma

From Eq.(2) divided by Eq.(1)

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

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

numpy.arctan2() returns values in [-pi, pi).
The value is converted to [0, 2pi) for the convenience of the later data processing.

Numpy Implemantation v0.1..v0.3

polarAzimuthCalc_171209.py
import numpy as np
import sys

'''
v0.3 Dec. 17, 2017
    - refactor > rename functions by adding "_rad"
    - refactor > add function comment
v0.2 Dec. 10, 2017
    - fix bug: calc_gamma_rad() was too complicated
        _+ use arctan2()
    - add Test_group_run_random()
v0.1 Dec. 09, 2017
    - add Test_group_run_axes()
    - add calc_gamma_rad()
    - add calc_beta_rad()
'''


def calc_beta_rad(pvec):
    '''
    polar angle [0, pi]
    '''
    return np.arccos(pvec[2])  # arccos:[0, pi]


def calc_gamma_rad(pvec):
    '''
    azimuth angle [0, 2pi]
    '''
    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_rad(elem)
        gamma_rad = calc_gamma_rad(elem)

        fmt = "{0} \tbeta:{1:.5f} gamma:{2:.5f}"
        msg = fmt.format(elem, beta_rad, gamma_rad)
        print(msg)


def Test_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 = Test_calc_xyz(bt_in, gm_in)
            bt_out = calc_beta_rad(wrk)
            gm_out = calc_gamma_rad(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()

Run example

The above code executes the Test_group_run_random().
In the Test_group_run_random(), $\beta$ and $\gamma$ are taken evenly in the each dimension.

$ python3 polarAzimuthCalc_171209.py 
+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

There are several errors.

  • (0, 0, 1) returns gamma errors
    • Since the value of gamma is degenerated
  • error of 4.4e-16 in the last two lines
    • sys.float_info_epsilon is 2.220e-16
    • the errors (4.4e-16) is twice of the sys.float_info_epsilon

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.

Maybe, pitfalls...