# DIFFUSION LIMITED AGGREGATION # nacteme grafickou knihovnu a numericku knihovnu import matplotlib as mpl from matplotlib import pyplot as plt import numpy as np # definujeme rozmer mrizky a maximalni pocet castic roz = 200 mass_max = 8000 # pole pro spiralni 1D index, bezpecne vetsi nez je potreba roz_1d = 12*roz**2 # pole, do kterzch budeme ukladat prepoctene kartezske souradnice # definujeme vektory pomahajici pri prepoctu z hexa na ctvercovou mrizku xs = np.zeros(mass_max,dtype=float) ys = np.zeros(mass_max,dtype=float) xhelp = [np.sqrt(3.)/2.,0.,-np.sqrt(3.)/2.,-np.sqrt(3.)/2.,0.,np.sqrt(3.)/2.] yhelp = [-0.5,-1.,-0.5,0.5,1.,0.5] # pole nul, do ktereho budeme ukladat stavy bunek pole = np.zeros(roz_1d,dtype=np.int) # inicializujeme bunky, hmotnost a nejvetsi vzdalenost od stredu pole[0] = 1 mass = 1 rmax = 0 while mass < mass_max: # pocatecni vzdalenost castice r = rmax+1 # pricist k hranici v radialnim smeru, za kterou castice nesmi utect rkill = rmax # vyber nahodny prvek z sestiuhelniku o rozmeru r i = int(6*r*np.random.random() + 1) nn = 0 # pohyb difundujici castice while nn == 0: i_full = i+3*r*(r-1) # pozice v poli 'pole' # nejsme-li ve vrcholu sestiuhelniku if (i % r) != 0: # nalezneme sousedy if i == 1: ileft = 6*r idownleft = 3*r*(r-1) else: ileft = i-1 idownleft = (i*(r-1))//r + 3*(r-1)*(r-2) # dvojite lomitko vraci floor pole_nn = [i+1+3*r*(r-1),ileft+3*r*(r-1),idownleft, (i*(r-1))//r+3*(r-1)*(r-2)+1, (i*(r+1))//r+3*r*(r+1), (i*(r+1))//r+3*r*(r+1)+1] obsazeno = np.where(pole[pole_nn] == 1) nn = len(obsazeno[0]) # je-li nejaky soused obsazen, prilep se if nn > 0: pole[i_full] = 1 rmax = max([rmax,r]) mass += 1 # jinak difunduj dal else: step = int(np.random.random()*6) # radialni sestup if (step == 2 or step == 3): r -= 1 # radialni vzestup elif (step == 4 or step == 5): r += 1 i = pole_nn[step] - 3*r*(r-1) # pokud jsme ve vrcholu sestiuhelniku else: if i == 6*r: iright = 1 iupright = 1 else: iright = i+1 iupright = i*(r+1)//r+1 if i == 1: ileft = 6 else: ileft = i-1 pole_nn = [iright+3*r*(r-1),ileft+3*r*(r-1), i*(r-1)//r+3*(r-1)*(r-2),i*(r+1)//r+3*r*(r+1), iupright+3*r*(r+1),i*(r+1)//r-1+3*r*(r+1)] obsazeno = np.where(pole[pole_nn] == 1) nn = len(obsazeno[0]) if nn > 0: pole[i_full] = 1 rmax = max([rmax,r]) mass += 1 else: step = int(np.random.random()*6) if step == 2: r -= 1 elif ((step == 3 or step == 4) or step == 5): r += 1 i = pole_nn[step] - 3*r*(r-1) # pokud castice utekla, zacneme znovu if r > (rmax+rkill): r = rmax+1 i = int(6*r*np.random.random()+1) continue x0 = r*np.sin(2.*np.pi/(6.*r)*(i - (i % r))) y0 = r*np.cos(2.*np.pi/(6.*r)*(i - (i % r))) x = x0 + (i % r)*xhelp[min(int(i//r),5)] y = y0 + (i % r)*yhelp[min(int(i//r),5)] xs[mass-1] = x ys[mass-1] = y # vykreslime castice pomoci sestiuhelnikovych symbolu dlaplot = plt.scatter(xs,ys,c='k',s=30000/(roz**2),marker=(6,0,30)) plt.xlim(-roz, roz) plt.ylim(-roz, roz) plt.axes().set_aspect('equal') plt.xlabel('x') plt.ylabel('y') plt.show()