 |
DLA(Diffusion-limited aggregation）とは日本語では拡散律速凝集といい、ブラウン運動する粒子が、核となる粒子に結合することで粒子の集合体（クラスタ）が成長する過程のことです。その集合体を凝集体といいます。雪の結晶や雷の成長などに見られ、美しいフラクタル構造となっていることが知られています。ブラウン運動は６章で詳しく勉強したので、その応用としてDLAを描いてみましょう。
なお、この章のプログラムは、16ビットのTurboCでは、配列のメモリが不足するため、RasoberryPIのUNIXで実行します。
はじめの例は元になる核は地表とし、空中の１点を発した雪が降りそそぎ、地表に到達すると、下から結合され、徐々に積み上がり、樹氷のようになる姿を描きます。雪は直線的に下に落ちるのではなく、右に左にゆらゆらと落下します。その様子はランダムウォークそのものです。
プログラムのポイント
シミュレーションプログラムをs1110.cに示します。
行15のように、画面に相当する状態配列を用意します。値は1あるいは0で、その位置にすでに雪があるかないかを示すこととします。初期状態は行26のように、地表のみ、1としておきます。
雪は、行31のように、画面中央トップから降りそそぐとしました。行39から40が雪はゆらゆら降るさまを示したもので、ランダムウォークです。つぎのタイミングで、粒子の位置を、x , yとし、右に行く割合が40%、左に行く割合が40%、残り20%が真下に落ちるとしています。雪が下に位置を変えるたびに、行41でその位置の状態がすでに雪で占められているかどうか調べます。雪がすでにあれば、行46のように一つ手前 x0, y0 にもどり、その位置の状態を、1に変えます。これが凝集体に捕捉された処理です。
適当な高さに成長したら、処理を終えます。
アウトプット
図１１．１．１に示します。
| s1110.c 樹氷 |
/* S1110.C
* DLA 樹氷
* (C) H.Ishikawa 2018/11/10
*/
#include "compati.h"
#include "window.h"
#define X_MAX 640
#define Y_MAX 480
int main(void){
int x0,y0; //雪の現在地
int x,y; //次の位置
int s[X_MAX + 10][Y_MAX + 10] = {}; //状態配列 1：雪がある 0：雪がない
float random;//乱数
/*グラフィクの準備*/
opengraph();
SetWindow(0, 0, 0, X_MAX, Y_MAX, 0,HIGHT, WIDTH,0);
srand(17);
x = 0.0;
y = 0.0;
x0 = 0.0;
y0 = 0.0;
for (x = 0; x <= X_MAX ; x ++) {
s[x][0] = 1; //初期状態 地上を状態1に
}
while(y0 < Y_MAX - 20){//高さがY_MAX - 20まで成長したらとめる
//雪の起点
x = X_MAX / 2;
y = Y_MAX;
while (y > 0){
random = RAND();
if(random > 0.6){//雪は右へ
if (x < X_MAX) { x = x + 1; }else{ x = X_MAX;}
}else if (random < 0.4){雪は左へ
if (x > 0) { x= x - 1;}else{ x = 0;}
}
y = y - 1;//雪は下に
if (s[x][y] == 0){//この位置に雪がないとき
PutPixel(0, (double)x, (double)y ,1);//雪の軌跡を青で描く
x0 = x;
y0 = y;
}else{//この位置にすでに雪があれば一つ前の位置を状態1に
s[x0][y0] = 1;//捕捉
PutPixel(0, (double)x0, (double)y0 ,15);//白をプロット
break;//次の雪を降らせる
}
}
}
getch();
closegraph();
return 0;
}
|
| 図１１．１．１ s1110.c 樹氷 |
 |
|
次の例はDLAとして有名な凝集体です。金属イオンが溶けた溶液に電極を入れると、金属が析出し、徐々に植物のような形に育っていきます。これを金属樹、金属葉とよばれます。これもフラクタル構造になることが知られています。
プログラムのポイント
下の図１１．２．１に従って凝集体のシミュレーションを説明します。まず核を一つ画面中央におきます。粒子は中央から等距離の円周上から、ランダムに発生させます。そこから2次元のランダムウォークをさせます。その位置x ,yの状態を調べ空いていれば（0であれば）ランダムウォークを続けます。もしすでに粒子で占められていれば（1であれば）、ひとつ戻って、隣接する位置x0 ,y0の状態を1に遷移させます。この粒子はここで捕捉されましたので、次の粒子を発生させ処理を繰り返します。
プログラムS1120.Cは、凝集体のシミュレーションです。x0 ,y0の状態を1に変化させるときに、白をプロットします。中央の核に粒子が捕捉され、凝集がはじまります。
 |
| 図１１．２．１ 凝集体のシミュレーションの方法 文献（１６） |
| s1120.c 凝集体のシミュレーション |
/* s1120.c
* 凝集体のシミュレーション
* (C) H.Ishikawa 2018/11/20
*/
#include <math.h> // sin() cos() M_PI
#include "compati.h"
#include "window.h"
#define X_MAX 480
#define Y_MAX 480
#include <time.h>
int main(void){
int x0,y0; //粒子現在位置
int x,y; //粒子次の位置
int s[X_MAX + 10][Y_MAX + 10] = {};//状態 1は粒子が固着
long t = 0;
int color;
double th;
double r;
/*グラフィクの準備*/
opengraph();
SetWindow(0, 0, 0, X_MAX, Y_MAX, (WIDTH - HIGHT) / 2, HIGHT, (WIDTH + HIGHT) / 2,0);
/*メイン*/
srand(107);
x = 0;
y = 0;
x0 = 0;
y0 = 0;
s[X_MAX / 2][Y_MAX / 2] = 1;//中央に核をおく
while(t <= 15000){
/*粒子発生*/
th = 2.0 * M_PI *RAND();
x = (int)(X_MAX /2.0 * cos(th) + X_MAX /2.0) ;
y = (int)(Y_MAX /2.0 * sin(th) + Y_MAX /2.0) ;
while (1){
/*ランダムウォーク*/
r = RAND();
if(r>0.67){//右へ
if (x < X_MAX) { x = x + 1; }else{ x = X_MAX;}
}else if (r < 0.33){//左へ
if (x > 0) { x= x - 1;}else{ x = 0;}
}//このほかはとどまる
r = RAND();
if(r>0.67){//上へ
if (y < Y_MAX) { y = y + 1; }else{ y = Y_MAX;}
}else if (r < 0.33){//下へ
if (y > 0) { y= y - 1;}else{ y = 0;}
}//このほかはとどまる
/*判定*/
if (s[x][y] != 1){
PutPixel(0, (double)x, (double)y,1);
x0 = x;
y0 = y;
}else{/*捕捉*/
s[x0][y0] = 1;
t = t + 1;
PutPixel(0, (double)x0, (double)y0 ,15);
break;/*次の粒子へ*/
}
}
}
getch();
closegraph();
return 0;
}
|
アウトプット
s1120.c 実行します。結果を図１１．２．２に示します。処理時間約10分
| 図１１．２．２ s1120.c 凝集体のシミュレーション |
 |
|
忍石（しのぶいし）とは、石灰岩・頁岩(けつがん)などの割れ目に二酸化マンガンなどの樹枝状結晶が付着して、シダ類の化石のように見えるものです。江戸時代の石のコレクター、木内石亭の著書「雲根志」には、忍草石として登場しており珍重されました。
この形で成長する結晶は多く、冬の窓に付く霜の雪片もこの一種です。
| 図１１．３．１ しのぶ石の例 （Wikipedia） |
 |
|
プログラムのポイント
しのぶ石のシミュレーションプログラムはs1130.c のとおりです。
初期化は行26のように地表のみ、1としておきます。粒子は、行30のように画面全体の任意の位置からランダムに発生するとします。粒子のランダムウォークは2次元で、1/3の割合で右、左、とどまるとします。上下方向も同様です。凝集体に捕捉される処理はｓ1120.ｃと同様です。
| s1130.c しのぶ石 |
/* S1130.C
* しのぶ石
* (C) H.Ishikawa 2018/11/10
*/
#include "compati.h"
#include "window.h"
#define X_MAX 640
#define Y_MAX 480
int main(void){
long x0,y0; //現在地
long x,y; //次の位置
long s[X_MAX + 1][Y_MAX + 1] = {}; //状態配列 1：粒子がある 0：粒子がない
float r;//乱数
/*グラフィクの準備*/
opengraph();
SetWindow(0, 0, 0, X_MAX, Y_MAX, 10,HIGHT-10, WIDTH-10,10);
srand(1003);
x = 0.0;
y = 0.0;
x0 = 0.0;
y0 = 0.0;
//初期化
for (x = 0; x <= X_MAX ; x ++) {
s[x][0] = 1; //地上のところどころを状態1に
}
while(y0 < Y_MAX - 20){//高さがY_MAX - 20まで
//粒子の起点は画面全体
x = (int)(RAND() * X_MAX);
y = (int)(RAND() * Y_MAX);
while (y > 0){
//２次元のランダムウォーク
r = RAND();
if(r>0.67){//右へ
if (x < X_MAX) { x = x + 1; }else{ x = X_MAX;}
}else if (r < 0.33){//左へ
if (x > 0) { x= x - 1;}else{ x = 0;}
}//このほかはとどまる
r = RAND();
if(r>0.67){//上へ
if (y < Y_MAX) { y = y + 1; }else{ y = Y_MAX;}
}else if (r < 0.33){//下へ
if (y > 0) { y= y - 1;}else{ y = 0;}
}//このほかはとどまる
if (s[x][y] == 0){//この位置に粒子がないとき
PutPixel(0, (double)x, (double)y ,1);//粒子の軌跡
x0 = x;
y0 = y;
}else{//この位置にすでに粒子があれば一つ前の位置を状態1に
s[x0][y0] = 1;
PutPixel(0, (double)x0, (double)y0 ,15);//白をプロット
break;//次の粒子を発生させる
}
}
}
getch();
closegraph();
return 0;
}
|
アウトプット
s1130.c を実行します。結果を図１１．３．２に示します。これを描くのに、約30分かかります。
| 図１１．３．２ s1130.c しのぶ石シミュレーション |
 |
|
１．核をどのように初期条件としてあたえるかによって、さまざまなDLAを描くことができます。核を円周上に置き、粒子のスタートは、中央としたとき、どのようなDLAが描けるか、トライしなさい。
（解答）
| a1110.c 円周上に核 |
/* a1110.c
* DLA 円周上に核
* (C) H.Ishikawa 2018/11/20
*/
#include <math.h> // sin() cos()
#include "compati.h"
#include "window.h"
#define X_MAX 480
#define Y_MAX 480
int main(void)
{
int x0,y0; //粒子現在位置
int x,y; //粒子次の位置
int s[X_MAX + 10][Y_MAX + 10] = {};//状態 1は粒子が固着
long t = 0;
int color;
double th;
double r;
/*グラフィクの準備*/
opengraph();
SetWindow(0, 0, 0, X_MAX, Y_MAX, (WIDTH - HIGHT) / 2, HIGHT, (WIDTH + HIGHT) / 2,0);
/*メイン*/
srand(13);
x = 0;
y = 0;
x0 = 0;
y0 = 0;
//初期化 円周上に核
for (th = 0 ; th < 2.0 * M_PI; th = th + 0.0001){
x = (int)(X_MAX /2.0 * cos(th) + X_MAX /2.0) ;
y = (int)(Y_MAX /2.0 * sin(th) + Y_MAX /2.0) ;
s[x][y] = 1;
}
while(t < 25000){
/*粒子発生は画面中央*/
x = (int)(X_MAX /2.0) ;
y = (int)(Y_MAX /2.0) ;
while (1){
/*ランダムウォーク*/
r = RAND();
if(r>0.67){//右へ
if (x < X_MAX) { x = x + 1; }else{ x = X_MAX;}
}else if (r < 0.33){//左へ
if (x > 0) { x= x - 1;}else{ x = 0;}
}//このほかはとどまる
r = RAND();
if(r>0.67){//上へ
if (y < Y_MAX) { y = y + 1; }else{ y = Y_MAX;}
}else if (r < 0.33){//下へ
if (y > 0) { y= y - 1;}else{ y = 0;}
}//このほかはとどまる
/*判定*/
if (s[x][y] != 1){
PutPixel(0, (double)x, (double)y,1);
x0 = x;
y0 = y;
}else{/*捕捉*/
s[x0][y0] = 1;
t = t + 1;
PutPixel(0, (double)x0, (double)y0 ,15);
break;/*次の粒子へ*/
}
}
}
getch();
closegraph();
return 0;
}
|
| 図１１．a．１ a1110.c の結果 |
 |
|
２．核を正方形の４辺 に置き、粒子のスタートは、画面全体としたとき、どのようなDLAが描けるか、トライしなさい。
（解答）
| a1120.c 正方形の４辺に核を置く |
/* a1120.c
* DLA ４辺に核
* (C) H.Ishikawa 2018/11/27
*/
#include <math.h> // sin() cos()
#include "compati.h"
#include "window.h"
#define X_MAX 480
#define Y_MAX 480
int main(void)
{
int x0,y0; //粒子現在位置
int x,y; //粒子次の位置
int s[X_MAX + 10][Y_MAX + 10] = {};//状態 1は粒子が固着
long t = 0;
int color;
double th;
double r;
/*グラフィクの準備*/
opengraph();
SetWindow(0, 0, 0, X_MAX, Y_MAX, (WIDTH - HIGHT) / 2 +10, HIGHT-10, (WIDTH + HIGHT) / 2 -10 ,10);
/*メイン*/
srand(103);
x = 0;
y = 0;
x0 = 0;
y0 = 0;
//初期化 正方形の４辺に核を置く
for (x = 0; x <= X_MAX ; x ++) {
s[x][0] = 1;
}
for (x = 0; x <= X_MAX ; x ++) {
s[x][Y_MAX] = 1;
}
for (y = 0; y <= Y_MAX ; y ++) {
s[0][y] = 1;
}
for (y = 0; y <= Y_MAX ; y ++) {
s[X_MAX][y] = 1;
}
while(t <100000){
//粒子の起点は画面全体
x = (int)(RAND() * X_MAX);
y = (int)(RAND() * Y_MAX);
while (1){
/*ランダムウォーク*/
r = RAND();
if(r>0.67){//右へ
if (x < X_MAX) { x = x + 1; }else{ x = X_MAX;}
}else if (r < 0.33){//左へ
if (x > 0) { x= x - 1;}else{ x = 0;}
}//このほかはとどまる
r = RAND();
if(r>0.67){//上へ
if (y < Y_MAX) { y = y + 1; }else{ y = Y_MAX;}
}else if (r < 0.33){//下へ
if (y > 0) { y= y - 1;}else{ y = 0;}
}//このほかはとどまる
/*判定*/
if (s[x][y] != 1){
PutPixel(0, (double)x, (double)y,0);
x0 = x;
y0 = y;
}else{/*捕捉*/
s[x0][y0] = 1;
t = t + 1;
//経過時間にしたがって色を変える
PutPixel(0, (double)x0, (double)y0 ,(int)(t/6000 + 1));
break;/*次の粒子へ*/
}
}
}
getch();
closegraph();
return 0;
}
|
| 図１１．a．２ a1120.c の結果 |
 |
|
|
 |
