[python] 流線描画

#流線
直交座標系$\boldsymbol{x} = (x, y)$，速度$\boldsymbol{u} = (u, v)$で表される流れ場の，ある時刻$t$における<b>流線（streamline）</b>は数学的に常微分方程式
$$
\frac{\text{d} x}{\text{d} u(x,y,t)} = \frac{\text{d} y}{\text{d} v(x,y,t)}
$$
で定義されます．この記事で紹介する流線はこの常微分方程式を解いて描くわけではなく，単純に速度ベクトルをつないだものを描画しているものと思ってください．

#matplotlib組み込みのstreamplotで描く

例として2次元定常Taylor-Green渦：
$$
\begin{aligned}
u &= \sin x \cos y \\\\
v &= -\cos x \sin y
\end{aligned}
$$
を使います．
matplotlib組み込みの[```streamplot```](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.streamplot.html)で簡単に流線を描くことができます．

```python

import matplotlib.pyplot as plt
import numpy as np
from numpy import sin, cos, pi, sqrt

n = 1000
x = np.linspace(0, 2*pi, n)
y = np.linspace(0, 2*pi, n)
X, Y = np.meshgrid(x, y)

# Taylor-Green Vortices
u = sin(X)*cos(Y)
v = - cos(X)*sin(Y)
speed = sqrt(u**2 + v**2)

# Plot
plt.figure(1)
plt.streamplot(X, Y, u, v, density=3, color='k', arrowstyle='-', linewidth=0.6)
plt.contourf(X, Y, speed, 100, cmap='viridis')
plt.show()
```

流線が描けました．背景は流速です．
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/119868/612ce3d6-06b1-8ced-00c8-44027b2e9c22.png" width=90%>


#Line Integral Convolution (LIC) で描く
こっちが本題．上のmatplotlib組み込みの```streamplot```でも十分綺麗なのですが，もっと美しく（芸術的に）描きたいときは<b>Line Integral Convolution (LIC，線積分畳み込み)</b> という手法を用いると流線をとても美しく描画できます．[wikipedia（英語）](https://en.wikipedia.org/wiki/Line_integral_convolution)に簡単な説明があります．アルゴリズムは[Cabral and Leedom, 1993](https://dl.acm.org/doi/pdf/10.1145/166117.166151?casa_token=A78GoJh369IAAAAA%3AfWlns5xYUUWiQC09j8IIMBVDZq4jYTRt_uJX2bQhc4V5ScK31VPFtDIXplbXd6dQro2TZ89ipnXo)参照．LICを用いて描いた流線の例として，[Phillips et al., 2015](https://iopscience.iop.org/article/10.1088/1748-3190/10/5/056020/meta)や[Van Der Kindere et al., 2019](https://link.springer.com/article/10.1007/s00348-019-2678-5)などがあります．
[SciPy Cookbook](https://scipy-cookbook.readthedocs.io/items/LineIntegralConvolution.html)にサンプルコードがあるのでありがたく使わせてもらいます．```setup.py```をビルドすると使えるようになります．
リンク先の```lic_demo.py```の前半を書き換えます．背景に流速を表示させたいので最後の行でLICと流速をかけています．

```python
Vector = np.stack([u, v], axis=2)
Vector = np.array(Vector, dtype=np.float32)
texture = np.random.rand(n, n).astype(np.float32)

plt.bone()
frame = 0
dpi = 1000
video = False
if video:
    kernellen = 31
    for t in np.linspace(0,1,16*5):
        kernel = np.sin(np.arange(kernellen)*np.pi/kernellen)*(1+np.sin(2*np.pi*5*(np.arange(kernellen)/float(kernellen)+t)))
        kernel = kernel.astype(np.float32)
        image = lic_internal.line_integral_convolution(Vector, texture, kernel)
        frame += 1
else:
    kernellen = 31
    kernel = np.sin(np.arange(kernellen)*np.pi/kernellen)
    kernel = kernel.astype(np.float32)
    image = lic_internal.line_integral_convolution(Vector, texture, kernel)

speed_LIC = image[1:-1, 1:-1]*speed[1:-1, 1:-1]
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/119868/d606a4e7-d6c3-c074-14a3-2747cbe70322.png" width=90%>

とても綺麗ですね．

#ちなみにMathematicaだともっと簡単
Mathematicaには組み込みで[```LineIntegralConvolutionPlot```](https://reference.wolfram.com/language/ref/LineIntegralConvolutionPlot.html)がある！

```mathematica
LineIntegralConvolutionPlot[{{Sin[x] Cos[y], -Sin[y] Cos[x]}, {"noise", 1000, 1000}}, {x, 0, 2 Pi}, {y, 0, 2 Pi},  
ColorFunction -> "BlueGreenYellow", RasterSize -> 300]
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/119868/8e3e9cda-a7bd-f2fa-d250-af7b7f541ca2.png" width=80%>


めちゃくちゃ間単にできた...