Buddhabrot 描画プログラムの考え方
計算式
Buddhabrot 画像描画には，以下に示す Mandelbrot 集合を求める時の漸化式を用いる．
\begin{equation}
\begin{cases}
x_{n+1}=x_n{}^2-y_n{}^2+a \\
y_{n+1}=2 x_n y_n+b
\end{cases}
\end{equation}
計算手順
- 計算・描画領域を定める．基本は，$-2 \leqq x \leqq 2$ ，$-2 \leqq y \leqq 2$ とする．
- 計算領域を有限個の微小領域に分割する．
- 乱数を用いて座標 (a,b) を定める．
- Mandelbrot集合を求める上述漸化式の値 $(x_{n+1},y_{n+1})$ を計算する．
- 計算された座標 $(x_{n+1},y_{n+1})$ が，分割された有限個の微小領域のどこに属するかを記憶する．
- 発散が確定 $(4 < x_{n+1}^2+y_{n+1}^2)$ するまで上記操作を繰り返す．
- 発散しない場合の値は破棄する．
- これにより，発散する場合について，分割された微小領域に何個の値が属するかを示すデータが得られる．
- 微小領域に存在するデータの個数により配色を定め，微小領域を塗りつぶすことにより Buddhabrot 画像が描画できる．
プログラミングの考え方
- 画像を支配する要素は，計算・描画領域，領域分割数，乱数発生回数，発散確定までの計算繰り返し数である． ここでは乱数発生回数は 1,000,000 回，発散確定までの計算繰り返し数は 10,000 回に固定し，計算・描画範囲，領域分割数，数値データ出力ファイル名をコマンドライン引数で入力する．
- 数値出力データファイル名をコマンドラインで入力するようにしたのは，計算繰り返し数を変更したい場合，100万回単位でデータを重ね合わせられるようにするためである．
- 乱数により (a,b) を定めるため，データの重ね合わせを有効にするため，random_seed を用いてプログラム起動毎に seed を変更するようにする．
- 作図は GMT で行うが，GMT スクリプトは Fortran プログラムには含めず，別途作成している．これもデータ重ね合わせの結果利用を考慮してのことである．
- 計算された座標 $(x_{n+1},y_{n+1})$ が，分割された有限個の微小領域のどこに属するかを記憶する．
- 配色は，白黒の濃淡を用いることとし，発散確定までの繰り返し数が少ない場合を黒，多い場合を白としている．
- GMT が吐き出す eps 画像は時としてファイルサイズが大きくなり扱い辛いので，最終成果は ImageMagick を用いて jpg 画像に変換している．
プログラム
Programs
| Filename | Description |
|---|---|
| f90_buddhabrot_gmt.txt | Plot coordinate calculation program |
| f90_add_data.txt | Plot data superposition program |
プロット座標計算用スクリプト( Fortran プログラム実行)
スクリプトでは，コンパイル後，プロット座標とグレースケールでの配色を定める数値を計算します．
gfortran -o f90_buddhabrot_gmt f90_buddhabrot_gmt.f90 ./f90_buddhabrot_gmt 200 200 0 0 2 fout200_1.txt ./f90_buddhabrot_gmt 200 200 0 0 2 fout200_2.txt
上記スクリプトでは，同じ作図領域にたいし，領域分割数を200とした場合を2ケース計算しています． 乱数の種をプログラム起動毎に買えているので，若干事なる画像が得られます．
Fortranプログラム実行用コマンド
./f90_buddhabrot_gmt jw ih x0 y0 rr filename
| jw | 作図の水平方向分割数(整数値) |
| ih | 作図の鉛直方向分割数(jwと同じ整数値) |
| x0 | 作図の中心座標のx値(実数) |
| y0 | 作図の中心座標のy値(実数) |
| rr | 作図半径(実数) |
| filename | 作図用数値データ出力ファイル名(txt) |
データ重ね合わせ実行スクリプト( Fortranプログラム実行)
スクリプトでは，コンパイル後，GMT による作図データの重ね合わせを行います．
gfortran -o f90_add_data f90_add_data.f90 ./f90_add_data 200 200 fout200_1.txt fout200_2.txt fout200_2m.txt
上記スクリプトでは，領域分割数 200 x 200 のデータが格納されたファイル fout200_1.txt と fout200_2.txt を重ね合わせ fout200_2m.txt を作成します．
GMT実行スクリプト
下記スクリプトの概要は以下のとおり．
- echo を用いて，GMT 用 cpt ファイルを作成します．
- 領域分割数 200 x 200，乱数発生回数 100 万回のデータファイル fout200_1.txt を用いて GMT により作図を行います．
- 領域分割数 200 x 200，乱数発生回数 100 万回のデータファイル fout200_2.txt を用いてGMTにより作図を行います．
- 領域分割数 200 x 200，データファイルの香さね合わせを行い乱数発生回数 200 万回としたデータファイル fout200_2m.txt を用いて GMT により作図を行います．
# ******************************************************** # Make color pallet file # z1 R1 G1 B1 z2 R1 G2 B2 (z1/z2: range of value) #********************************************************* echo "0 0 0 0 500 125 125 125" > _col.cpt echo "500 125 125 125 1000 255 255 255" >> _col.cpt echo "F 255 255 255" >> _col.cpt echo "B 0 0 0" >> _col.cpt range=0/200/0/200 scale=2/2 inp_dat=fout200_1.txt fig=_gmt.eps psxy $inp_dat -R$range -JX$scale -SC0.01 -C_col.cpt -m -P > $fig convert -trim -density 300 -rotate 180 $fig fig_buddha_0200_1.jpg inp_dat=fout200_2.txt fig=_gmt.eps psxy $inp_dat -R$range -JX$scale -SC0.01 -C_col.cpt -m -P > $fig convert -trim -density 300 -rotate 180 $fig fig_buddha_0200_2.jpg inp_dat=fout200_2m.txt fig=_gmt.eps psxy $inp_dat -R$range -JX$scale -SC0.01 -C_col.cpt -m -P > $fig convert -trim -density 300 -rotate 180 $fig fig_buddha_0200_2m.jpg
なお配色を決める際参考となる，各ファイル内の微小要素内のデータ個数の最大値を検索する awk スクリプトを以下に示しておきます．
gawk 'BEGIN{cmax=0}$0!=">"{if(cmax<$3)cmax=$3}END{print cmax}' fout200_1.txt
gawk 'BEGIN{cmax=0}$0!=">"{if(cmax<$3)cmax=$3}END{print cmax}' fout200_2.txt
gawk 'BEGIN{cmax=0}$0!=">"{if(cmax<$3)cmax=$3}END{print cmax}' fout200_2m.txt
作図事例
描画領域分割 200 x 200
| fig_buddha_0200_1.jpg | fig_buddha_0200_2.jpg | fig_buddha_0200_2m.jpg |
| 乱数発生100万回 | 乱数発生100万回 | 乱数発生200万回 |
 |
 |
 |
描画領域分割 400 x 400
| fig_buddha_0400_1m.jpg | fig_buddha_0400_2m.jpg | fig_buddha_0400_3m.jpg |
| 乱数発生100万回 | 乱数発生200万回 | 乱数発生300万回 |
 |
 |
 |
| fig_buddha_0400_4m.jpg | fig_buddha_0400_5m.jpg | |
| 乱数発生400万回 | 乱数発生500万回 | |
 |
 |
