Python + gnuplot でジュリア集合を描く
かなりまえのエントリーでマンデルブロ集合を描くっていのをやりましたけど、
今度はジュリア集合を描いてみよっていうことで、
前につくったプログラムをちょっと変えてジュリア集合を描いてみました。
また例によって、
- Pyrhonでプロットデータ出力
- gnuplotでグラフ描画&画像ファイル出力
っていうプロセスでやってみました。
ジュリア集合とは?
という漸化式で定義されるがで発散しない複素数zの初期値の集合のこと。
(cは複素数の定数)
つまり、今度はcではなくを変えていってグラフをプロットする、と。
cは適当に複素数を決めておけば、いいんだけども、
このcの値によっていろいろと変わったグラフがプロットされるのがこいつ。
ソースコード
この前のマンデルブロ集合をプロットするときにつかったソースコードを少し変えただけです。
ちょっとした説明は前回をみればたぶん大丈夫。
【julia.py】
# -*- coding: utf-8 -*- import sys def juliaSet(z, c, K, LOOPMAX): #発散まで回したループ数を返す n = 0 while (n<LOOPMAX and abs(z)<K): z = z**2 + c n += 1 return n def getArg(): #コマンドライン引数用関数 argv = sys.argv argc = len(argv) if (argc != 10): #引数がない場合は終了 print ('Usage: python %s filename c_real c_image' 'minReal maxReal minImage maxImage step loop') %argv[0] quit() return argv def main(): argv = getArg() writer = open(argv[1], 'w') loop = int(argv[9]); K = 100.0 step = float(argv[8]) x_min = float(argv[4]); x_max = float(argv[5]) y_min = float(argv[6]); y_max = float(argv[7]) c = float(argv[2]) + (float(argv[3]) * 1j) x = x_min while (x < x_max): y = y_min while (y < y_max): z0 = x + (y * 1j) n = juliaSet(z0, c, K, loop) writer.write(str(z0.real)+'\t'+str(z0.imag)+'\t'+str(n)+'\n') y += step writer.write('\n') x += step print "finish data making." if __name__ == '__main__': main()
【juliaplot.sh】
#!/bin/sh #引数チェック if [ $# -ne 9 ]; then echo "Usage: $0 filname c_real c_image minReal maxReal minImage maxImage step loop" 1>&2 exit 1 fi #pythonによるジュリア集合データ出力 python julia.py $1.data $2 $3 $4 $5 $6 $7 $8 $9 # #gnuplotによるプロット echo "set term png; set output \"$1\"; set grid; set pm3d map; set size square; splot \"$1.data\"" | gnuplot #openはmacで使える open $1 exit 0
実行結果
まずは、としたときのプロット。
プロット範囲は実部・虚部共に1.6~-1.6で、0.06毎に、ループ上限7000回で動かしてみます。
$ ./juliaplot.sh julia01 -0.74759 0.083587 -1.6 1.6 -1.6 1.6 0.006 7000
個人的にはジュリア集合っていうとこんな感じのなんだけれども、
パラメータを変えていくとけっこう面白い形が見えてきます。
でプロットしてみる。
$ ./juliaplot.sh julia02 0.250166 -0.000004 -1.5 1.5 -1.5 1.5 0.005 3000
さっきとは全く違った形に。
最後に、でプロット。
$ ./juliaplot.sh julia03 -0.125475 0.7494 -1.4 1.4 -1.4 1.4 0.006 700
まとめ的なものは前回と一緒なので省略。
いい感じの暇つぶしにはなりそう。
わりと時間がかかる処理なので、ベンチマーク代わりにつかったり、CUDAで実装して性能を比較してみたりするにはいいかもしれないですね。