三次元直交座標と極座標との変換〜pythonでの実装〜

今回は様々な分野で用いる三次元の直交座標(デカルト座標)と極座標の変換を概説し，その後pythonによる実装を試みます。
現在書き進めている音響理論基礎でも後々利用することになりますので，もしそちらをご覧頂いている方も目を通して頂ければ良いかと思います。

では早速本題です。

直交座標と極座標の座標変換

まず今回考える状況は次の図のようになっています。

二次元極座標では方位角が$\theta$なのですが三次元では方位角が$\phi$で仰角の方が$\theta$になっているので注意です。

これについて，「自分は方位角が$\theta$の方が慣れてるから」と言って入れ替えると他の文献やサイトとの整合性を取るのが少し面倒になりますので，個人的にはマジョリティに合わせて覚えるのが無難だと考えています。

この場合の極座標$(r, \theta, \phi)$から直交座標$(x, y, z)$への座標変換公式は次のようになります。

$$\begin{cases}
{x= r \sin \theta \cos \phi} \\
{y=r \sin \theta \sin \phi} \\
{z=r \cos \theta}
\end{cases}$$

直交座標$(x, y, z)$から極座標$(r, \theta, \phi)$への変換は次のとおりです。

$$\begin{cases}
{r= \sqrt{x^{2} + y^{2} + z^{2}}} \\
{\theta = \cos^{-1} \frac{z}{\sqrt{x^{2} + y^{2} + z^{2}}}} \\
{\phi=\tan^{-1} \frac{y}{x}}
\end{cases}$$

基底ベクトルの表現とベクトル場の変換

極座標での基底ベクトル$(\boldsymbol{e}_{r}, \boldsymbol{e}_{\theta}, \boldsymbol{e}_{\phi})$は直交座標の基底ベクトル$(\boldsymbol{e}_{x}, \boldsymbol{e}_{y}, \boldsymbol{e}_{z})$を用いて次のように表現(変換)されます。

$$\begin{cases}
{\boldsymbol{e}_{r} = \sin \theta \cos \phi \boldsymbol{e}_{x} + \sin \theta \sin \phi \boldsymbol{e}_{y} + \cos \theta \boldsymbol{e}_{z}} \\
{\boldsymbol{e}_{\theta} = \cos \theta \cos \phi \boldsymbol{e}_{x} + \cos \theta \sin \phi \boldsymbol{e}_{y} – \sin \theta \boldsymbol{e}_{z}} \\
{\boldsymbol{e}_{\phi} = – \sin \phi \boldsymbol{e}_{x} + \cos \phi \boldsymbol{e}_{y}}
\end{cases}$$

若しくはよく以下のように表記されます。

\begin{align}
\boldsymbol{e}_{r} &= \left ( \begin{array}{l} {\sin \theta \cos \phi} \ {\sin \theta \sin \phi} \ {\cos \theta}\end{array} \right ), \notag \\
\boldsymbol{e}_{\theta} &= \left ( \begin{array}{l} {\cos \theta \cos \phi} \ {\cos \theta \sin \phi} \ {- \sin \theta}\end{array} \right ), \notag \\
\boldsymbol{e}_{\phi} &= \left ( \begin{array}{l} {- \sin \phi} \ {\cos \phi} \ {0}\end{array} \right ) \notag \\
\end{align}

これを用いて極座標の基底で表されるベクトル場$\boldsymbol{W} = W_{r} \boldsymbol{e}_{r} + W_{\theta} \boldsymbol{e}_{\theta} + W_{\phi} \boldsymbol{e}_{\phi}$を直交座標の基底で表されるベクトル場に変換してみましょう。変換先のベクトル場を$\boldsymbol{V} \equiv V_{x} \boldsymbol{e}_{x} + V_{y} \boldsymbol{e}_{y} + V_{z} \boldsymbol{e}_{z}$とすると次のように書けます。

$$\begin{cases}
{V_{x} = \sin \theta \cos \phi W_{r} + \cos \theta \cos \phi W_{\theta} – \sin \phi W_{\phi}} \\
{V_{y} = \sin \theta \sin \phi W_{r} + \cos \theta \sin \phi W_{\theta} + \cos \phi W_{\phi}} \\
{V_{z} = \cos \theta W_{r} – \sin \theta W_{\theta}}
\end{cases}$$

このようにベクトル場も簡単に変換できます。※直交座標から極座標へのベクトル場の変換はあまり利用しないと思い省いています

pythonによる実装

ここまでの内容をpythonの関数で実装してみましょう。
まずは極座標$(r, \theta, \phi)$から直交座標$(x, y, z)$への座標変換です。
numpyをnpとしてインポートしてご利用ください。

# polar coordinates -&gt; cartesian coordinates
def polar2cartesian(position):
    num = position.shape[0]
    newPosition = np.empty([num,3], dtype=np.float64)
    newPosition[:,0] = position[:,0] * np.sin(position[:,1]) * np.cos(position[:,2])
    newPosition[:,1] = position[:,0] * np.sin(position[:,1]) * np.sin(position[:,2])
    newPosition[:,2] = position[:,0] * np.cos(position[:,1])
    return newPosition

式そのままですね。

次に直交座標$(x, y, z)$から極座標$(r, \theta, \phi)$への変換です。

# cartesian coordinates -&gt; polar coordinates
def cartesian2polar(position):
    num = position.shape[0]
    newPosition = np.empty([num,3], dtype=np.float64)
    newPosition[:,0] = np.linalg.norm(position, axis=1)
    newPosition[:,1] = np.arccos(position[:,2] / newPosition[:,0])
    newPosition[:,2] = np.arctan2(position[:,1], position[:,0])
    nan_index = np.isnan(newPosition[:,1])
    newPosition[nan_index,1] = 0
    return newPosition

ここで0割りによってnanが発生するインデックスについては0に置き換えるようにしています。
最後に極座標のベクトル場$\boldsymbol{W} = W_{r} \boldsymbol{e}_{r} + W_{\theta} \boldsymbol{e}_{\theta} + W_{\phi} \boldsymbol{e}_{\phi}$を直交座標のベクトル場$\boldsymbol{V} \equiv V{x} \boldsymbol{e}_{x} + V_{y} \boldsymbol{e}_{y} + V_{z} \boldsymbol{e}_{z}$に変換するプログラムです。

# cartesian coordinates -&gt; polar coordinates
# vector field of polar coordinates -&gt; vector field of cartesian coordinates
def polar2cartesian_vec(W, theta, phi):
    V = np.empty_like(W)
    V[0] = W[0]*np.sin(theta)*np.cos(phi) + W[1]*np.cos(theta)*np.cos(phi) - W[2]*np.sin(phi)
    V[1] = W[0]*np.sin(theta)*np.sin(phi) + W[1]*np.cos(theta)*np.sin(phi) + W[2]*np.cos(phi)
    V[2] = W[0]*np.cos(theta) - W[1]*np.sin(theta)
    return V

実装は簡単ですね。