コッホ曲線（Mathematicaショートフラクタルシリーズ）

## 今回のテーマ

Mathematicaショートフラクタルシリーズ、第3回は、フラクタルの説明で最も使われているであろう、[コッホ曲線](https://ja.wikipedia.org/wiki/%E3%82%B3%E3%83%83%E3%83%9B%E6%9B%B2%E7%B7%9A)です。

![koch.png](https://qiita-image-store.s3.amazonaws.com/0/6117/8e513449-746f-37d1-3877-3ac6fb055475.png)

## ほぼ完成形

```mathematica:koch_original.nb
generator = {0, Pi/3, -Pi/3, 0};
f[list_] := list + # & /@ generator
inclinationList = Flatten[Nest[f, generator, 5]];
cPoints = Join[{0}, Accumulate[E^(I inclinationList)]];
Graphics[Line[Transpose[{Re[cPoints], Im[cPoints]}]]]
```

## 偏角関数

今回のコッホ曲線描画アルゴリズムは、
僕が学部生のときに発表した研究内容に基づきます。

自主課題研究「[偏角関数を用いた曲線の研究](http://butchi.jp/documents/b3-research.pdf)」(2006) ([PPT発表スライド](http://butchi.jp/documents/b3-research.ppt))
<img width="850" alt="koch_inclination.png" src="https://qiita-image-store.s3.amazonaws.com/0/6117/7f119b3b-69e1-f91a-e2c1-ee1ebe8e6265.png">

この図は、コッホ曲線の第0段階から第2段階までを左から走査して横軸を長さ、縦軸を角度とした偏角関数を示したものです。

偏角関数の
第0段階は $0$
第1段階は $0, \frac\pi3, -\frac\pi3, 0$
第2段階は $0, \frac\pi3, -\frac\pi3, 0, \frac\pi3, \frac{2\pi}3, 0, \frac\pi3, \dots$

という列と捉えられます、$\frac\pi3$ で除算すれば
第1段階 $0$
第2段階 $0, 1, -1, 0$
第3段階 $0, 1, -1, 0, 1, 2, 0, 1, \dots$

となり見やすくなります。
（Mathematicaで実装する際は最初から $\frac\pi3$ が掛かっててもいいし、整数列を作って単位を掛けてもいいです。Cとか一般的な言語だと後者が賢明でしょうが、今回は前者です。）

フラクタルのパーツと偏角関数の関係は以下の図でより伝わるかと思います。

![koch_angular-function_image.png](https://qiita-image-store.s3.amazonaws.com/0/6117/74e115d7-e404-6321-9873-636fb28e1da2.png)

## 生成過程

```
generator = {0, Pi/3, -Pi/3, 0};
```

これが第1段階です。

```
f[list_] := list + # & /@ generator
```

これは偏角リストにジェネレーターの各定数を足し、

```
{
  {0, Pi/3, -Pi/3, 0} + 0,
  {0, Pi/3, -Pi/3, 0} + Pi/3,
  {0, Pi/3, -Pi/3, 0} + -Pi/3,
  {0, Pi/3, -Pi/3, 0} + 0
}
```

と計算します。（リスト同士の加算は次元を合わせる必要がありますが、リストと定数の加算では各要素に定数が足されます。）

```
inclinationList = Flatten[Nest[f, generator, 5]];
```

で第6段階まで進みます。
最終的に一気に平坦化するので、その時点で一連の偏角リストが得られます。

そして次のように偏角から複素平面上の座標を構築します。
オイラーの公式より、 $e^{a i}$ が1辺先となります。
（先頭に原点の``0``を忘れずに）

```
cPoints = Join[{0}, Accumulate[E^(I inclinationList)]];
```

[複素数の魔法](http://qiita.com/butchi_y/items/a4a7e992870d75e5ee7a)で曲線をプロットします。

```
Graphics[Line[Transpose[{Re[cPoints], Im[cPoints]}]]]
```

## ショートコード化

- ワンライナー化（2重の無名関数をうまく作れなかったのでやむを得ず``Function``を用意）
- ``Accumulate``と``Join``を併用するより``Plus``を``FoldList``した方が端的

結果、115文字のプログラムができました。

```mathematica:koch_short.nb
Graphics@Line@Transpose@{Re[#],Im[#]}&@FoldList[Plus,0,E^(Flatten@Nest[Function[l,l+#&/@{0,Pi/3,-Pi/3,0}],{0},5]I)]
```

<blockquote class="twitter-tweet" lang="ja"><p lang="en" dir="ltr"><a href="https://twitter.com/butchi_y">@butchi_y</a> (More info: <a href="https://t.co/HBUuaZOB3w">https://t.co/HBUuaZOB3w</a>) <a href="https://twitter.com/hashtag/wolframlang?src=hash">#wolframlang</a> <a href="https://t.co/qj9AvcqEFN">pic.twitter.com/qj9AvcqEFN</a></p>&mdash; Tweet-a-Program (@wolframtap) <a href="https://twitter.com/wolframtap/status/676890608229867520">2015, 12月 15</a></blockquote>
<script async src="//platform.twitter.com/widgets.js" charset="utf-8"></script>


## デモ

ちなみに、これをもっと一般化したデモをWolframサイトに掲載しています。
これまた有名なフラクタル図形を多数再現できるほか、manualモードで角度も自由に変えられます。
僕のフラクタル研究の集大成のひとつですね。

[Angular Function Fractals - Wolfram Demonstrations Project](http://demonstrations.wolfram.com/AngularFunctionFractals/)
[![Terdragon](http://demonstrations.wolfram.com/AngularFunctionFractals/HTMLImages/index.en/popup_1.jpg)](http://demonstrations.wolfram.com/AngularFunctionFractals/)


## 予告

明日は[ドラゴン曲線](https://ja.wikipedia.org/wiki/%E3%83%89%E3%83%A9%E3%82%B4%E3%83%B3%E6%9B%B2%E7%B7%9A)です。