pythonで多次元の極座標変換

#はじめに
3次元以上のベクトルに対して直交座標から極座標へ変換するメソッドが見つからなかったので自分で作成しました。せっかくなので公開します。pythonのライブラリで直交座標から極座標へ変換するメソッドが既にある場合には紹介していただけると嬉しいです。
なお、このメソッドを用いて生じた問題については一切責任は負いません。

#直交座標と極座標の関係
まずは直交座標と極座標の関係を数式で表現します。多次元の極座標変換を参考にしました。
$n \in \mathbb{N}$として、$(x_1, \cdots, x_n) \in \mathbb{R}^n, (r, \theta_1, \cdots, \theta_{n-1}) \in \mathbb{R}_{\geq 0}\times [-\pi, \pi)^{n-1}$とします。このとき、直交座標と極座標の変換公式は、

\begin{align}
x_1 &= r{\rm cos}\theta_1 {\rm cos}\theta_2 \cdots {\rm cos}\theta_{n-2}{\rm cos}\theta_{n-1} \\
x_2 &= r{\rm cos}\theta_1 {\rm cos}\theta_2 \cdots {\rm cos}\theta_{n-2}{\rm sin}\theta_{n-1} \\
x_3 &= r{\rm cos}\theta_1 {\rm cos}\theta_2 \cdots {\rm sin}\theta_{n-2} \\
\vdots \\
x_n &=  r{\rm sin}\theta_1
\end{align}

です。これを$(r, \theta_1, \cdots, \theta_{n-1})$について解くと、

\begin{align}
r &= \sqrt{\sum_{i=1}^n x_i^2} \\
\theta_1 &= {\rm arctan}\left(\frac{x_{n}{\rm sin}\theta_2}{x_{n-1}}\right) \\
\theta_2 &= {\rm arctan}\left(\frac{x_{n-1}{\rm sin}\theta_3}{x_{n-2}}\right) \\
\vdots \\
\theta_{n-1} &= {\rm arctan}\left(\frac{x_2}{x_1} \right)
\end{align}

となります。この変換公式を使って、直交座標から極座標に変換するメソッドをpythonで作成します。

#実装
##方針
引数が直交座表現のベクトル、出力が極座標表現のベクトルとなるメソッドを作成します。なお、引数、返り値ともにNumPy.arrayとします。
座標系の変換をするにあたり、以下の場合分けを行います。

1. ノルムが0のとき
   1. 入力されたベクトルの次元数が1以下のとき
   2. それ以外
2. 入力されたベクトルの次元数が1のとき
3. 入力されたベクトルの次元数が2のとき
4. それ以外

上記のそれぞれの場合について、返り値を以下のように定めます。

1-1. 文字列「2次元以上のベクトルを入力してください」を表示し、Noneを返す。
1-2. 入力されたベクトルの同じ大きさのゼロベクトルを返す。
2. 文字列「2次元以上のベクトルを入力してください」を表示し、Noneを返す。
3. 極座標ベクトル $\left(\sqrt{x_1^2+x_2^2},\ {\rm arctan}\left(\frac{x_2}{x_1}\right) \right)$ を返す。
4. 極座標ベクトル $\left(\sqrt{x_1^2+x_2^2},\ {\rm arctan}\left(\frac{x_{n}{\rm sin}\theta_2}{x_{n-1}}\right),\cdots, \ {\rm arctan}\left(\frac{x_2}{x_1}\right) \right)$ を返す。

##実装例
上記の場合分けに基づいて直交座標から極座標に変換するメソッドは下記です。

get_pc.py

#import
import numpy as np
from math import atan2, sin

#直交座標⇒極座標変換メソッド
def get_pc(cartesian_coordinate_vec):
    
    cartesian_coordinate_vec = np.array(cartesian_coordinate_vec) #NumPy arrayへの変換
    dim = len(cartesian_coordinate_vec) #入力ベクトルの次元数を取得
    r = np.linalg.norm(cartesian_coordinate_vec) #入力ベクトルのノルムを取得
        
        
    if r == 0: #上記1に該当
        if dim <= 1: #上記1-1に該当
            print("2次元以上のベクトルを入力してください") #2次元以上のベクトルしか極座標変換できない
        else: #上記1-2に該当
            return(np.zeros(shape=dim)) #引数と同じ大きさのゼロベクトルを返す 
    
    elif dim == 1: #上記2に該当
        print("2次元以上のベクトルを入力してください") #2次元以上のベクトルしか極座標変換できない
    
    elif dim == 2: #上記3に該当
        polar_coordinate_vec = np.array([]) #出力ベクトルを入れるベクトルを用意
        
        polar_coordinate_vec = np.append(polar_coordinate_vec, r) #出力ベクトルにノルムを格納
        polar_coordinate_vec = np.append(polar_coordinate_vec, atan2(cartesian_coordinate_vec[1], cartesian_coordinate_vec[0])) 
        
        return(polar_coordinate_vec)
    
    else: #上記4に該当
        polar_coordinate_vec = np.array([]) #出力ベクトルを入れるベクトルを用意
        polar_coordinate_vec = np.append(polar_coordinate_vec, atan2(cartesian_coordinate_vec[1], cartesian_coordinate_vec[0])) 
        
        for i in range(2, len(cartesian_coordinate_vec)):
            tmp = atan2(cartesian_coordinate_vec[i]*sin(polar_coordinate_vec[i-2]), cartesian_coordinate_vec[i-1])
            polar_coordinate_vec = np.append(polar_coordinate_vec, tmp)
        
        polar_coordinate_vec = np.append(polar_coordinate_vec, r)
        polar_coordinate_vec = polar_coordinate_vec[::-1]
        
        return(polar_coordinate_vec)

##実行例
実際に上記のメソッドを実行してみました。

test0 = [1, 1]
print(get_pc(test0))

このとき、表示結果は(1.414, 0.7853)となり、おおよそ$(\sqrt{2}, \pi/4)$となったので問題ないと判断しました。他の実行例も下表にまとめておきます。

引数	表示	返り値
[]	2次元以上のベクトルを入力してください	None
[0]	2次元以上のベクトルを入力してください	None
[1]	2次元以上のベクトルを入力してください	None
[0, 0]	-	[0, 0]
[1, 1]	-	[1.4142, 0.7853]
[1, 2, 3]	-	[3.7417, 0.9303, 1.1071]
[1, -1, 1, -1]	-	[ 2.000, 0.5236, -2.5261, -0.7854]

#まとめ
今回は多次元ベクトルに対して直交座標から極座標に変換するメソッドを作成しました。しかし、自作のメソッドは不安なので、どなたか3次元以上のベクトルに対しても直交座標から極座標に変換するメソッドをご存知の方がいらっしゃたら教えてください。

#参考文献
多次元の極座標変換