# these are code snipets that perform one time step in a lattice gas CA
# although the CA is on a hexagonal grid, we code it as a cubic latice, with 2 'extra' connections on the diagonals
# this requires keeping track of even and odd rows (cols are not changed)

# for simplicity we model a rectangular region, in addition we require even number nodes in both x and y

import numpy as np

lenx = 1000
leny = 1000
# fill cell with a random salting of integers, 0 to 63
cell = np.random.randint(0,64,size=(leny,lenx))                                # current cell grid
celln= np.zeros_like(cell)                                                     # next time step grid

rot = [9,18,36,9]                                                              # used in the randomizing step

# take one time step over the interior cells
# this encodes the entire 'transition table', by using binary numbers
for j in range(1,lenx-1):
    for i in range(2,leny-1,2):   # even lines
        celln[i,j] = ( cell[i,j+1]&8 + cell[i-1,j]&16 + cell[i-1,j-1]&32 +     # this moves all 6 surounding particles into the central cell
                       cell[i,j-1]&1 + cell[i+1,j-1]&2 + cell[i+1,j]&4  )
        if   celln[i,j] == 9:                                                  # these 3 are the 2 particle interactions
            celln[i,j] = rot[np.random.randint(1,3)]
        elif celln[i,j] == 18:
            celln[i,j] = rot[np.random.randint(2,4)]
        elif celln[i,j] == 36:
            celln[i,j] = rot[np.random.randint(0,2)]
        elif celln[i,j] == 21:                                                 # these are the 3 particle rotational collisions
            celln[i,j] = 42
        elif celln[i,j] == 42:
            celln[i,j] = 21
                    
    for i in range(1,leny-1,2):   # odd lines        
        celln[i,j] = ( cell[i,j+1]&8 + cell[i-1,j+1]&16 + cell[i-1,j]&32 + 
                       cell[i,j-1]&1 + cell[i+1,j]&2 + cell[i+1,j+1]&4  )
        if   celln[i,j] == 9:
            celln[i,j] = rot[np.random.randint(1,3)]
        elif celln[i,j] == 18:
            celln[i,j] = rot[np.random.randint(2,4)]
        elif celln[i,j] == 36:
            celln[i,j] = rot[np.random.randint(0,2)]
        elif celln[i,j] == 21:
            celln[i,j] = 42
        elif celln[i,j] == 42:
            celln[i,j] = 21