import sys
import numpy as np
import math
import random
import matplotlib.pyplot as plt

def kernel1():
    kernel = []
    vec1 = [0.05, 0.2, 0.05]
    vec2 = [0.20,-1.0, 0.20]
    vec3 = [0.05, 0.2, 0.05]
    kernel.append(vec1)
    kernel.append(vec2)
    kernel.append(vec3)
    
    return kernel
    
kernel = kernel1()
print("kernel", kernel)

kernel [[0.05, 0.2, 0.05], [0.2, -1.0, 0.2], [0.05, 0.2, 0.05]]

def laplace(x, y, L, M, kernel):
    laplace = 0
    
    for i in range(3):
        for j in range(3):
            a = x + i - 1
            b = y + j - 1
            if (a >= 0 and b >= 0 and a < L and b < L):
                val = M[a][b] * kernel[i][j]
                laplace += val
    return laplace

def init_grid(val, x0, x1, y0, y1, M):
    for x in range(x0,x1,1):
        for y in range(y0,y1,1):
            M[x][y] = val
            
def init_grid_random(val, x0, x1, y0, y1, M):
    for x in range(x0,x1,1):
        for y in range(y0,y1,1):
            M[x][y] = val*random.uniform(0, 1)

def copyM(L, M):
    Copy = np.zeros((L,L))
    
    for x in range(L):
        for y in range(L):
            Copy[x][y] = M[x][y]
    return Copy

def brusselator(L, T, Da, Db, ka, kb, kernel, A, B):
    inc_t = 0.05
    inc_x = 0.5
    
    #pattern formation conditions
    if kb >= 1 + ka*ka:
        print("kb < 1+ka^2  failed", kb, 1+ka*ka)
        sys.exit()
    if kb <= (1 + ka*np.sqrt(Da/Db))*(1 + ka*np.sqrt(Da/Db)):
        print("kb > (1 + ka*sqrt(DA/DB)^2 failed", kb, (1 + ka*np.sqrt(Da/Db))*(1 + ka*np.sqrt(Da/Db)))
        sys.exit()
 
    t = 0
    n = 0
    while(t < T):
        thisA = copyM(L, A)
        thisB = copyM(L, B)
        
        for x in range(L):
            for y in range(L):
                Aval = thisA[x][y]
                Bval = thisB[x][y]
                
                # the reaction term
                addA  = ka - (1+kb)*Aval + Aval*Aval*Bval
                addB  = kb*Aval - Aval*Aval*Bval

                # the diffusion term
                addA += inc_x*DA*laplace(x, y, L, thisA, kernel)
                addB += inc_x*DB*laplace(x, y, L, thisB, kernel)
                                   
                if Aval+addA*inc_t < 0:
                    break
                if Bval+addB*inc_t < 0:
                    break

                A[x][y] += addA*inc_t
                B[x][y] += addB*inc_t
        
        
        if n % 50 == 0:
            print("time", t)
            plt.matshow(A)
            plt.colorbar()
            plt.show()


        t += inc_t
        n += 1
        
        
T = 500
L = 100

# My setup
DA = 5
DB = 12
ka = 3
kb = 9.9
A = np.zeros((L,L))
B = np.zeros((L,L))
init_grid_random(ka, 0, L-1, 0, L-1, A)
brusselator(L, T, DA, DB, ka, kb, kernel, A, B)

time 0

time 2.499999999999999

time 4.99999999999999

time 7.499999999999981

time 10.000000000000007

time 12.500000000000043

time 15.000000000000078

time 17.500000000000114

time 20.00000000000015

time 22.500000000000185

DA = 2.8
DB = 22.4
ka = 1.5
kb = 2.4
A = np.zeros((L,L))
B = np.zeros((L,L))
init_grid(1, 48, 52, 48, 52, A)
brusselator(L, T, DA, DB, ka, kb, kernel, A, B)

time 0

time 2.499999999999999

time 4.99999999999999

time 7.499999999999981

time 10.000000000000007

time 12.500000000000043

time 15.000000000000078

time 17.500000000000114

time 20.00000000000015

time 22.500000000000185