Pythonで弾道軌道凝集ーお餅が空から降ってくるー

弾道軌道凝集とは

習字紙に墨汁を一滴垂らすと、黒いシミがじわじわと広がっていきます。このとき、墨汁によって黒くなった場所のぎざぎざとした境界線は、少しづつ形を変えながら広がっていきます。これを「成長する荒れた界面」といいます。
この「成長する荒れた界面」を計算機シミュレーションで模するために作られたモデルの1つが弾道軌道凝集です。

まず粒子(お餅)が空から格子上の地面(お皿の上)に向かって真っすぐ降ってきます。その粒子は地面に着くとそのままそこに落ち着きます。すでに粒子が置かれた場所に落ちるとその上に重なります。また、他の粒子の上にいなくても真横に他の粒子がいるとそこにくっついてしまいます。

これを繰り返していくことで、密な部分と疎な部分ができフラクタルのような模様が表れます。

案外すぐできる

初期位置を選んで、真横と真下に粒子があるかチェックしながら落とすだけなので案外すぐできました。

import numpy as np
import random
import matplotlib.pyplot as plt
import time

class BTA():
    def __init__(self):
        self.L = 500
        self.field = np.zeros([self.L,self.L], dtype=int)
        
    def init(self):
        ini_list = [[i,self.L-1] for i in range(self.L)]
        return random.choice(ini_list)
    
    def contact(self):
        x,y = self.init()
        for i in range(y+1):
            y_i = y - i
            #print(x,y_i)
            if x == 0:
                if self.field[x,y_i-1]==1 or 
                         self.field[self.L-1,y_i]==1 or self.field[1,y_i]==1:
                    self.field[x,y_i] = 1
                    break
            elif x == self.L-1:
                if self.field[x,y_i-1]==1 or 
                         self.field[x-1,y_i]==1 or self.field[0,y_i]==1:
                    self.field[x,y_i] = 1
                    break
            else:                
                if self.field[x,y_i-1]==1 or 
                         self.field[x-1,y_i]==1 or self.field[x+1,y_i]==1:
                    self.field[x,y_i] = 1
                    break
            if y_i == 0:
                #print("bottom")
                self.field[x,0] = 1
                break
            
    def procedure(self):
        for i in range(self.L**2):
            self.contact()
            
        plot_arr = np.array(np.where(self.field == 1))
        plt.figure(figsize=(8,8))
        plt.ylim(-1,self.L)
        plt.plot(*plot_arr,".",ms=1)
        #plt.show()
        plt.savefig("Ballistic_trajectory_aggregation.png", 
                           bbox_inches="tight", pad_inches=0.0)

bta = BTA()
bta.procedure()