はじめに

　Matplotlibライブラリを利用して、円のグラフを作成します。

【目次】

はじめに
円周の作図
おわりに

円周の作図

　Matplotlibライブラリを利用して、2次元空間(平面)上に円周(circumference)のグラフを作成します。また、円座標系(circular coordinates)をグラフで確認します。

　利用するライブラリを読み込みます。

# 利用ライブラリ
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

座標の計算

　まずは、円周上の点の座標計算を確認します。

　角度と半径を指定して、円周上の点(ノルムが指定した半径になる2次元ベクトル)を計算します。

# 半径を指定
r = 14.0

# 度数法の角度を指定
t = 180.0

# 弧度法の角度(ラジアン)を計算
t *= np.pi / 180.0
print(t)

# 円周上の点の座標を計算
x = r * np.cos(t)
y = r * np.sin(t)

# 座標を格納
p = np.array([x, y])
print(p)
print(np.linalg.norm(p))

3.141592653589793
[-1.40000000e+01  1.71450552e-15]
14.0

　度数法における角度を $\theta^{\circ}$ 、弧度法における角度(ラジアン)を $\theta$ で表すことにします。度数法と弧度法は $\theta = \frac{2 \pi}{360^{\circ}} \theta^{\circ}$ 、 $\theta^{\circ} = \frac{360}{2 \pi} \theta$ の関係です。つまり、 $\theta^{\circ} = 0^{\circ}$ のとき $\theta = 0$ 、 $\theta^{\circ} = 180^{\circ}$ のとき $\theta = \pi$ 、 $\theta^{\circ} = 360^{\circ}$ のとき $\theta = 2 \pi$ です。
　 $\pi$ は円周率で、np.piで扱えます。

　半径を $r$ 、角度(x軸とのなす角)を $\theta$ とします。オブジェクト名はr, tとします。
　半径が $r$ の円周上の点(ノルムが $r$ の2次元ベクトル)を $\mathbf{p} = (x, y)$ とすると、x軸・y軸の座標は、次の式で計算できます。

$\displaystyle \begin{cases} x = r \cos \theta \\ y = r \sin \theta \end{cases} \tag{1}$

　直交座標系の座標を $(x, y)$ 、円座標系の座標を $(r, \theta)$ で表します。
　式(1)は、円座標系から直交座標系への変換式です。

　 $\mathbf{p}$ のユークリッドノルム $\|\mathbf{p}\| = \sqrt{x^2 + y^2}$ は $r$ になります。ユークリッドノルムは、np.linalg.norm()で計算できます。

　あるいは、2次元ベクトルを指定して、球面上の点(ノルムが指定した半径になるベクトル)に変換します。

# ベクトルを指定
p = np.array([10.0, 17.0])
#p = np.array([x, y]) # 確認用
print(np.linalg.norm(p))

# 半径を指定
r = 14.0

# ノルムをrに変換(円周上の座標を計算)
p *= r / np.linalg.norm(p)
print(p)
print(np.linalg.norm(p))

19.72308292331602
[ 7.09828177 12.06707901]
14.0

　任意の2次元ベクトル $\tilde{\mathbf{p}} = (\tilde{x}, \tilde{y})$ を指定して、ノルムで割り正規化(ノルムを1に変換)します。さらに、半径 $r$ を掛けるとノルムが $r$ のベクトルになります。

$\displaystyle \mathbf{p} = r \frac{ \tilde{\mathbf{p}} }{ \|\tilde{\mathbf{p}}\| } = r \frac{ \tilde{\mathbf{p}} }{ \sqrt{\tilde{x}^2 + \tilde{y}^2} }$

　この計算により、2次元空間上の点 $\tilde{\mathbf{p}}$ と角度 $\theta$ が同じである球面上の点(ベクトル $\tilde{\mathbf{p}}$ と同じ方向であるノルムが $r$ のベクトル) $\mathbf{p}$ が得られます。ベクトルのスカラ倍なのでベクトルの向きは変わりません。

　円周上の点から半径と角度を計算します。

# 半径を計算
r = np.linalg.norm(p)
print(r)

# 弧度法の角度(ラジアン)を計算
sgn_y = 1.0 if p[1] >= 0.0 else -1.0
t = sgn_y * np.arccos(p[0] / np.linalg.norm(p))
print(t)

# 度数法の角度を計算
t *= 180.0 / np.pi
print(t)

14.0
1.0390722595360908
59.53445508054011

　半径 $r$ と角度 $\theta$ は、次の式で計算できます。

$\displaystyle \begin{align*} &\begin{cases} r = \|\mathbf{p}\| = \sqrt{x^2 + y^2} \\ \theta = \mathrm{sgn}(y) \arccos \Bigl( \frac{z}{\|\mathbf{p}\|} \Bigr) = \mathrm{sgn}(y) \arccos \Bigl( \frac{x}{\sqrt{x^2 + y^2}} \Bigr) \\ \end{cases} \tag{2}\\ &\mathrm{sgn}(y) = \begin{cases} 1 & (y \geq 0) \\ -1 & (y < 0) \end{cases} \end{align*}$

　 $\arccos x$ は逆コサイン関数(コサイン関数 $\cos x$ の逆関数)で、np.arccos()で計算できます。y軸の値 $y$ の符号を $\mathrm{sgn}(y)$ で表しています。
　式(2)は、直交座標系から円座標系への変換式です。

　ここでは、2通りの座標計算方法を確認しました。以降は、1つ目の方法を用います。

円の描画

　次は、円周全体の座標を計算して、円のグラフを作成します。

　1周期分のラジアンを作成して、円周上の座標を計算します。周期については「変数と各軸の値の関係」で確認します。

# 半径を指定
r = 1.0

# ラジアンを作成
t = np.linspace(start=0.0, stop=2.0*np.pi, num=81)
print(t[:5])
print(len(t))

# 円周の座標を計算
x = r * np.cos(t)
y = r * np.sin(t)
print(x[:5].round(3))
print(y[:5].round(3))

[0.         0.07853982 0.15707963 0.23561945 0.31415927]
81
[1.    0.997 0.988 0.972 0.951]
[0.    0.078 0.156 0.233 0.309]

　 $0 \leq \theta \leq 2 \pi$ の範囲のラジアンを作成してtとします。
　 $\theta$ の範囲(最大値と最小値の差)が $2 \pi$ であれば、他の値であっても同じグラフになります。
　tを用いて、座標の計算式(1)を計算して、x軸・y軸の値をそれぞれx, yとします。

　確認のため、点ごとにノルムを計算します。また、色付け用の配列を作成します。

# 2つの軸の値を結合
xy = np.stack([x, y], axis=1)
print(xy[:5])
print(xy.shape)
print(np.linalg.norm(xy, axis=1)[:5])

# 2つの軸の値と0ベクトルを結合
xyo = (np.stack([x, y, np.zeros_like(x)], axis=1) + r) * 0.5 / r
print(xyo[:5])

[[1.         0.        ]
 [0.99691733 0.0784591 ]
 [0.98768834 0.15643447]
 [0.97236992 0.23344536]
 [0.95105652 0.30901699]]
(81, 2)
[1. 1. 1. 1. 1.]
[[1.         0.5        0.5       ]
 [0.99845867 0.53922955 0.5       ]
 [0.99384417 0.57821723 0.5       ]
 [0.98618496 0.61672268 0.5       ]
 [0.97552826 0.6545085  0.5       ]]

　ノルムの計算用に、x, yをそれぞれ1列として結合した配列xyを作成します。
　np.linalg.norm()のaxix引数に1を指定して行(点)ごとにノルムを計算すると、どの点も値がrになるのを確認できます。

　散布図の装飾に関して、各点の座標 $(x, y)$ を正規化した値をRGB値として使うことで、点ごとの座標に応じて色付けられます。
　x, yと全てが0の列を結合した配列xyoを作成します。どの列(軸)を0の列にするかによって、配色が変わります。
　各軸の値は最小値が $-r$ 、最大値が $r$ です。そこで、各点に $r$ を足すことで $0$ から $2 r$ の値になり、さらに $2 r$ で割ることで $0$ から $1$ の値に変換できます。x軸(1つ目の列)の値が大きい( $x$ が $r$ に近く $y$ が $0$ に近い)ほど赤色、y軸(2つ目の列)の値が大きいほど緑色、3つ目の列の値が大きいほど青色に近い色になります。

　円のグラフを作成します。

# タイトル用の文字列(座標の計算式)を作成
coord_label = '$x = r\ \\cos \\theta$\n'
coord_label += '$y = r\ \\sin \\theta$'

# 円を作図
fig, ax = plt.subplots(figsize=(7, 7), facecolor='white')
ax.plot(x, y, color='black') # 曲線
#ax.scatter(x, y, c=y) # 散布図:(y軸の値により色付け)
#ax.scatter(x, y, color=xyo) # 散布図:(xy座標の値により色付け)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('$r=' + str(r) + '$\n' + coord_label, loc='left')
fig.suptitle('circle', fontsize=20)
ax.grid()
ax.set_aspect('equal')
plt.show()

　Axes.plot()で曲線(折れ線グラフ)、Axes.scatter()で散布図として円を描画できます。

　綺麗な円を描画するには、Axes.set_aspect('equal')を使ってアスペクト比を1に設定します。

　cmap引数にカラーマップ名を指定して、y軸の値などに応じてグラデーションで色付けられます。
　散布図に関して、color引数に「各点の座標 $(x, y)$ を正規化した値をRGB値」として指定して、点ごとの座標に応じて色付けられます。

　円は、原点からの距離(ノルム)が $r$ の点なのが分かります。

　以上で、円を描画できました。次からは、変数(角度)と2つの軸の値の関係を確認します。

変数と各軸の値の関係

　角度とx軸・y軸それぞれの値の関係をグラフで確認します。

　1周期分以上のラジアンを作成して、円周の座標を計算します。

# 半径を指定
r = 1.0

# ラジアンを作成
t = np.linspace(start=-2.0*np.pi, stop=2.0*np.pi, num=200)
print(t[:5])
print(len(t))

# 円周の座標を計算
x = r * np.cos(t)
y = r * np.sin(t)
print(x[:5].round(3))
print(y[:5].round(3))

[-6.28318531 -6.22003772 -6.15689013 -6.09374253 -6.03059494]
200
[1.    0.998 0.992 0.982 0.968]
[0.    0.063 0.126 0.188 0.25 ]

　「円の描画」のときと同様にして、円周の座標を計算します。ただしこの例では、周期性が分かるように、ラジアン $\theta, \phi$ の範囲をそれぞれ $-2 \pi$ から $2 \pi$ とします。

　 $\theta$ 軸のラジアン目盛用の配列を作成します。グラフ自体には不要な処理です。

# 半円の目盛の数(分母の値)を指定
denom = 6.0

# 目盛の通し番号(分子の値)を作成
tick_vals = np.arange(start=-2.0*denom, stop=2.0*denom+0.1)
print(tick_vals[:5])

# 目盛用のラジアンを計算
tick_t_vals = tick_vals / denom * np.pi
print(tick_t_vals[:5].round(2))

# 目盛ラベルを作成
tick_label_vals = [
    '$' + ('-' if nomer < 0.0 else '') + '\\frac{' + str(abs(nomer)) + '}{' + str(int(denom)) + '} \pi$' for nomer in tick_vals.astype('int16')
]
print(tick_label_vals[:3])

[-12. -11. -10.  -9.  -8.]
[-6.28 -5.76 -5.24 -4.71 -4.19]
['$-\\frac{12}{6} \\pi$', '$-\\frac{11}{6} \\pi$', '$-\\frac{10}{6} \\pi$']

　軸目盛ラベルを $\frac{i}{n} \pi$ の形で表示することにします。この例では2周期分を描画するので、 $i$ を $-2 n$ から $2 n$ の整数 $i = -2 n, -2 (n-1), \dots , -1, 0, 1, \dots, 2 (n-1), 2 n$ とします。 $n$ をdenom、 $i$ をtick_valsとして値を指定します。
　軸目盛の座標用のラジアン $t = \frac{i}{n} \pi$ を計算してtick_t_valsとします。

　 $\theta$ と $x$ の関係のグラフを作成します。

# x軸の値を曲線で作図
fig, ax = plt.subplots(figsize=(9, 3), facecolor='white')
ax.plot(t, x)
ax.set_xticks(ticks=tick_t_vals, labels=tick_label_vals) # θ目盛ラベル
ax.set_xlabel('$\\theta$')
ax.set_ylabel('$x$')
ax.set_title('$r=' + str(r) + '$', loc='left')
fig.suptitle('$x = r\ \\cos \\theta$', fontsize=20)
ax.grid()
ax.set_aspect('equal')
plt.show()

　 $\theta = t, t + 2 \pi$ の $x$ の値が一致するのが分かります。 $2 \pi$ ごとに周期的に推移します。

　 $\theta$ と $y$ の関係のグラフを作成します。

# y軸の値を曲線で作図
fig, ax = plt.subplots(figsize=(9, 3), facecolor='white')
ax.plot(t, y)
ax.set_xticks(ticks=tick_t_vals, labels=tick_label_vals) # θ目盛ラベル
ax.set_xlabel('$\\theta$')
ax.set_ylabel('$y$')
ax.set_title('$r=' + str(r) + '$', loc='left')
fig.suptitle('$y = r\ \\sin \\theta$', fontsize=20)
ax.grid()
ax.set_aspect('equal')
plt.show()

　こちらも $2 \pi$ ごとに周期します。

　ここでは、変数(角度)と2つの軸それぞれの値の関係を確認しました。次は、全ての軸の値との関係を確認します。

変数と座標の関係

　最後に、円周上に変数(角度)に関する軸目盛を描画して、変数と座標の関係をグラフで確認します。

　円周の座標を計算します。

# 半径を指定
r = 1.0

# ラジアンを作成
t = np.linspace(start=0.0, stop=2.0*np.pi, num=100)
print(t[:5])
print(len(t))

# 円周の座標を計算
x = r * np.cos(t)
y = r * np.sin(t)

[0.         0.06346652 0.12693304 0.19039955 0.25386607]
100

　「円の描画」のときのコードで、円周の座標を計算します。

　 $\theta$ 軸目盛用の配列を作成します。

# 半円の目盛の数(分母の値)を指定
denom = 6.0

# 目盛の通し番号(分子の値)を作成
tick_vals = np.arange(start=0.0*denom, stop=2.0*denom)
print(tick_vals)

# 目盛用のラジアンを計算
tick_t_vals = tick_vals / denom * np.pi
print(tick_t_vals.round(2))

# 目盛ラベルを作成
tick_label_vals = [
    '$' + ('-' if nomer < 0.0 else '') + '\\frac{' + str(abs(nomer)) + '}{' + str(int(denom)) + '} \pi$' for nomer in tick_vals.astype('int16')
]
print(tick_label_vals[:3])

# 目盛用の座標を計算
tick_x_vals = r * np.cos(tick_t_vals)
tick_y_vals = r * np.sin(tick_t_vals)
print(tick_x_vals.round(2))
print(tick_y_vals.round(2))

[ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11.]
[0.   0.52 1.05 1.57 2.09 2.62 3.14 3.67 4.19 4.71 5.24 5.76]
['$\\frac{0}{6} \\pi$', '$\\frac{1}{6} \\pi$', '$\\frac{2}{6} \\pi$']
[ 1.    0.87  0.5   0.   -0.5  -0.87 -1.   -0.87 -0.5  -0.    0.5   0.87]
[ 0.    0.5   0.87  1.    0.87  0.5   0.   -0.5  -0.87 -1.   -0.87 -0.5 ]

　「変数と各軸の値の関係」のときと同様にして、ラジアン目盛用の値を作成します。この例では1周期分を描画するので、 $i$ を $0$ から $2 n$ の整数 $i = 0, 1, \dots, 2 n$ とします。
　さらに、軸目盛用のラジアンtick_t_valsを用いて、目盛の描画位置の座標を計算してtick_*_valsとします。

　なす角 $\theta$ の角マークと角ラベル用の配列を作成します。

# 角マーク用のラジアンを作成
angle_t_vals = np.linspace(start=0.0, stop=1.0/denom*np.pi, num=50)

# 角マーク用の座標を計算
d = 0.1
angle_x_vals = d * np.cos(angle_t_vals)
angle_y_vals = d * np.sin(angle_t_vals)

# 角ラベル用の座標を計算
d = 0.15
angle_x = d * np.cos(0.5/denom*np.pi)
angle_y = d * np.sin(0.5/denom*np.pi)

　この例では、なす角 $\theta$ を示すために、軸目盛1つ分の角マークを表示することにします。
　そこで、角マーク用のラジアン $0 \leq t \leq \frac{\pi}{n}$ をangle_t_valsとして、円弧の座標を計算してangle_***_valsとします。

　角ラベルは、角マークの中点に表示することにします。
　そこで、 $r \sin \frac{\pi}{2 n}$ と $r \cos \frac{\pi}{2 n}$ を計算してangle_***とします。

　円周上に変数(角度)軸を重ねて表示したグラフを作成します。

# タイトル用の文字列(座標の計算式)を作成
coord_label = '$x = r\ \\cos \\theta$\n'
coord_label += '$y = r\ \\sin \\theta$'

# グラフサイズ用の値を指定
axis_size = r * 1.5

# 円周上の変数(ラジアン)目盛を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white')
ax.plot(x, y, color='black') # 円
ax.plot([[0, 0], [1, 0]], [[0, 0], [0, 1]],
        color='black', linewidth=2) # x,y軸
ax.plot(np.stack([np.zeros_like(tick_vals), tick_x_vals], axis=1).T, 
        np.stack([np.zeros_like(tick_vals), tick_y_vals], axis=1).T, 
        color='gray', linestyle=':') # θ目盛グリッド線
for i in range(len(tick_vals)):
    ax.annotate(xy=[tick_x_vals[i], tick_y_vals[i]], 
                text='|', rotation=tick_vals[i]*30.0+90.0, ha='center', va='center') # θ目盛線
    d = 1.15
    ax.text(tick_x_vals[i]*d, tick_y_vals[i]*d, 
            s='$\\frac{' + str(int(tick_vals[i])) + '}{' + str(int(denom)) + '} \pi$', 
            size=15, ha='center', va='center') # θ目盛ラベル
ax.plot(angle_x_vals, angle_y_vals, 
        color='black', lw=1) # θの角マーク
ax.text(angle_x, angle_y, 
        s='$\\theta$', size=15, ha='center', va='center') # θの角ラベル
ax.text(0.5*r, 0.0, 
        s='$r$', size=15, ha='center', va='bottom') # 半径rラベル
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('$r=' + str(r) + '$\n' + coord_label, loc='left')
fig.suptitle('circular coordinates', fontsize=20)
ax.grid()
ax.set_xlim(left=-axis_size, right=axis_size)
ax.set_ylim(bottom=-axis_size, top=axis_size)
ax.set_aspect('equal')
plt.show()

　目盛ラベルや角ラベルをAxes.text()で描画します。
　x軸・y軸・z軸を黒色の線分として描画します。

　x軸線から反時計回りの角度(なす角)が $\theta$ に対応しているのが分かります。時計回りの角度は負の値で表現できます。

　円周上の点を変化させたアニメーションで確認します。

・作図コード(クリックで展開)

# フレーム数を指定
frame_num = 120

# ラジアンとして利用する値を作成
t_n = np.linspace(start=-2.0*np.pi, stop=2.0*np.pi, num=frame_num+1)[:frame_num]

# グラフオブジェクトを初期化
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white')
fig.suptitle('circular coordinates', fontsize=20)

# 作図処理を関数として定義
def update(i):
    
    # 前フレームのグラフを初期化
    plt.cla()

    # i番目のラジアンを取得
    t = t_n[i]
    
    # 円周上の点の座標を計算
    point_x = r * np.cos(t)
    point_y = r * np.sin(t)
    
    # 角マーク用のラジアンを作成
    angle_t_vals = np.linspace(start=0.0, stop=t, num=100)

    # 角マーク用の座標を計算
    d = 0.1
    angle_x_vals = d * np.cos(angle_t_vals)
    angle_y_vals = d * np.sin(angle_t_vals)

    # 角ラベル用の座標を計算
    d = 0.15
    angle_x = d * np.cos(0.5*t)
    angle_y = d * np.sin(0.5*t)
    
    # タイトル用の文字列(座標の計算式)を作成
    coord_label = '$r=' + str(r) + '$\n'
    coord_label += '$\\theta^{\circ}='+str((t*180.0/np.pi).round(1))+'^{\circ}, ' + '\\theta='+str((t/np.pi).round(2))+'\pi$\n'
    coord_label += '$x = r\ \\cos \\theta, y = r\ \\sin \\theta$'
    
    # 円周上の点を作図
    ax.plot(x, y, color='black') # 円
    ax.plot([[0.0, 0.0], [1.0, 0.0]], [[0.0, 0.0], [0.0, 1.0]],
            color='black', linewidth=2) # x,y軸
    ax.plot(np.stack([np.zeros_like(tick_vals), tick_x_vals], axis=1).T, 
            np.stack([np.zeros_like(tick_vals), tick_y_vals], axis=1).T, 
            color='gray', linestyle=':') # θ目盛グリッド線
    for i in range(len(tick_vals)):
        ax.annotate(xy=[tick_x_vals[i], tick_y_vals[i]], 
                    text='|', rotation=tick_vals[i]*30.0+90.0, ha='center', va='center') # θ目盛線
        d = 1.15
        ax.text(tick_x_vals[i]*d, tick_y_vals[i]*d, 
                s='$\\frac{' + str(int(tick_vals[i])) + '}{' + str(int(denom)) + '} \pi$', 
                size=15, ha='center', va='center') # θ目盛ラベル
    ax.scatter(point_x, point_y, 
               s=50, 
               label='('+str(point_x.round(2))+', '+str(point_y.round(2))+')', zorder=10) # 円周上の点
    ax.plot([0.0, point_x], [0.0, point_y]) # 原点と円周上の点の線分
    ax.plot(angle_x_vals, angle_y_vals, 
            color='black', lw=1) # θの角マーク
    ax.text(angle_x, angle_y, 
            s='$\\theta$', size=15, ha='center', va='center') # θの角ラベル
    ax.text(0.5*r, 0.0, 
            s='$r$', size=15, ha='center', va='bottom') # 半径rラベル
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_title(coord_label, loc='left')
    ax.grid()
    ax.legend(loc='upper left')
    ax.set_xlim(left=-axis_size, right=axis_size)
    ax.set_ylim(bottom=-axis_size, top=axis_size)
    ax.set_aspect('equal')

# gif画像を作成
ani = FuncAnimation(fig=fig, func=update, frames=frame_num, interval=100)

# gif画像を保存
ani.save('circle_axis_2d.gif')