From GeoMod
from visual import *
from raster_map import *
''' set dimensions of grid'''
nrows = 41
ncols = 41
'''set parameters'''
dx = 100000.0 #meters
dz = 32000.0 #meters
dt_E = 1000 #year
K_E = 50000.0 #meters per year
S_E = K_E * dt_E / (dx*dx)
'''Determine if S exceeds the stability threshold'''
#if S > 0.0:
#exit
'''Initial conditions: create elevation array'''
Z = cr_top = raster_map(dz*ones((nrows,ncols)), dx=dx) #elevation array
color_scale = color_map(32000.0, 33000.0, 100.0, colmin=vector(1,0,0), colmax=vector(0,1,0))
Z.data[20,20]=dz+2000.0 #high elevation near center
'''Boundary conditions'''
Z.data[:,0]=dz #South
Z.data[:,nrows-1]=dz #North
Z.data[0,:]=dz #West
Z.data[ncols-1,:]=dz #East
'''Draw map of elevation'''
print Z.data
Z.line_3d(color_map=color_scale,scale=100.0)
nsteps = 1000
for nt in range(1, nsteps+1):
Z_new = Z.data[1:nrows-1, 1:ncols-1] + S * (Z.data[2:nrows, 1:ncols-1] - 4 * Z.data[1:nrows-1, 1:ncols-1] +
Z.data[0:nrows-2, 1:ncols-1] + Z.data[1:nrows-1, 2:ncols] +
Z.data[1:nrows-1, 0:ncols-2])
'''update'''
Z.data[1:nrows-1,1:ncols-1] = Z_new
Z.update_line_3d()
print nt, Z_new,
