セルオートマトンでシミュレーションをしよう(8)

樹枝状析出

モデルの詳細は「樹枝状析出」のコンテンツを参照。ここでは、具体的なプログラムを考える。 粒子は一番上の段のランダムな位置に出現し、ランダムウォークで動いてゆく。成長の「基盤」として、一番下の段には初めから粒子を敷き詰めておけば、成長の規則は「すでに静止している粒子に触れたら、そこで停止」だけでよい。”触れる”範囲は周囲８セルとする（上下左右の4セルだけにしてもよいが、その場合は少々工夫が必要）。

以下はプログラム例である 

#include<stdio.h>
#include<eggx.h>
#include<stdlib.h>
#define NX 100
#define NY 50
#define EDGE 4
main()
{ 
  int s[NX][NY];
  int x, y;
  int i,j,k,l;
  int particle,step,rseed;
  int win;

/* preparation for graphics */
  win=gopen(NX*EDGE,NY*EDGE);
  window(win,0,0,NX-1,NY-1);
  gsetbgcolor(win,"white");
  newcolor(win,"blue");
  layer(win,0,1); 

/* set initial state */
  for(i=0; i<NX; ++i){
    for(j=1; j<NY; ++j){
     s[i][j] = 0;
  }}
  for(i=0; i<NX; ++i){ s[i][0] =1;}


/* input step and rseed */
  scanf("%d%d",&particle,&rseed);
  srandom(rseed);

/* display initial state */
  gclr(win);
  for(i=0; i<NX-1; ++i){
    for(j=0; j<NY-1; ++j){
      if(s[i][j] == 1)fillrect(win,i,j,1,1);
  }}
  copylayer(win,1,0);

 for(k=0; k<particle; ++k){

/* start one particle */
    x = random()%NX;
    y = NY-2;

   for(;;){

 /* start one step */
     x += (random()%3-1);
     x = (x+NX)%NX;
     y += (random()%3-1);
     if(y == NY-1) { break; }
     else if(
        s[(x+NX-1)%NX][y+1] == 1 || s[x][y+1] == 1 || s[(x+NX+1)%NX][y+1] ==1 ||
        s[(x+NX-1)%NX][y]==1     || s[(x+NX+1)%NX][y] ==1 ||
        s[(x+NX-1)%NX][y-1] == 1 || s[x][y-1] == 1 || s[(x+NX+1)%NX][y-1] ==1 )
       {s[x][y] = 1; break;} 
  
/* display one particle */
     gclr(win);
     for(i=0; i<NX-1; ++i){
      for(j=0; j<NY-1; ++j){
        if(s[i][j] == 1)fillrect(win,i,j,1,1);
     }}
     if(y != NY-1)fillrect(win,x,y,1,1);
     copylayer(win,1,0);
     msleep(1);
/* end of one step */
  }
 }
}

特に画面表示はあまり効率のよいプログラムになっていないが、これで動くのでよいことにする。また、このように1ステップごとに表示していると計算の進行が非常に遅くなる。成長の様子を加速して見たければ、粒子が静止するまで画面表示しないことにしたほうがよい。

粒子のスタート位置を最上段ではなくその一段下にしてあるのは、プログラム上のテクニカルな理由で、こうしないとs[x][y+1]が存在しないセルを指す場合が生じてしまう。なお、粒子が最上段に上がってしまったら、その時点でその粒子は捨てて、次の粒子に移るようにしている（したがって、多くの粒子が無駄になる）。 

     x = (x+NX)%NX;
というのはかなりトリッキーだが、周期的境界条件の処理である。もちろんここはif文でもかまわないのだが、このやりかたはその下の条件判断部分で 
  s[(x+NX-1)%NX][y+1] == 1
というように応用できる。なぜNXを足すのか考えてみよ。 

   for(;;){
という部分、何やってるのという感じだが、これは「永久ループ」である。break文が実行されない限り、いつまででもループを回り続けるという危険なプログラムである

上のプログラムは粒子が完全にランダムウォークするので、なかなか着地してくれない。横にはランダムウォークするが縦には１ステップにつき１段ずつ落ちる（ランダムではない）という別のモデルを考えることもできるだろう。その場合、移動部分は 

     x += (random()%3-1);
     x = (x+NX)%NX;
     y -= 1;
とすればよい。ただし、このふたつは別のモデルなので振る舞いも微妙に異なる
次へ