gnuplot は, プログラムちっくなお絵かきするのに優れています.
あまりご存じないかもしれませんが, gnuplot では条件分岐, ループ, 関数等のプログラミングをするには十分な機能を備えています.
再帰呼び出しをすることもできます.
再帰呼び出しと言えば, フラクタル図形ですね.
gnuplot でフラクタル図形を描くことはとても簡単です.
シェルピンスキーのギャスケット
とても簡単に描けます.
reset number = 8 sizex = 900 sizey = 900 wid = 0.95 pad = (1 - wid) / 2 * 100; padbottom = (1 - sizex * wid / 2 * sqrt(3) / sizey) / 2 set terminal jpeg size sizex,sizey set output "sierpinski.jpeg" set lmargin pad set rmargin pad set tmargin 0 set bmargin 0 set origin 0,padbottom set border 0 set notitle set nokey unset xtics unset ytics set xrange [0:1] set yrange [0:1.0/sizex*sizey/wid] set style arrow 1 nohead linewidth 1 sierpinski(n, x, y, dx) \ = n >= number ? \ sprintf("set object polygon from %f,%f to %f,%f to %f,%f to %f,%f fillstyle solid;", x, y, x + dx, y, x + dx / 2, y + dx / 2 * sqrt(3), x, y) : \ sierpinski(n + 1, x, y, dx / 2). \ sierpinski(n + 1, x + dx / 2, y, dx / 2). \ sierpinski(n + 1, x + dx / 4, y + dx / 4 * sqrt(3), dx / 2); eval(sierpinski(0, 0.0, 0.0, 1.0)) plot -100
再帰的に関数を呼び, sprintf で生成した大量の arrow を文字列結合演算子「.」で結合しまくり, 最後に eval すれば描画できます.
条件分岐には, 三項演算子を使います.
シェルピンスキーのカーペット
こちらは, 直ぐに再帰爆発するから, あまりnumberを増やせない.
reset number = 5 sizex = 800 sizey = 800 wid = 0.95 pad = (1 - wid) / 2 * 100; padbottom = (1 - sizex * wid / sizey) / 2 set terminal jpeg size sizex,sizey set output "carpet.jpeg" set lmargin pad set rmargin pad set tmargin 0 set bmargin 0 set origin 0,padbottom set border 0 set notitle set nokey unset xtics unset ytics set xrange [0:1] set yrange [0:1.0/sizex*sizey/wid] set style arrow 1 nohead linewidth 1 carpet(n, x, y, d) \ = n >= number ? \ sprintf("set object rect from %f,%f to %f,%f fc rgb '#000000' fs solid;", x, y, x + d, y + d) : \ carpet(n + 1, x, y, d / 3) . \ carpet(n + 1, x + d / 3, y, d / 3) . \ carpet(n + 1, x + 2 * d / 3, y, d / 3) . \ carpet(n + 1, x, y + d / 3, d / 3) . \ carpet(n + 1, x + 2 * d / 3, y + d / 3, d / 3) . \ carpet(n + 1, x, y + 2 * d / 3, d / 3) . \ carpet(n + 1, x + d / 3, y + 2 * d / 3, d / 3) . \ carpet(n + 1, x + 2 * d / 3, y + 2 * d / 3, d / 3); eval(carpet(0, 0.0, 0.0, 1.0)) plot -100
コッホ曲線
回転行列が計算できたら簡単に描けます.
reset number = 6 sizex = 1200 sizey = 400 wid = 0.9 pad = (1 - wid) / 2 * 100; padbottom = (1 - sizex * wid / 6 * sqrt(3) / sizey) / 2 set terminal jpeg size sizex,sizey set output "koch.jpeg" set lmargin pad set rmargin pad set tmargin 0 set bmargin 0 set origin 0,padbottom set border 0 set notitle set nokey unset xtics unset ytics set xrange [0:1] set yrange [0:1.0/sizex*sizey/wid] set style arrow 1 nohead linewidth 1 rotx(x,y,t) = cos(t) * x - sin(t) * y roty(x,y,t) = sin(t) * x + cos(t) * y koch(n, x, y, dx, dy) \ = n >= number ? \ sprintf("set arrow from %f,%f to %f,%f as 1;", x, y, x+dx, y+dy) : \ koch(n + 1 \ , x \ , y \ , dx / 3 \ , dy / 3). \ koch(n + 1 \ , x + dx / 3 \ , y + dy / 3 \ , rotx(dx / 3, dy / 3, pi / 3) \ , roty(dx / 3, dy / 3, pi / 3)). \ koch(n + 1 \ , x + dx / 3 + rotx(dx / 3, dy / 3, pi / 3) \ , y + dy / 3 + roty(dx / 3, dy / 3, pi / 3) \ , rotx(dx / 3, dy / 3, - pi / 3) \ , roty(dx / 3, dy / 3, - pi / 3)). \ koch(n + 1 \ , x + 2 * dx / 3 \ , y + 2 * dy / 3 \ , dx / 3 \ , dy / 3); eval(koch(0, 0.0, 0.0, 1.0, 0.0)) plot -100
高木曲線
難しそうに思える高木曲線も, 意外と簡単なプログラムで描けます.
これも, number を変化させて遊んでみて下さい.
reset number = 9 sizex = 400 sizey = 400 wid = 0.95 pad = (1 - wid) / 2 * 100; padbottom = (1 - sizex * wid / 2 * sqrt(3) / sizey) / 2 set terminal jpeg size sizex,sizey set output "takagi.jpeg" set lmargin pad set rmargin pad set tmargin 0 set bmargin 0 set origin 0,padbottom set border 0 set notitle set nokey unset xtics unset ytics set xrange [0:1] set yrange [0:1.0/sizex*sizey/wid] set style arrow 1 nohead linewidth 1 takagi(n, x, y, dx, dy) \ = n >= number ? \ sprintf("set arrow from %f,%f to %f,%f as 1;", x, y, x + dx, y + dy) : \ takagi(n + 1, x, y, dx / 2, dy / 2 + dx / 2). \ takagi(n + 1, x + dx / 2, y + dy / 2 + dx / 2, dx / 2, dy / 2 - dx / 2); eval(takagi(0, 0.0, 0.0, 1.0, 0.0)) plot -100
ドラゴン曲線
こちらもシンプルなシステム.
reset number = 12 sizex = 600 sizey = 600 wid = 0.95 pad = (1 - wid) / 2 * 100; padbottom = (1 - sizex * wid / sizey) / 2 set terminal jpeg size sizex,sizey set output "dragon.jpeg" set lmargin pad set rmargin pad set tmargin 0 set bmargin 0 set origin 0,padbottom set border 0 set notitle set nokey unset xtics unset ytics set xrange [0:1] set yrange [0:1.0/sizex*sizey/wid] set style arrow 1 nohead linewidth 1 hilbert(n, x, y, dx, dy) \ = n >= number ? \ sprintf("set arrow from %f,%f to %f,%f as 1;", x, y, x + dx, y + dy) : \ hilbert(n + 1, x, y, (dx - dy) / 2, (dy + dx) / 2) . \ hilbert(n + 1, x + dx, y + dy, - (dx + dy) / 2, (dx - dy) / 2); eval(hilbert(0, 0.2, 0.4, 0.7, 0.0)) plot -100
ヒルベルト曲線
これが, 一番アタマを捻ったし, 時間がかかった.
隣と接続するのが結構難しい.
でも, きちんと動いたプログラムを見なおしたら大したことなかった.
手で図をきちんと描けば分かる.
reset number = 4 sizex = 600 sizey = 600 wid = 0.95 pad = (1 - wid) / 2 * 100; padbottom = (1 - sizex * wid / sizey) / 2 set terminal jpeg size sizex,sizey set output "hilbert.jpeg" set lmargin pad set rmargin pad set tmargin 0 set bmargin 0 set origin 0,padbottom set border 0 set notitle set nokey unset xtics unset ytics set xrange [0:1] set yrange [0:1.0/sizex*sizey/wid] set style arrow 1 nohead linewidth 1 hilbert(n, x, y, dx, dy, dx2, dy2, dx3, dy3) \ = n >= number ? \ sprintf("set arrow from %f,%f to %f,%f as 1;", x + (dx - dy) / 4, y + (dy + dx) / 4, x + (3 * dx - dy) / 4, y + (3 * dy + dx) / 4) . \ sprintf("set arrow from %f,%f to %f,%f as 1;", x + (3 * dx - dy) / 4, y + (3 * dy + dx) / 4, x + (3 * dx - 3 * dy) / 4, y + (3 * dy + 3 * dx) / 4) . \ sprintf("set arrow from %f,%f to %f,%f as 1;", x + (dx - 3 * dy) / 4, y + (dy + 3 * dx) / 4, x + (3 * dx - 3 * dy) / 4, y + (3 * dy + 3 * dx) / 4) . \ sprintf("set arrow from %f,%f to %f,%f as 1;", x + (dx - 3 * dy) / 4, y + (dy + 3 * dx) / 4, x + (dx - 3 * dy) / 4 + dx2 / 2, y + (dy + 3 * dx) / 4 + dy2 / 2) . \ sprintf("set arrow from %f,%f to %f,%f as 1;", x + (dx - dy) / 4, y + (dy + dx) / 4, x + (dx - dy) / 4 + dx3 / 2, y + (dy + dx) / 4 + dy3 / 2) : \ hilbert(n + 1, x + dx / 2, y + dy / 2, - dy / 2, dx / 2, dx3 / 2, dy3 / 2, dx / 2, dy / 2). \ hilbert(n + 1, x + dx / 2, y + dy / 2, dx / 2, dy / 2, - dy / 2, dx / 2, 0.0, 0.0). \ hilbert(n + 1, x + dx / 2 - dy / 2, y + dy / 2 + dx / 2, dx / 2, dy / 2, -dx / 2, -dy / 2, 0.0, 0.0). \ hilbert(n + 1, x - dy, y + dx, dy / 2, - dx / 2, 0.0, 0.0, dx2 / 2, dy2 / 2); eval(hilbert(0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0)) plot -100