Help improve code efficiency

[benbrown03]

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import colors
import random
import os
import collections
 
 
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
  contaminated = set()
  contaminated.add((cell_source_x, cell_source_y))
  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 ))
 
    neighbors = collections.defaultdict(int)
    for source_x, source_y in contaminated:
      for off_x, off_y in [(0, 1), (1, 1), (1, 0), (1, -1), (0, -1), (-1, -1), (-1, 0), (-1, 1)]:
        neighbor_x, neighbor_y = source_x + off_x, source_y + off_y
        if 0 <= neighbor_x < numCells and 0 <= neighbor_y <= numCells:
          neighbors[(neighbor_x, neighbor_y)] += 1
 
    for coordinates, count in neighbors.items():
      possible_x, possible_y = coordinates
      # *** Model Developement 2 ***
      p = (count/8)*((K[possible_x][possible_y])/Kmax) #gives probability
      n = random.uniform(0,1)
      if p>n:
        Cnew[possible_x][possible_y]=1#make this =p to give range of intensity?
        contaminated.add(coordinates)
 
 
 
    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)

that looks to be loads faster thank you so much for your help
i get this error come up after running it a few times

Traceback (most recent call last):
File "main.py", line 90, in <module>
Cnew[possible_x][possible_y]=1#make this =p to give range of intensity?
IndexError: index 100 is out of bounds for axis 0 with size 100