DLA sticking prob figure 8.py

# -*- coding: utf-8 -*-
"""
Created on Tue Nov 22 16:09:36 2022

@author: ppyma4
"""
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(1234) #same random seed each time

N = 1000 #size of lattice
Np = 3000 #numer of particles
cluster = 4
xycentre = N//2

relkill = 2.3 #value required to set new kill radius
relinj = 2 #value to required to set new injection radius

sprob = np.array([80, 60, 40, 20]) #sticking probability, 80%, 60%, 40% and 20%
title = np.array(['80%', '60%', '40%', '20%'])

fig = plt.figure(figsize = (14, 9))

for k in range(cluster):
    
    R = 6 #initial injection radius
    kill = 8 #initila kill radius

    lattice = np.zeros((N, N)) #new lattice
    lattice[xycentre, xycentre] = 1 #particle at center of lattice[x, y]
    
    xarray = ([])
    yarray = ([])
    xarray = np.append(xarray, xycentre)
    yarray = np.append(yarray , xycentre)
    
    for i in range(Np): #particle
        n = 0 #start while loop
    
        randomtheta = np.random.uniform(0, 2*np.pi) #spawns random particle on injection circle
        xs = round(R*np.cos(randomtheta)) + xycentre
        ys = round(R*np.sin(randomtheta)) + xycentre

        while n == 0:   
            random = np.random.randint(0, 4)
            if random == 0: #random walk of particle
                xs = xs + 1
                xs = xs*(lattice[xs, ys] == 0) + (xs-1)*(lattice[xs, ys] == 1)#prevent particles moving inside one another by removing step if it has moved into one
            elif random == 1:
                ys = ys + 1
                ys = ys*(lattice[xs, ys] == 0) + (ys-1)*(lattice[xs, ys] == 1)
            elif random == 2:
                xs = xs - 1
                xs = xs*(lattice[xs, ys] == 0) + (xs+1)*(lattice[xs, ys] == 1)
            elif random == 3:
                ys = ys - 1
                ys = ys*(lattice[xs, ys] == 0) + (ys+1)*(lattice[xs, ys] == 1)

            position = (xs-xycentre)**2 + (ys-xycentre)**2 #not rooted as makes code run faster, kill is squared to compensate this

            if (position >= kill**2): #reinject particle back into injection radius if touches kill circle
                randomtheta = np.random.uniform(0, 2*np.pi)
                xs = round(R*np.cos(randomtheta)) + xycentre
                ys = round(R*np.sin(randomtheta)) + xycentre
            elif lattice[xs+1, ys] == 1 or lattice[xs-1, ys] == 1 or lattice[xs, ys+1] == 1 or lattice[xs, ys-1] == 1: #sticking of particle
                if np.random.randint(0,100) <= sprob[k]: #sticking probability
                    
                    lattice[xs, ys] = 1 #set point in lattice to equal one
                        
                    position =  np.sqrt(position)
                
                    xarray = np.append(xarray, xs)
                    yarray = np.append(yarray ,ys)
                
                    if R < position*relinj: #create new injection radius if stuck particle is furthest from the seed
                        R = position*relinj
                        kill = position*relkill
                    
                    n = 1 #end while loop     
    
    plt.subplot(2, 2, k+1)

    plt.plot(yarray, xarray, 'o', markersize = 0.85, color = 'white', markeredgecolor = 'black', markeredgewidth = 1, zorder = 0)
    plt.plot(yarray[0], xarray[0], 'o', markersize = 1, color = 'red', zorder = 10)
    
    plt.axis('equal')
    plt.xlabel('$y$', fontsize = 17)
    plt.ylabel('$x$', fontsize = 17)
    plt.tick_params(axis = 'both', labelsize = 17)
    plt.title(title[k], fontsize = 17)
    
plt.subplots_adjust(wspace = 0.31, hspace = 0.48)