完成版「 対数らせん」

2023.03.02 「代数らせんアニメ」を元に、「対数らせんアニメ」を作りました。曲率 θ と曲率半径 r の関係は、r = a * exp( b * θ ) です。下のプログラムでは、a は 1 程度、b は 0.いくつか 程度でお試しください。b を小さくすると、らせん中心部の空洞が大きくなります。アニメとソースコードを示します。良ければ「正しい代数らせん」のコードと比較してみてください。

# spiral-logarithmic.py
# らせんを巻いていく。
# これは「対数らせん」(ベルヌーイのらせん)　　2023.03.02  by Kero

import sys
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation

# 描画画面設定。
fig, ax = plt.subplots(figsize=(16,4))
ax.set_aspect('equal') # 正方形の方眼紙。
artists = [] # パラパラ漫画設定。最初は白紙。

# 回転のマトリックス。
# 座標の列ベクトル群に左から掛けると図が原点中心に角度rdだけ右回転。
def A(rd):
    Am = np.array([[np.cos(rd), -np.sin(rd), 0],
                  [np.sin(rd),  np.cos(rd), 0],
                  [0., 0., 1.]])
    return Am
# x軸で反転させるマトリックス。
B = np.array([[1,0,0],[0,-1,0],[1,1,1]])

# 「対数らせん」関数（r= a・exp(b・θ)
# 原点に対して、回転角 ki、r= a・exp(b・θ） の点を生成。
def rasenB(a, b, ki):
    r = a * np.exp(b * ki)
    x = r * np.sin(ki)
    y = r * np.cos(ki)
    z = np.array([[x],[y],[1]]) # らせんにこの新しい点を追加。
    return z

# らせんを描く準備。
su = input("何回巻きますか？：")
aa = input("r= a・exp(b・θ)の定数 a : ")
bb = input("r= a・exp(b・θ)の定数 b : ")
xm = input("x軸長さ(空リターンで自動設定) : ")
ym = input("y軸長さ(空リターンで自動設定) : ")

suu = float(su) # 巻数の文字列から実数に変換。
tmk = 2 * np.pi * suu # トータルの巻角度。
kzm = tmk / 100 # トータル巻角度のキザミを100個とする。

b = float(bb)
a = float(aa) / np.exp(b*2*np.pi)

mak2 = 0. # 回転した長さ。最初はゼロ。
ras = np.array([[0],[0],[1]]) # 原点を中心に「原らせん」を生成。

title = "Logarithmic spiral (r= a*exp(b*θ) )"
title = title + ", turns= " + su
title = title + ", constant a= " + aa
title = title + ", constant b= " + bb

# らせんを描く。
for ki in np.linspace(0, tmk, 100):
    # 【対数らせん】= rasenB
    z = rasenB(a, b, ki)
    
    ras = np.append(ras, z, axis=1) # 原らせんに新たな点を追加。
    Am = A(ki) # アフィン変換の回転マトリックスを作る。角度ki右回転。
    ras2 = np.dot(Am, ras) # かいてーん！！
    z2 = ras2[:, -1:] # らせん終端（最外側）座標を取り出す。
    ras3 = ras2 - z2 # 終端を原点に平行移動。
    mak2 = mak2 + a * np.exp(b * ki) * kzm # 回転した長さを積立て貯金。
    mak3 = np.array([[mak2],[0],[0]]) # それをベクトルにして。
    ras4 = ras3 + mak3 # 回転した長さ分だけ右に平行移動。
    ras5 = np.dot(B,ras4) # x軸で反転。
    
    ras6 = np.delete(ras5, 0, axis=1) # らせんの最初の座標を削除。

    artist= ax.plot(ras6[0], ras6[1], c='b') # グラフ表示。
    artists.append(artist) # 一コマをパラパラ漫画の１ページに収録。

# xy軸の設定。入力値を調べる。数値以外なら自動設定。
if xm.isnumeric():
    xmax = float(xm)
else:
    xmax = np.max(ras5[0]) * 1.1
if ym.isnumeric():
    ymax = float(ym)
else:
    ymax = np.max(ras5[1]) * 1.1

plt.xlim(0,xmax) # x軸の長さは、最終らせんのx最大値の1割増し。
plt.ylim(0,ymax) # y軸の長さは、最終らせんのy最大値の1割増し。
plt.text(0, 1.1 * ymax, title)

# 以下はアニメの常套句。spiral-logarithmic.gif に保存。
anim = ArtistAnimation(fig, artists, interval=50, repeat=True)
anim.save('spiral-logarithmic.gif', writer='pillow')
        
plt.show()
plt.close()

sys.exit()