# 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()