Help improve code efficiency

[benbrown03]

if anyone could help me improve the efficiency of this code it would be greatly appreciated. I would like it to only loop over the newly contaminated cells at each time step instead of looping over the whole grid. Any other ideas welcome.

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import colors
import random
import os


print("=== Cellular Automata for Radioactive Waste Disposal ===")

# ======= User Inputs

L = 20.0 # km

numCells = 100 # This needs to remain the same for!

totalTime = (40 * 12) # number of time steps (in months)

print_every_n_steps = 480 # don't make this too small!
print_time_every = 200

source_x = 0.2 # x - coordinate location of source of radio active waste
source_y = 0.2 # y - coordinate location of course  of radio active waste

city_x = 12
city_y = 16
cell_city_x = int(numCells * (city_x / L))
cell_city_y = int(numCells * (city_y / L))

verbosity = 0

#Permeability Maximum Variable

# *** Model Developement 3 ***
Kmax = 3.5#Permeability Maximum Variable
#K = np.loadtxt("perm.txt", dtype='f', delimiter=' ') 
Ktmp = np.loadtxt("perm.txt")

K=Ktmp.astype(np.float)

cmap = mpl.cm.jet
norm = mpl.colors.Normalize(vmin=0, vmax = 1)
fig, ax = plt.subplots()
ax.imshow(K, cmap=cmap, norm=norm)

ax.grid(which='major', axis='both', linestyle='-', color='k', linewidth=2)
ax.axes.get_xaxis().set_visible(False)
ax.axes.get_yaxis().set_visible(False)
results_dir = os.path.join(os.getcwd(),'solutions/')
if not os.path.isdir(results_dir):
  os.makedirs(results_dir)
    
plt.savefig(results_dir + "perm.png")

cell_source_x = int(numCells * (source_x / L))
cell_source_y = int(numCells * (source_y / L))
count = 0

for n in range(0, 20):
  C = np.zeros((numCells	,numCells ))

  C[cell_source_x,	cell_source_y] = 1; #1 is contaminated
  iteration = str(n+1)
  print("Iteration:" + iteration)
  #C[cell_city_x][cell_city_y] = 1 #will show city 

  for t in range(0,totalTime+1):
    
    if (t % print_time_every == 0) : 
      print("Time = " + str(t))

    Cnew = np.zeros((numCells	,numCells ))

    for i in range(0,numCells):

      for j in range(0,numCells):

        if (C[i][j] == 0) :#and i+1 j+1 ==0 etc

          if(verbosity > 0):
            print("Cell " + str(i) + " , " + str(j))

          # *** Model Developement 1 **
          #sij = contaminated neighbours
          
          sij = 0
          
  #using 'if' and 'elif' to assess the 9 different cases for a given cell (i,j)

          if(i!=0 and i!=numCells-1 and j!=0 and j!=numCells-1):
            sij= C[i][j+1] + C[i][j-1] + C[i-1][j] + C[i-1][j-1] + C[i-1][j+1] + C[i+1][j] + C[i+1][j-1] + C[i+1][j+1]
          
          elif(i==0 and j!=0 and j!=numCells-1):#mid left
            sij= C[i][j+1] + C[i][j-1] + C[i+1][j+1] + C[i+1][j] + C[i+1][j-1]

          elif(i==numCells-1 and j!=0 and j!=numCells-1):#mid right
            sij= C[i][j+1] + C[i][j-1] + C[i-1][j+1] + C[i-1][j] + C[numCells-1][numCells-1]
          
          elif(i!=0 and i!=numCells-1 and j==numCells-1):#mid top
            sij= C[i-1][j] + C[i-1][j-1] + C[i][j-1] + C[i+1][j-1] + C[i+1][j]
            
          elif(i!=0 and i!=numCells-1 and j==0):#mid bottom
            sij= C[i-1][j] + C[i-1][j+1] + C[i][j+1] + C[i+1][j+1] + C[i+1][j]

          elif(i == 0 and j==0) : #bottom left
            sij = C[0][1] + C[1][0] + C[1][1]
          
          elif(i==0 and j==numCells-1):#top left
            sij = C[1][numCells-1] + C[0][numCells-2] + C[1][numCells-2]
          
          elif(i==numCells-1 and j==0):#bottom right
            sij = C[numCells-1][1] + C[numCells-2][0] + C[numCells-2][1]
          
          else:#top right
            sij = C[numCells-1][numCells-2] + C[numCells-2][numCells-1] + C[numCells-2][numCells-2]
          
          
          # *** Model Developement 2 ***
          p = (sij/8)*((K[i][j])/Kmax) #gives probability
          n = random.uniform(0,1)
          if p>n:
            Cnew[i][j]=1#make this =p to give range of intensity?


        else :

          Cnew[i][j] = 0
            
    
    C = Cnew; # Set new contamination Matrix to old C.
    #this is how it goes through, adding on each time
    # Plot grid and save to file

    if (t % print_every_n_steps == 0) : 
      cmap = mpl.cm.jet
      norm = mpl.colors.Normalize(vmin=0, vmax = 1)
      fig, ax = plt.subplots()
      ax.imshow(C, cmap=cmap, norm=norm)

      ax.grid(which='major', axis='both', linestyle='-', color='k', linewidth=2)
      ax.axes.get_xaxis().set_visible(False)
      ax.axes.get_yaxis().set_visible(False)
      results_dir = os.path.join(os.getcwd(),'solutions/')
      if not os.path.isdir(results_dir):
        os.makedirs(results_dir)
      
      plt.savefig(results_dir + "Run:" + (iteration) + "_solution_time_" + str(t) + ".png")

  city = C[cell_city_x][cell_city_y]    
  if city == 1:
    print("City Dead")
    count += 1
  else:
    print("City Alive")
  print(count)