ちょこっとコーディング：幾何図形や曲線を描く

# はじめに
　突然ですが，数式や幾何にミステリアスな魅力を感じることはありませんか？
 　かっこいい幾何図形，きれいな曲線を描画したいとき，Web情報[^1][^2]やPython Matplotlibが頼りです．スキマ時間の「*ちょこっとコーディング*」，気楽なちょこっとプログラミングとしても手頃なお題です．図形や空間認識能力を心地よく刺激する効果も期待できます．コードが描くきれいな図形，それを作り出す数式の意味に自然に興味が出たら，とても良い感じです．某情報学部の各フロアにちりばめられた「あの曲線」も含まれます．

# いろいろな曲線
　図形の名前でどんな形状か，ビジュアルが頭の中に浮かんできますか？

- Lemniscate
- Heart Curve
- Cardioid
- Lissajous Curve
- Suddle Curve
- Rose Curve
- Golden Spiral
- Penrose Triangle
- Strange Attractor (Lorenz)

# 動作環境
- MacOS sonoma 14.1.2
- Python3.11.6
- pip 23.3.2
- VSCode 1.85.1
- JupyterLab 3.4.5
- numpy                         1.24.4
- japanize-matplotlib           1.1.3
- matplotlib                    3.5.3
- matplotlib-inline             0.1.6
- plotly                        5.13.0

必要に応じて，上記のライブラリをpipでインストールしておきましょう：
```
    $ python3 -m pip install update pip
    $ python3 -m pip install --upgrade pip
    $ python3 -m pip install numpy matplotlib plotly  
```

# 図形とコード
　以下，図形の視覚化例と，そのコード例を示します：
 
## Lemniscate
![Fig-Lemniscate.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635215/db553c5e-1ced-a3ef-fd9a-11217ada4f91.png)

```python
import numpy as np
import matplotlib.pyplot as plt

theta = np.linspace(0, 2 * np.pi, 1000)

def lemniscate(a,theta):
    x = a * np.sqrt(2)*np.cos(theta)/(1 + np.sin(theta)**2)
    y = x * np.sin(theta)
    return x , y

x,y = lemniscate(2,theta)

plt.plot(x, y,linewidth=12,color='gold')
plt.axis('equal')
plt.title('Lemniscate Curve')
plt.show()
```

## Heart Curve

![FIg-Heart.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635215/edb732a4-bf70-6e0c-6944-46deccc7b19c.png)

```python
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 2 * np.pi, 100)
x = 16 * np.sin(t) ** 3
y = 13 * np.cos(t) - 5 * np.cos(2 * t) - 2 * np.cos(3 * t) - np.cos(4 * t)

plt.plot(x, y, linewidth=16, color='gold')
plt.axis('equal')
plt.title('Heart Curve')
plt.show()
```

## Cardioid

![FIg-Cardioid.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635215/b8874a9d-335b-7775-991d-e1159e2a3d3a.png)

```python
import numpy as np
import matplotlib.pyplot as plt

theta = np.linspace(0,2*np.pi,300)
r =  3 + 5 * np.cos(theta)

plt.polar(theta, r, 'gold',linewidth=10)
plt.show()
```


## Lissajous

![FIg-Lissajous.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635215/2f897c40-515c-3586-f46c-4daf071793ae.png)

```python
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 2*np.pi, 300)
x = np.sin(t + np.pi/4)
y = np.cos(2*t)

theta = 12
theta = np.deg2rad(theta) 
R = np.array([[np.cos(theta), -np.sin(theta)], 
              [np.sin(theta), np.cos(theta)]])
x_kp, y_kp = np.dot(R, [x, y])

plt.figure(figsize=(4,4))
plt.plot(x_kp, y_kp, 'orange', linewidth=16)
plt.axis('equal')
plt.show()
```
![Fig-Kaishi.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635215/338d7300-6f99-788e-008c-163b559fbe00.png)

```python
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0, 2*np.pi, 600)
x = np.cos(2*t + np.pi/2)
y = np.sin(3*t)

theta = 0
theta = np.deg2rad(theta) 
R = np.array([[np.cos(theta), -np.sin(theta)], 
              [np.sin(theta), np.cos(theta)]])
x_kp, y_kp = np.dot(R, [x, y])

plt.figure(figsize=(4,4))
plt.plot(x_kp, y_kp, 'orange', linewidth=8)
plt.axis('equal')
plt.show()
```


## Suddle Point
![FIg-Suddle.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635215/8e4dff2b-1cd5-043d-e716-496e6b0d037c.png)


```python
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

X = np.linspace(-5, 5, 50)
Y = np.linspace(-5, 5, 50)
X, Y = np.meshgrid(X, Y)
Z = X**2 - Y**2

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot_surface(X, Y, Z, color='gold', cmap='copper')

plt.show()
```

## Suddle Point (plotly)

- Plotlyにより，対話的に視点を変えたり，拡大縮小回転しながら眺めることができます．（JupyterLab, Colaboratory経由でブラウザ表示）

![FIg-Suddle-plotly.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635215/35aa9d19-1c72-ec7d-4dd9-68626d3021a9.png)

```python
import numpy as np
import plotly.graph_objects as go

def f(x,y):
    return x**2 - y**2

x = np.linspace(-5, 5, 50) 
y = np.linspace(-5, 5, 50) 
X, Y = np.meshgrid(x, y) 

Z = f(X, Y) 

fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_layout(title='Saddle Point for x^2 - y^2') 
fig.show()
```

## Rose Curve
![FIg-Rose.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635215/36f34734-90a6-3a94-f7d0-16afebe08704.png)

```python
import numpy as np
import matplotlib.pyplot as plt

theta = np.linspace(0, 2*np.pi, 1000)

N = 5
r = 7 * np.cos(N * theta)

plt.polar(theta, r, 'gold',linewidth=7)
plt.show()
```

## Golden Spiral
![FIg-Spiral.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635215/967d95ed-290f-6e5a-da31-6567230e2e79.png)

```python
import matplotlib.pyplot as plt
import numpy as np

def polar_to_cartesian(r, theta):
    x = r * np.cos(theta)
    y = r * np.sin(theta)
    return x, y

theta_min = 0
theta_max = 6 * np.pi
theta_step = 0.02
theta = np.arange(theta_min, theta_max, theta_step)

phi = (1 + np.sqrt(5)) / 2  # golden ratio
r = phi ** (theta / np.pi)

x, y = polar_to_cartesian(r, theta)

plt.figure(figsize=(5,5))
plt.plot(x,y,'gold',linewidth=7)
plt.axis('equal')
plt.show()
```


## Penrose Triangle（Penrose三角形[^3]）

![FIg-Penrose.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635215/a7943d51-4226-d8a7-9b2d-365b269a38a5.png)
```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import math
# 重心（中心）に対して，一定角度で回転した座標
def rotate_th_degrees(th,X,Y):
    N = len(X)
    def center(x,y):
        xs,ys = sum(x),sum(y)
        return xs/N, ys/N
    rth = math.radians(th)
    
    Rx,Ry = [],[]
    G = center(X,Y)
    for i in range(N):
        x = X[i] - G[0]
        y = Y[i] - G[1]
        xt = x * math.cos(rth) - y * math.sin(rth)
        yt = x * math.sin(rth) + y * math.cos(rth)
        Rx += [xt + G[0]]
        Ry += [yt + G[1]]
    return Rx,Ry
    
# Penroseの不可能三角形
def penrose():
    k = math.sqrt(3)/2
    S = 2.0
    L = S/8
    H = L * k
    x0 = S/2
    y0 = k*x0
    def p(x,y):
        return list( [ L*x-x0, H*y-y0 ] )
    p1 = p(1.0,0.0)
    p2 = p(7.0,0.0)
    p3 = p(7.5,1.0)
    p4 = p(4.5,7.0)
    p5 = p(3.5,7.0)
    p6 = p(0.5,1.0)
    p7 = p(4.5,5.0)
    p8 = p(3.0,2.0)
    p9 = p(2.0,2.0)
    p10 = p(5.0,2.0)
    p11 = p(5.5,1.0)
    p12 = p(4.0,4.0)
    plist = [p1,p2,p3,p4,p5,p6,p1,p2,p7,p8,p9,p4,p9,p10,p12,p11,p6] 
    return [ p[0] for p in plist ],[ p[1] for p in plist ]

x,y = penrose()
px,py = rotate_th_degrees(180,x,y)

fig = plt.figure()
fig, ax = plt.subplots()
ax.plot(x,y,linewidth=4)
ax.plot(px,py,linewidth=4)
ax.set_aspect('equal')
plt.show()
```

- Penrose三角形を，重心に対して45度ずつ，回転表示
![FIg-R-Penrose.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635215/a69d0ec2-46d2-b3e4-9f99-6f0a154cfe73.png)

```python
fig = plt.figure()
fig, ax = plt.subplots()

x,y = penrose()
ax.plot(x,y,linewidth=2)
ax.set_aspect('equal')    

for th in range(45,360,45):
    x,y = rotate_th_degrees(th,x,y)
    ax.plot(x,y,linewidth=2)

plt.show()
```


## Strange Attractor (Lorenz)[^4]
![FIg-Lorenz.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3635215/0ded06c9-44c8-805f-c81f-a5a9a6dd1239.png)

```Python
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import colors

def lorenz(x, y, z, s=8, r=22, b=8/3):
    xd = s*(y - x)
    yd = r*x - y - x*z
    zd = x*y - b*z
    return xd, yd, zd

steps,dt = 18000,0.005
xs = np.empty(steps+1)
ys = np.empty(steps+1)
zs = np.empty(steps+1)
xs[0], ys[0], zs[0] = 1, 1, 1

for i in range(steps):
    xd, yd, zd = lorenz(xs[i], ys[i], zs[i])
    xs[i+1] = xs[i] + xd*dt
    ys[i+1] = ys[i] + yd*dt
    zs[i+1] = zs[i] + zd*dt

ax = plt.figure(figsize=(8,8), dpi=300).add_subplot(projection='3d')
ax.view_init(elev=25, azim=120)
cn = colors.Normalize(ys.min(),ys.max())

for i in range(steps):
    ax.plot(xs[i:i+2], ys[i:i+2], zs[i:i+2], 
            color=plt.cm.viridis(cn(xs[i])), lw=0.3)

ax.set_title('Lorenz "strange" Attractor', size=8)
plt.show()
```

# おわりに

　幾何や数式に興味を抱くきっかけと，スキマ時間に*ちょこっとコーディング*のおすすめ部材として，きれいな曲線，幾何図形を描いてみるのも，ちょっとおしゃれな趣味ではないでしょうか？


#  参考情報：
- いろいろな曲線の確認 (http://izumi-math.jp/S_Yoshida/matome/sc_iroironakyokusen.pdf)
- 媒介変数表示された有名な曲線７つ(https://manabitimes.jp/math/898)
- Penroseの不可能三角形　(https://ja.wikipedia.org/wiki/ペンローズの三角形)
- ローレンツ・アトラクタ（https://ja.wikipedia.org/wiki/アトラクター）


[^1]: いろいろな曲線の確認 (http://izumi-math.jp/S_Yoshida/matome/sc_iroironakyokusen.pdf)
[^2]: 媒介変数表示された有名な曲線７つ(https://manabitimes.jp/math/898)
[^3]: Penroseの不可能三角形　(https://ja.wikipedia.org/wiki/ペンローズの三角形)
[^4]: アトラクタ（https://ja.wikipedia.org/wiki/アトラクター）