# -*- coding: utf-8 -*-
"""
Created on Tue Apr  6 20:45:49 2021
Example of using a random walk CA to mimic heat flow acroos a conductivity boundary
@author: neil
"""

import matplotlib.pyplot as plt
import numpy as np

size = 10
rightBC = 0
leftBC  = 1000

# create a region, use integers, and make a new and an old region
r    = np.ones([size,size])
rnew = np.zeros([size])
rold = np.zeros([size,size])
r[:,0] = leftBC
r[:,-1]= rightBC

# this is a bit of an odd iteration, instead of just iterating on a single vector
# we iterate over a matrix that is as deep as the vector, each row is the next 
# iteration. When the bottom of the matrix is reached, we just move to the top
# of the matrix and continue.  This allows viewing the iteration progress as a
# 2D plot [commented out below]
iters =10000
while iters > 0:
    iters -=1
    for i in range(1,size):
        ran = np.random.uniform()
        d   = r[i-1,0]*.5
        rnew[1] = ran*d
    
        for j in range(1,size-1):
            ran = np.random.uniform()
            d   = r[i-1,j]*.5
            rnew[j]   = rnew[j] - d
            rnew[j-1] = rnew[j-1] + ran*d
            rnew[j+1] = rnew[j+1] + (1-ran)*d
        
        r[i,1:-1] = .5*r[i-1,1:-1]+.25*rnew[1:-1] + .5*rold[i-1,1:-1]
        rnew = np.zeros([size])
    for j in range(0,size):
        r[0,j] = np.sum(r[:,j])/size
    rold = .8*rold + .2*r

r = rold.copy()            
x = [0,2,4,6,7,8,9,10,11,12]
# a contour plot.  Here we make the contours stand out as black
fig4 = plt.figure()
ax4  = fig4.add_subplot(111)
ax4.plot(x,r[0,:],'r-*')
ax4.plot([6,6],[0,1000],"g:")
ax4.text(3,550,"K=1")
ax4.text(9,550,"K=2")
ax4.set_ylabel("Temperature")
ax4.set_xlabel("Distance")
ax4.set_title("CA temperature field with conductivity contrast")
# cs = ax4.contour(r,colors ='k')     # contour just plots the contours
# ax4.contourf(r)                     # contourf fills the contours with colors
# ax4.clabel(cs,inline=1,fontsize=12) # clabel puts labels on the contours
# ax4.invert_yaxis()
# ax4.set_ylabel("Depth")
# ax4.set_xlabel("Width")
# ax4.set_title("Temperature field as a contour")