import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors as c
import csv
import time

plt.style.use('_mpl-gallery-nogrid')


#Data makes up a grid with grid cell size (total length)/(len(elevData) and 
#grid dimensions len(elevData) x len(elevData[0])
#get actual elevations from file
'''
results = []
with open("subRegion2-5406.csv1801x1286") as csvfile:
	reader = csv.reader(csvfile, quoting=csv.QUOTE_NONNUMERIC) #change contents to floats
	for row in reader: # each row is a list
		results.append(row)
elevData = results
'''

#Array for Data:
elevData = np.array([[None]*5]*5) # in m
rainfallData = np.array([[None]*len(elevData[0])]*len(elevData)) #in µm (just due to odd round off errors) 
absorptionData = np.array([[None]*len(elevData[0])]*len(elevData)) # in µm/min
absorptionDataStatic = np.array([[None]*len(elevData[0])]*len(elevData)) # in µm/min
totalWaterAbsorbed = 0 # in m^3

#tile dimension
#tileDim = 133.5826*1000/5406 # in m 
#(full interval is 133.5826km wide there are 10812 pieces in full file so 10812/2=5406 pieces so tileDim = 133.5826*1000/5406)
#tileDim=5.0/float(len(elevData)) # in m
tileDim=10 # in m

rain=float(input("How many centimeters of rain should fall?"))

#filling the array with some random math to make it more fun

'''
i=0
while i < len(elevData):
	k=0
	while(k< len(elevData[0])):
		elevData[i][k] = math.trunc((20*math.cos(k*math.pi/4)**2+40*i+0.5*i**2+k*1.3))
		rainfallData[i][k] = (1.5*math.sin(k)*math.sin(i) + 1.5)*rain/3 #in cm over entire interval (measures total rainfall)
		absorptionData[i][k] = 5/(10+elevData[i][k]/15) #in cm/min
		absorptionDataStatic[i][k] = 5/(10+elevData[i][k]/15)  #in cm/min
		k+=1
	i+=1


i=0
while(i < len(elevData)):
	k=0
	while(k< len(elevData[0])):
		#elevData[i][k] = -5*math.exp(-i/4+2)+1500+k
		rainfallData[i][k] = rain #in cm over entire interval (measures total rainfall)
		absorptionData[i][k] = 0 #5/(10+elevData[i][k]/15) #in cm/min
		absorptionDataStatic[i][k] = 0 #5/(10+elevData[i][k]/15)  #in cm/min
		k+=1
	i+=1


#cool filling
i=0
while(i < len(elevData)):
	k=0
	while k < len(elevData[0]):
		r=(6*i/(len(elevData)-1)-3) # gives i on the scale -3 to 3
		f=(6*k/(len(elevData[0])-1)-3) #similarly gives k on the scale -3 to 3
		elevData[i][k] = (math.exp(-1*(r*r*f*f)) - math.exp(-1*(r*r+f*f))+1)*50 #offset to make absorption more consistent
		rainfallData[i][k] = rain #in µm over entire interval (measures total rainfall)
		absorptionData[i][k] = 0 #in cm/min
		absorptionDataStatic[i][k] = 0 #in cm/min
		k+=1
	i+=1

'''
i=0
while(i < len(elevData)):
	k=0
	while k < len(elevData[0]):
		r=(5*i/(len(elevData)-1)) #gives i on a 0 to 5 range
		f=(5*k/(len(elevData[0])-1)) #gives k on a 0 to 5 range
		elevData[i][k] = .2*r+.2*f#+3000 #offset to make absorption more consistent
		rainfallData[i][k] = rain #in µm over entire interval (measures total rainfall)
		absorptionData[i][k] = 0 #in cm/min
		absorptionDataStatic[i][k] = 0	#in cm/min
		k+=1
	i+=1
'''
#filling for accurate ABQ Area
i=0
while(i<len(elevData)):
	k=0
	while(k<len(elevData[0])):
		rainfallData[i][k] = rain
		if(elevData[i][k]<=1440):
			absorptionData[i][k] = .02413 #B HSG Rating in cm/min
		elif(elevData[i][k]<=1520):
			absorptionData[i][k] = .06011333 #A HSG Rating in cm/min
		elif(elevData[i][k]<=1630):
			absorptionData[i][k] = .013335 #Average of B and C Rating in cm/min
		else:
			absorptionData[i][k] = .00254 # C HSG Rating in cm/min
		#if within heavily urbanized area
		#if enitre interval 
		#if coords are within 35-35.2 and -106.3--105.9 of whole interval 34.4-35.6 -106.7--105.5
		#if i/(len(elevData)-1) is within 1/3 and 2/3
		#if k/(len(elevData[0])-1) is within 1/3 and 1/2
		if(i/(len(elevData)-1)>=1/3 and i/(len(elevData)-1)<=2/3):
			if(k/(len(elevData[0])-1)>=1/3 and k/(len(elevData[0])-1)<=1/2):
				if(elevData[i][k]>=1520): #if in river isn't city
					absorptionData[i][k] = 1/2*absorptionData[i][k] #cut absorption in half
		
		#if subset the these proportions for heavy urbanization
		if(i/(len(elevData)-1)<=.7767 and k/(len(elevData[0])-1)<=.5393):
			absorptionData[i][k] = 1/2*absorptionData[i][k] #cut absorption in half
		
		absorptionDataStatic[i][k] = absorptionData[i][k] #in cm/min
		k+=1
	i+=1
'''

#Initializing
printAbsorb = False
printExtras = False
printExtras2 = False
printResults = False
printFiles = False
percentGoal = 1

#choose which way to calculate
newMethod = True
oldMethod = not(newMethod)
recalculation = True

if(input("Initializing Data? (y/n)") == "y"):
	print("Elevation Data: ",elevData)
	print("Min: ",np.min(elevData))
	print("Max: ",np.max(elevData))
	print("Rainfall Data: ",rainfallData)
	



printExtensivityLevel = int(input("What level of extra data is requested? (0 to 2) (none to a lot)"))
if(printExtensivityLevel == 0):
	printExtras = False
	printExtras2 = False
elif(printExtensivityLevel == 1):
	printExtras = True
	printExtras2 = False
elif(printExtensivityLevel == 2):
	printExtras =  True
	printExtras2 = True
	
if(input("Print preliminary results? (y/n)")=="y"):
	printResults = True
else:
	printResults = False

if(input("Print files? (y/n)")=="y"):
	printFiles = True
else:
	printFiles = False


timescale = int(input("How long should the rainfall be spread out over? (minutes)"))
addition = int(input("How long should the model run after the rain has stopped? (minutes)"))

timescale += addition

#One limitation of the my algorithmn is that it only allows water to travel one pixel per tick at the maximum which could prove quite misrepresentative when having huge numbers of pixels and limited ticks (thus the need to really kick up the tick number when using large datasets)
accuracy = int(input("Into how many pieces should the program break up the interval? (more pieces --> more accurate) (pieces > total min of simulation recommended)"))

#Potentially precalculate water flow per cell here

currentSurfaceWater = np.array([[0.0]*len(elevData[0])]*len(elevData))
futureSurfaceWater = np.array([[0.0]*len(elevData[0])]*len(elevData)) #to be added next tick

#Initialized for use later in loops
problemAreas = np.array([[0]*len(elevData[0])]*len(elevData))
erosionAreas = np.array([[0]*len(elevData[0])]*len(elevData))
waterMovement = np.array([[0.0]*len(elevData[0])]*len(elevData))
waterDeficitArr = np.array([[0.0]*len(elevData[0])]*len(elevData))


stepSize = timescale/accuracy #for rain and absorption
print("Step size: ",stepSize, "min")
print("tileDim: ", tileDim, "m")


#To speed up calculation time the proportions of flowig water will be precalculated
flows = np.array([[[0.0]*len(elevData[0])]*len(elevData)]*8)
flowsTemp = np.array([[[0.0]*len(elevData[0])]*len(elevData)]*8)

#the first set of brackets in flows gives the direction in which the flow is to be done eg: flows[0] gives the flows north
# 0 is north, up, or decreasing in i
# 1 is northeast, diagonally up right, or increasing k and decreasing i
# 2 is east, right, or increasing in k
# 3 is southeast, diagonally down right, or increasing k and i
# 4 is south, down, or increasing i
# 5 is southwest, diagonally down left, or decreasing k and increasing i
# 6 is west, left, or decreasing k
# 7 is northwest, diagonally up left, or decreasing k and i

#important to note that these directions might not reflect the real north of the data
#they are really just to help visualize the directions of traversal and checking for this array.

	# Actual fluid dnyamics equations give (after isolated and integrated) that total volumetric flow of liquid on a plane through a rectange (l*d). Source: https://ocw.mit.edu/courses/12-090-introduction-to-fluid-motions-sediment-transport-and-current-generated-sedimentary-structures-fall-2006/72c0d3ee9656652db67315787b47c7a8_ch4.pdf
	# Kinda common knowledge ish (?)
	# is given by:
	# 2/3*d^3*9.8m/s^2*1g/cm^3*sin(α)/(1.308mPa*s)*l= 2/3*9.8/1.308*m*g/(s^3*cm^3*g/ms^2)*l=2/3*9.8/1.308*m^2/(s*cm^3)*(100cm/m)^3*l=2/3*9.8*100^3/1.308/sm*l
	# Which can be simplified to 
	#				4994903*sin(α)*d^3*l
	#				--------------------   gives cm^3/s
	#						s*m
	# This is an issue as my model works in depth not volume, so to fix this, dividing this change in volume by the area it affects (tileDim*tileDim),
	# gives the change in depth for the whole model. 
	#				4994903*sin(α)*d^3*l
	#				--------------------   
	#					s*m*tileDim^2
	# Finding l is also a bit strange: l could feasibly be tileDim for linearly adjacent tiles but what happens then for diagonals? 
	# To avoid counting borders as passing water from two different calculations the entire perimeter (4*tileDim) will be broken up among the adjacent 
	# tiles (numbering 8) so each tile will recieve 1/2*tileDim as its l to run water through and thus 1/2*tileDim = l can be substituted in:
	#				4994903*sin(α)*d^3*1/2*tileDim		2497451.5*sin(α)*d^3
	#				-----------------------------	=	-------------------- 
	#					   s*m*tileDim^2				   s*m*tileDim
	# Since however, d is in cm we should convert the meters in the bottom to cm so we get 24974.515 as the leading coefficient
	# Likewise we need to convert the tileDim which is in m to cm thus the coefficient is 249.74515
	# Since we are getting cm/s out of this equation we need to multiply by the stepSize with a conversion from s to min to get the real water depth change per step so the final coefficient is 14984.709

# flows gives the precalculable part of this: the 14984.709*sin(α)*stepSize
#												  ------------------------------
#															 tileDim
# Which will be consequently multiplied by the variable part during the looped section of the algorithm: the depth of water cubed 

startF = time.time()
i=0
while(i<len(flows[0])): #loops over the tiles
	k=0
	while(k<len(flows[0][0])):
		home = elevData[i][k]
		
		#consider northern tile					
		if(i!=0): #check if exists													 # - 0 -
			if(elevData[i-1][k] >= home):#conversion cm --> m						 # - X -
				flows[0][i][k] = 0 #if its higher just stop							 # - - -
			else:
				α=math.atan((home - elevData[i-1][k])/tileDim)
				flows[0][i][k]= ((14984.709)*math.sin(α)*stepSize)/tileDim
			if(flows[0][i][k]<(10**-6)):
				flows[0][i][k]=0

		#consider northeastern tile
		if(i!=0 and k!=len(elevData[0])-1): #check if exists						 # - - 1
			if(elevData[i-1][k+1] >= home):											 # - X -
				flows[1][i][k] = 0													 # - - -
			else:
				α=math.atan((home - elevData[i-1][k+1])/(math.sqrt(2)*tileDim)) 
				#multiplied by √2 because √2*tileDim is the true distance of the diagonal
				flows[1][i][k]= ((14984.709)*math.sin(α)*stepSize)/tileDim
			if(flows[1][i][k]<(10**-6)):
				flows[1][i][k]=0

		#consider eastern tile														 # - - -
		if(k!=len(elevData[0])-1): #exists?											 # - X 2
			if(elevData[i][k+1] >= home):											 # - - -
				flows[2][i][k] = 0
			else:
				α=math.atan((home - elevData[i][k+1])/tileDim)
				flows[2][i][k]= ((14984.709)*math.sin(α)*stepSize)/tileDim
			if(flows[2][i][k]<(10**-6)):
				flows[2][i][k]=0

		#consider southeastern tile 
		if(k!=len(elevData[0])-1 and i!=len(elevData)-1):							 # - - -
			if(elevData[i+1][k+1]>=home):											 # - X -
				flows[3][i][k] = 0													 # - - 3
			else:
				α=math.atan((home - elevData[i+1][k+1])/(math.sqrt(2)*tileDim))
				flows[3][i][k]= ((14984.709)*math.sin(α)*stepSize)/tileDim
			if(flows[3][i][k]<(10**-6)):
				flows[3][i][k]=0
				
		#consider southern tile 
		if(i!=len(elevData)-1):														 # - - -
			if(elevData[i+1][k]>=home):												 # - X -
				flows[4][i][k] = 0													 # - 4 -
			else:
				α=math.atan((home - elevData[i+1][k])/tileDim)
				flows[4][i][k]= ((14984.709)*math.sin(α)*stepSize)/tileDim
			if(flows[4][i][k]<(10**-6)):
				flows[4][i][k]=0

		#consider southwestern tile 
		if(i!=len(elevData)-1 and k!=0):											 # - - -
			if(elevData[i+1][k-1]>=home):											 # - X -
				flows[5][i][k] = 0													 # 5 - -
			else:
				α=math.atan((home - elevData[i+1][k-1])/(math.sqrt(2)*tileDim))
				flows[5][i][k]= ((14984.709)*math.sin(α)*stepSize)/tileDim
			if(flows[5][i][k]<(10**-6)):
				flows[5][i][k]=0

		#consider western tile 
		if(k!=0):																	 # - - -
			if(elevData[i][k-1]>=home):												 # 6 X -
				flows[6][i][k] = 0													 # - - -
			else:
				α=math.atan((home - elevData[i][k-1])/tileDim)
				flows[6][i][k]= ((14984.709)*math.sin(α)*stepSize)/tileDim
			if(flows[6][i][k]<(10**-6)):
				flows[6][i][k]=0

		#consider northwestern tile 
		if(k!=0 and i!=0):															 # 7 - -
			if(elevData[i-1][k-1]>=home):											 # - X -
				flows[7][i][k] = 0													 # - - -
			else:
				α=math.atan((home - elevData[i-1][k-1])/(math.sqrt(2)*tileDim))
				flows[7][i][k]= ((14984.709)*math.sin(α)*stepSize)/tileDim
			if(flows[7][i][k]<(10**-6)):
				flows[7][i][k]=0
				
		k+=1
	i+=1
		
t=0

while(t < timescale):
	start = time.time()
	if(t >= timescale - addition):
		#print("rain done")
		i=len(elevData)
	else:
		i=0

	#complete rainfall + absorption
	while(i<len(elevData)):
		k=0
		while(k<len(elevData[0])):
			#rainfall in amount & times according to stepSize and time raining
			surfaceWater = rainfallData[i][k]/(timescale-addition)*stepSize
			
			#Print the surface water before the absorption if extra data of level 2 is requested
			if(printExtras2):
				print("surfaceWater: ", surfaceWater)
			
			#absorption (again times stepSize)
			surfaceWater -= absorptionData[i][k]*stepSize
			totalWaterAbsorbed += absorptionData[i][k]*stepSize/100*tileDim**2
			
			#print the surface water after the absorption as well
			if(printExtras2):
				print("surfaceWater2: ", surfaceWater)
			
			
			#Store surface water additions to current
			currentSurfaceWater[i][k]+=surfaceWater
			
			#Don't want negative values for this
			if(surfaceWater<0):
				currentSurfaceWater[i][k] = 0

			if(printExtras2):
				print("Array: i: ",i,", k: ",k,", value: ",currentSurfaceWater[i][k])
			
			k+=1
		i+=1

	
	#if results are requested then print 

	if(printExtras):
		print("current surface water: ")
		print(currentSurfaceWater)
		


	#complete flow
	i=0
	while(i<len(elevData)):
		k=0
		while(k<len(elevData[0])):
			#original approach
			if(oldMethod):
				doTop=True
				doBottom=True
				doRight=True
				doLeft=True
				
				#print extras if wanted for debugging
				if(printExtras2):
					print("I and K:", i,", ", k)
					print(i)
					print(k)
					print(len(elevData))
					print(len(elevData[0]))
				
				#ensure no checking values off the edge
				if(i==0):
					doTop = False
				if(i==(len(elevData)-1)):
					doBottom = False
				if(k==0):
					doLeft = False
				if(k==len(elevData[0])-1):
					doRight = False
				
				#print extras if wanted for debugging
				if(printExtras2):
					print("Top",doTop)
					print("Bottom",doBottom)
					print("Left",doLeft)
					print("Right",doRight)
					
				#convert booleans into an array of coordinates to check
				
				#add middle tile
				tempArr = np.reshape(np.array([i,k]),(1,2))
				directionArr = np.array([], dtype=int)
				
				
				if(doTop):																 # - X -
					numRows = len(tempArr)												 # - - -
					tempArr = np.reshape(np.append(tempArr,[i-1,k]),(numRows+1,2))		 # - - -
					directionArr = np.append(directionArr,0)
					#print("N",tempArr)
					if(doLeft):															 # X - -
						numRows = len(tempArr)											 # - - -
						tempArr = np.reshape(np.append(tempArr,[i-1,k-1]),(numRows+1,2)) # - - -
						directionArr = np.append(directionArr,7)
						#print("NW",tempArr)
					if(doRight):														 # - - X
						numRows = len(tempArr)											 # - - -
						tempArr = np.reshape(np.append(tempArr,[i-1,k+1]),(numRows+1,2)) # - - -
						directionArr = np.append(directionArr,1)
						#print("NE",tempArr)
				if(doBottom):															 # - - -
					numRows = len(tempArr)												 # - - -
					tempArr = np.reshape(np.append(tempArr,[i+1,k]),(numRows+1,2))		 # - X -
					directionArr = np.append(directionArr,4)
					#print("S",tempArr)
					if(doLeft):															 # - - -
						numRows = len(tempArr)											 # - - -
						tempArr = np.reshape(np.append(tempArr,[i+1,k-1]),(numRows+1,2)) # X - -
						directionArr = np.append(directionArr,5)
						#print("SW",tempArr) 
					if(doRight):														 # - - -
						numRows = len(tempArr)											 # - - -
						tempArr = np.reshape(np.append(tempArr,[i+1,k+1]),(numRows+1,2)) # - - X
						directionArr = np.append(directionArr,3)
						#print("SE",tempArr)
				if(doRight):															 # - - -
					numRows = len(tempArr)												 # - - X
					tempArr = np.reshape(np.append(tempArr,[i,k+1]),(numRows+1,2))		 # - - -
					directionArr = np.append(directionArr,2)
					#print("E",tempArr)
				if(doLeft):																 # - - -
					numRows = len(tempArr)												 # X - -
					tempArr = np.reshape(np.append(tempArr,[i,k-1]),(numRows+1,2))		 # - - -
					directionArr = np.append(directionArr,6)
					#print("W",tempArr)
				
				#convert coordinates into values from elevData
				elevArray = np.array([])
				waterArray = np.array([])
				n=0
				
				while(n<len(tempArr)):
					elevArray = np.append(elevArray,elevData[tempArr[n][0]][tempArr[n][1]])
					waterArray = np.append(waterArray,currentSurfaceWater[tempArr[n][0]][tempArr[n][1]])
					n+=1
					
				#print if requested
				if(printExtras):
					print("tempArray: ", tempArr)
					print("elevArray: ",elevArray)
					print("waterArrray: ",waterArray)
					
				#determine flow vectors/distribution for water
				
				#This was the original approach:
					#Determine the two lowest sides and send the water accordingly
					#Decided against as it proves inacurrant for extremes
					#a local maximum height where water should really flow off in all directions would be simulated as just flowing off in two
					#a local min would also be sending water uphill given the two lowest considered --> Would need additional controls
					#Following approach is more natural for what water truely does
				'''
				lowest = elevArray[0]
				lowest2 = elevArray[0]
				lowestCoord = np.array([tempArr[0],tempArr[1]])
				lowestCoord2 = np.array([tempArr[0],tempArr[1]])
				n=0
				
				# code finds lowest in whole
				while(n<len(elevArray)):
					#print("elev: ",elevArray[n])
					if(elevArray[n] < lowest):
						#print("lowest: ",lowestCoord)
						lowestCoord2 = np.array([lowestCoord]).copy()
						lowestCoord = np.array([tempArr[2*n],tempArr[2*n+1]])
						lowest2 = lowest
						lowest = elevArray[n]
					elif(elevArray[n] < lowest2):
						lowest2 = elevArray[n]
						lowestCoord2 = np.array([tempArr[2*n],tempArr[2*n+1]])
						
					n+=1
				
				print("Lowest and Second lowest")
				print(lowest)
				print(lowest2)
				print("Their Coordinates, respectively:")
				print(lowestCoord)
				print(lowestCoord2)
				'''
				# ^ now outdated array dimensions --> won't run
				
				
				
				#Different approach: 
				#Determine which are higher than the center and send water proportionally to the percentage of toptal drop
				
				homeElev = elevArray[0]
				
				totalDeficit = 0
				realWaterDeficit = 0
				coordsArr = np.array([], dtype= int)
				#if the tile is higher or equal to home +currentSurfaceWater declinesArr will not be overwritten thus if they are all set to homeElev +currentSurfaceWater then if the tile is
				#higher or equal to and the declinesArr value is not changed then the resulting rate from the equation is 0
				declinesArr = np.array([homeElev+currentSurfaceWater[i][k]/100]*8, dtype= float)
				totalRate = 0
				waterRateArr = np.array([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0])
				
				#iterate through elevArray and pick out lowers
				#print(directionArr)
				#print(elevArray)
				n=1
				while(n<len(elevArray)):
					if(elevArray[n]+waterArray[n] < homeElev+currentSurfaceWater[i][k]/100):
						totalDeficit += (homeElev+currentSurfaceWater[i][k]/100-elevArray[n]-waterArray[n])
						numRows = len(coordsArr)
						declinesArr[directionArr[n-1]] = elevArray[n]+waterArray[n] #stores real elevation (sum of actual physical and water hieght) (since elevArray begins with homeElev and directionArr does not it accesses one behind it
						coordsArr = np.reshape(np.append(coordsArr,[tempArr[n][0],tempArr[n][1]]),(numRows+1,2))
					n+=1
				
				
				if(printExtras2):
					print("Coords Array: ", coordsArr)
					print("Declines Array: ", declinesArr)
					print("Total Deficit: ", totalDeficit)
					print("Home Elevation: ", homeElev)
					print("Elevation Array: ", elevArray)
					print("Coords Array length: ",len(coordsArr)) 
	
				#iterate through array of declines and gauge and send water from currentSurfaceWater
				n=0
				while(n<len(declinesArr)):
					'''
					#(homeElev-declinesArr[1])/totalDeficit% of the current water goes to coordsArr[1] with adjustment 
					#for amount of time	 and for actual decline (water goes down a steeper decline faster) old approach
					#simplistic assumptions
						#assumes all water can get through in 1 min
						#assumes that once the slope of a decline has reached -1 units down per one unit across flow is uninhibited
					if(stepSize <=1):#<-- simulation devolves when the step size is increased above one as flow is then inhibited
						if((homeElev - declinesArr[n])/(tileDim)<=1):
							additive = stepSize*((homeElev - declinesArr[n])/(2*tileDim))
						else:
							additive = stepSize
					elif((homeElev - declinesArr[n])/50<=1):
						additive = (homeElev - declinesArr[n])/50
					else:
						additive = 1 #essentially lets all water through
					'''
					
					#to get the slope angle for the following calculation we need to determine if it is diagonally from the point of interest
					if(n==0 or n==2 or n==4 or n==6):  #essentially if it is just North, East, South, or West it is not diagonal
						#arctan of Δelev/Δx
						α=math.atan((homeElev+currentSurfaceWater[i][k]/100 - declinesArr[n])/tileDim) #into m
					else:
						#if it is diagonal the distance that it is really covering is √2*tileDim
						α=math.atan((homeElev+currentSurfaceWater[i][k]/100 - declinesArr[n])/(math.sqrt(2)*tileDim)) #into m
						
					#See flows definition for info on formula
					
					if(declinesArr[n]<homeElev): #if lower just on its own the use this
						realWaterRatePerStep = ((14984.709)*math.sin(α))*((currentSurfaceWater[i][k])**3)/(tileDim)*stepSize
					else: #if its lower only because of the water level being higher then calculate fow only by the water that is higher.
						realWaterRatePerStep = ((14984.709)*math.sin(α))*((homeElev*100-declinesArr[n]*100+currentSurfaceWater[i][k])**3)/(tileDim)*stepSize #into centimeters
					#realWaterRatePerStep=additive*currentSurfaceWater[i][k]
	
					waterRateArr[n] = realWaterRatePerStep
					totalRate += realWaterRatePerStep
					if(printExtras2):
						print("water rate at n = ",n,": ",realWaterRatePerStep)
	
					
					n+=1
				if(printExtras):
					print("rate array: ", waterRateArr)
	
				#It is necessary to pull the actual implimentation of this calculation out of the loop because it finds the maximum rate at which the water can leave the square so after that if there is more water than can leave by this model then all is good but if it were inside the loop and it is very steep and there isn't much water then a majority of water could leave when maybe the other sides had the same water flow max or a higher max. This is really just to distribute the water evenly across the directions.
				
			if(newMethod):
				waterRateArr = np.array([])
				d = currentSurfaceWater[i][k]
				totalRate = 0
				realWaterDeficit = 0
	
				n=0
				while(n<8):
					waterRateArr = np.append(waterRateArr,flows[n][i][k]*(d**3))
					totalRate += waterRateArr[n]
					n+=1
				#print("Total Rate 1: ", totalRate)
				#print("water rate array 1: ", waterRateArr)
				if(recalculation):
					#See if recalculation is needed:
					home = elevData[i][k]
					homeWaterCM = currentSurfaceWater[i][k] #in cm
					realHomeElev = home + homeWaterCM/100 # in m
					#check northern tile
					if(i!=0): #if exists
						realNorthElev = elevData[i-1][k]+currentSurfaceWater[i-1][k]/100 # in m
						#print("realHomeElev - realNorthElev ",realHomeElev - realNorthElev)
						#print("i and k: ", i," ",k)
						if(elevData[i-1][k]>=home): #if north is taller
							if(realNorthElev<realHomeElev): #if north is shorter in real heights (heights with water added)
								#recalculate
								α=math.atan((realHomeElev - realNorthElev)/tileDim) #all in m
								waterRateArr[0] = ((14984.709)*math.sin(α))*((realHomeElev-realNorthElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[0] 1: ",waterRateArr[0])
							#if it's taller or in the middle with real heights then there is still no flow which normal model does fine
							#(by in the middle I mean when the real height is great enough to impede flow coming in but not great enough to make any of its own flow) 
							#See the next if statement for example
						elif(elevData[i-1][k]<home): #if north is shorter
							if(realNorthElev>home and realNorthElev<realHomeElev): #but if the water reaches a point where it stops some of the water from flowing over
								#essentially if the water level of the other reaches above the ground level of the other then whatever water it 'covers' won't flow
								#recalculate
								α=math.atan((realHomeElev - realNorthElev)/tileDim) # all in m
								waterRateArr[0] = ((14984.709)*math.sin(α))*((realHomeElev-realNorthElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[0 2: ",waterRateArr[0])
							elif(realNorthElev>=realHomeElev): #if its taller in real elevation
								#recalculate
								waterRateArr[0]=0
								if(printExtras):
									print("waterRateArr[0] 3: ",waterRateArr[0])
							#elif(realNorthElev<home): if it were simply below
								#don't recalculate as the normal model does this just fine
					#check northeastern tile
					if(i!=0 and k!=len(elevData[0])-1): #if exists
						realNortheastElev = elevData[i-1][k+1]+currentSurfaceWater[i-1][k+1]/100 # in m
						if(elevData[i-1][k+1]>=home): #if taller
							if(realNortheastElev<realHomeElev): #ifshorter in real heights (heights with water added)
								#recalculate
								α=math.atan((realHomeElev - realNortheastElev)/(math.sqrt(2)*tileDim)) #all in m
								waterRateArr[1] = ((14984.709)*math.sin(α))*((realHomeElev-realNortheastElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[1] 1: ",waterRateArr[1])
							#if it's taller or in the middle with real heights then there is still no flow which normal model does fine
							#(by in the middle I mean when the real height is great enough to impede flow coming in but not great enough to make any of its own flow) 
							#See the next if statement for example
						elif(elevData[i-1][k+1]<home): #if shorter
							if(realNortheastElev>home and realNortheastElev<realHomeElev): #but if the water reaches a point where it stops some of the water from flowing over
								#essentially if the water level of the other reaches above the ground level of the other then whatever water it 'covers' won't flow
								#recalculate
								α=math.atan((realHomeElev - realNortheastElev)/(math.sqrt(2)*tileDim)) # all in m
								waterRateArr[1] = ((14984.709)*math.sin(α))*((realHomeElev-realNortheastElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[1] 2: ",waterRateArr[1])
							elif(realNortheastElev>=realHomeElev): #if taller in real elevation
								#recalculate
								waterRateArr[1]=0
								if(printExtras):
									print("waterRateArr[1] 3: ",waterRateArr[1])
							#elif(realNortheastElev<home): if it were simply below
								#don't recalculate as the normal model does this just fine
					#check easterern tile
					if(k!=len(elevData[0])-1): #if exists
						realEastElev = elevData[i][k+1]+currentSurfaceWater[i][k+1]/100 # in m
						if(elevData[i][k+1]>=home): #if taller
							if(realEastElev<realHomeElev): #if shorter in real heights (heights with water added)
								#recalculate
								α=math.atan((realHomeElev - realEastElev)/tileDim) #all in m
								waterRateArr[2] = ((14984.709)*math.sin(α))*((realHomeElev-realEastElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[2] 1: ",waterRateArr[2])
							#if it's taller or in the middle with real heights then there is still no flow which normal model does fine
							#(by in the middle I mean when the real height is great enough to impede flow coming in but not great enough to make any of its own flow) 
							#See the next if statement for example
						elif(elevData[i][k+1]<home): #if shorter
							if(realEastElev>home and realEastElev<realHomeElev): #but if the water reaches a point where it stops some of the water from flowing over
								#essentially if the water level of the other reaches above the ground level of the other then whatever water it 'covers' won't flow
								#recalculate
								α=math.atan((realHomeElev - realEastElev)/tileDim) # all in m
								waterRateArr[2] = ((14984.709)*math.sin(α))*((realHomeElev-realEastElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[2] 2: ",waterRateArr[2])
							elif(realEastElev>=realHomeElev): #if taller in real elevation
								#recalculate
								waterRateArr[2]=0
								if(printExtras):
									print("waterRateArr[2] 3: ",waterRateArr[2])
							#elif(realEastElev<home): if it were simply below
								#don't recalculate as the normal model does this just fine
					#check southeastern tile
					if(i!=len(elevData)-1 and k!=len(elevData[0])-1): #if exists
						realSoutheastElev = elevData[i+1][k+1]+currentSurfaceWater[i+1][k+1]/100 # in m
						if(elevData[i+1][k+1]>=home): #if taller
							if(realSoutheastElev<realHomeElev): #if shorter in real heights (heights with water added)
								#recalculate
								α=math.atan((realHomeElev - realSoutheastElev)/(math.sqrt(2)*tileDim)) #all in m
								#print("realHomeElev: ",realHomeElev, " realSoutheastElev: ",realSoutheastElev)
								waterRateArr[3] = ((14984.709)*math.sin(α))*((realHomeElev-realSoutheastElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[3] 1: ",waterRateArr[3])
							#if it's taller or in the middle with real heights then there is still no flow which normal model does fine
							#(by in the middle I mean when the real height is great enough to impede flow coming in but not great enough to make any of its own flow) 
							#See the next if statement for example
						elif(elevData[i+1][k+1]<home): #if shorter
							if(realSoutheastElev>home and realSoutheastElev<realHomeElev): #but if the water reaches a point where it stops some of the water from flowing over
								#essentially if the water level of the other reaches above the ground level of the other then whatever water it 'covers' won't flow
								#recalculate
								α=math.atan((realHomeElev - realSoutheastElev)/(math.sqrt(2)*tileDim)) # all in m
								waterRateArr[3] = ((14984.709)*math.sin(α))*((realHomeElev-realSoutheastElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[3] 2: ",waterRateArr[3])
							elif(realSoutheastElev>=realHomeElev): #if taller in real elevation
								#recalculate
								waterRateArr[3]=0
								if(printExtras):
									print("waterRateArr[3] 3: ",waterRateArr[3])
							#elif(realSoutheastElev<home): if it were simply below
								#don't recalculate as the normal model does this just fine
					#check Southern tile
					if(i!=len(elevData)-1): #if exists
						realSouthElev = elevData[i+1][k]+currentSurfaceWater[i+1][k]/100 # in m
						if(elevData[i+1][k]>=home): #if taller
							if(realSouthElev<realHomeElev): #if shorter in real heights (heights with water added)
								#recalculate
								α=math.atan((realHomeElev - realSouthElev)/tileDim) #all in m
								waterRateArr[4] = ((14984.709)*math.sin(α))*((realHomeElev-realSouthElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[4] 1: ",waterRateArr[4])
							#if it's taller or in the middle with real heights then there is still no flow which normal model does fine
							#(by in the middle I mean when the real height is great enough to impede flow coming in but not great enough to make any of its own flow) 
							#See the next if statement for example
						elif(elevData[i+1][k]<home): #if north is shorter
							if(realSouthElev>home and realSouthElev<realHomeElev): #but if the water reaches a point where it stops some of the water from flowing over
								#essentially if the water level of the other reaches above the ground level of the other then whatever water it 'covers' won't flow
								#recalculate
								α=math.atan((realHomeElev - realSouthElev)/tileDim) # all in m
								waterRateArr[4] = ((14984.709)*math.sin(α))*((realHomeElev-realSouthElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[4] 2: ",waterRateArr[4])
							elif(realSouthElev>=realHomeElev): #if its taller in real elevation
								#recalculate
								waterRateArr[4]=0
								if(printExtras):
									print("waterRateArr[4] 3: ",waterRateArr[4])
							#elif(realSouthElev<home): if it were simply below
								#don't recalculate as the normal model does this just fine
					#check southwestern tile
					if(i!=len(elevData)-1 and k!=0): #if exists
						realSouthwestElev = elevData[i+1][k-1]+currentSurfaceWater[i+1][k-1]/100 # in m
						if(elevData[i+1][k-1]>=home): #if taller
							if(realSouthwestElev<realHomeElev): #if shorter in real heights (heights with water added)
								#recalculate
								α=math.atan((realHomeElev - realSouthwestElev)/(math.sqrt(2)*tileDim)) #all in m
								waterRateArr[5] = ((14984.709)*math.sin(α))*((realHomeElev-realSouthwestElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[5] 1: ",waterRateArr[5])
							#if it's taller or in the middle with real heights then there is still no flow which normal model does fine
							#(by in the middle I mean when the real height is great enough to impede flow coming in but not great enough to make any of its own flow) 
							#See the next if statement for example
						elif(elevData[i+1][k-1]<home): #if shorter
							if(realSouthwestElev>home and realSouthwestElev<realHomeElev): #but if the water reaches a point where it stops some of the water from flowing over
								#essentially if the water level of the other reaches above the ground level of the other then whatever water it 'covers' won't flow
								#recalculate
								α=math.atan((realHomeElev - realSouthwestElev)/(math.sqrt(2)*tileDim)) # all in m
								waterRateArr[5] = ((14984.709)*math.sin(α))*((realHomeElev-realSouthwestElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[5] 2: ",waterRateArr[5])
							elif(realSouthwestElev>=realHomeElev): #if taller in real elevation
								#recalculate
								waterRateArr[5]=0
								if(printExtras):
									print("waterRateArr[5] 3: ",waterRateArr[5])
							#elif(realSouthwestElev<home): if it were simply below
								#don't recalculate as the normal model does this just fine
					#check Western tile
					if(k!=0): #if exists
						realWestElev = elevData[i][k-1]+currentSurfaceWater[i][k-1]/100 # in m
						if(elevData[i][k-1]>=home): #if taller
							if(realWestElev<realHomeElev): #if shorter in real heights (heights with water added)
								#recalculate
								α=math.atan((realHomeElev - realWestElev)/tileDim) #all in m
								waterRateArr[6] = ((14984.709)*math.sin(α))*((realHomeElev-realWestElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[6] 1: ",waterRateArr[6])
							#if it's taller or in the middle with real heights then there is still no flow which normal model does fine
							#(by in the middle I mean when the real height is great enough to impede flow coming in but not great enough to make any of its own flow) 
							#See the next if statement for example
						elif(elevData[i][k-1]<home): #if north is shorter
							if(realWestElev>home and realWestElev<realHomeElev): #but if the water reaches a point where it stops some of the water from flowing over
								#essentially if the water level of the other reaches above the ground level of the other then whatever water it 'covers' won't flow
								#recalculate
								α=math.atan((realHomeElev - realWestElev)/tileDim) # all in m
								waterRateArr[6] = ((14984.709)*math.sin(α))*((realHomeElev-realWestElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[6] 2: ",waterRateArr[6])
							elif(realWestElev>=realHomeElev): #if its taller in real elevation
								#recalculate
								waterRateArr[6]=0
								if(printExtras):
									print("waterRateArr[6] 3: ",waterRateArr[6])
							#elif(realWestElev<home): if it were simply below
								#don't recalculate as the normal model does this just fine
					#check northwestern tile
					if(i!=0 and k!=0): #if exists
						realNorthwestElev = elevData[i-1][k-1]+currentSurfaceWater[i-1][k-1]/100 # in m
						if(elevData[i-1][k-1]>=home): #if taller
							if(realNorthwestElev<realHomeElev): #if shorter in real heights (heights with water added)
								#recalculate
								α=math.atan((realHomeElev - realNorthwestElev)/(math.sqrt(2)*tileDim)) #all in m
								waterRateArr[7] = ((14984.709)*math.sin(α))*((realHomeElev-realNorthwestElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[7] 1: ",waterRateArr[7])
							#if it's taller or in the middle with real heights then there is still no flow which normal model does fine
							#(by in the middle I mean when the real height is great enough to impede flow coming in but not great enough to make any of its own flow) 
							#See the next if statement for example
						elif(elevData[i-1][k-1]<home): #if shorter
							if(realNorthwestElev>home and realNorthwestElev<realHomeElev): #but if the water reaches a point where it stops some of the water from flowing over
								#essentially if the water level of the other reaches above the ground level of the other then whatever water it 'covers' won't flow
								#recalculate
								α=math.atan((realHomeElev - realNorthwestElev)/(math.sqrt(2)*tileDim)) # all in m
								waterRateArr[7] = ((14984.709)*math.sin(α))*((realHomeElev-realNorthwestElev)**3)/(tileDim)*stepSize*1000000 # into cm
								if(printExtras):
									print("waterRateArr[7] 2: ",waterRateArr[7])
							elif(realNorthwestElev>=realHomeElev): #if taller in real elevation
								#recalculate
								waterRateArr[7]=0
								if(printExtras):
									print("waterRateArr[7] 3: ",waterRateArr[7])
							#elif(realNorthwestElev<home): if it were simply below
								#don't recalculate as the normal model does this just fine
					n=0
					while(n<8):
						waterRateArr[n]=waterRateArr[n]
						n+=1
				
					n=0
					totalRate =0
					while(n<8):
						totalRate+=waterRateArr[n]
						flowsTemp[n][i][k]=waterRateArr[n]
						n+=1
					#print("Total Rate 2: ", totalRate)
					#print("water rate array 2: ", waterRateArr)
				
			
			

			n=0
			if(totalRate==0):
				n=8
			while(n<8):
				#if the system can drain more than is there then divide by capacity to drain
				if(totalRate>currentSurfaceWater[i][k]):
					percentWater = (waterRateArr[n])/totalRate

				#if the system can't drain all the water there then simply drain what can be drained
				if(totalRate<=currentSurfaceWater[i][k]):
					#percent water equals the rate divided by what is there because then when it is multiplied by what is there to get the actual water 
					#it returns the actual rate
					if(currentSurfaceWater[i][k]!=0):
						percentWater = waterRateArr[n]/currentSurfaceWater[i][k]
					else:
						percentWater = 0

				waterToGo = percentWater*currentSurfaceWater[i][k]

				homeElev = elevData[i][k]+currentSurfaceWater[i][k]/100 #real elevation

				


				if(waterToGo!=0): #if tile doesn't exist water to go is zero
					if(n==0):
						#if(waterToGo>homeElev-elevData[i-1][k]-currentSurfaceWater[i-1][k]): #if it wants to flow more than the actual difference in height just do the difference in height
							#waterToGo=homeElev-elevData[i-1][k]-currentSurfaceWater[i-1][k]
						futureSurfaceWater[i-1][k] += waterToGo
						realWaterDeficit += waterToGo
						if(printExtras):
							print("Percent Water: ",percentWater)
							print("Tile: ", i-1," ",k)
							print("Future water", futureSurfaceWater[i-1][k])
					elif(n==1):
						#if(waterToGo>homeElev-elevData[i-1][k+1]-currentSurfaceWater[i-1][k+1]): #if it wants to flow more than the actual difference in height just do the difference in height
							#waterToGo=homeElev-elevData[i-1][k+1]-currentSurfaceWater[i-1][k+1]
						futureSurfaceWater[i-1][k+1] += waterToGo
						realWaterDeficit += waterToGo
						if(printExtras):
							print("Percent Water: ",percentWater)
							print("Tile: ", i-1," ",k+1)
							print("Future water", futureSurfaceWater[i-1][k+1])
					elif(n==2):
						#if(waterToGo>homeElev-elevData[i][k+1]-currentSurfaceWater[i][k+1]): #if it wants to flow more than the actual difference in height just do the difference in height
							#waterToGo=homeElev-elevData[i][k+1]-currentSurfaceWater[i][k+1]
						futureSurfaceWater[i][k+1] += waterToGo
						realWaterDeficit += waterToGo
						if(printExtras):
							print("Percent Water: ",percentWater)
							print("Tile: ", i," ",k+1)
							print("Future water", futureSurfaceWater[i][k+1])
					elif(n==3):
						#if(waterToGo>homeElev-elevData[i+1][k+1]-currentSurfaceWater[i+1][k+1]): #if it wants to flow more than the actual difference in height just do the difference in height
							#waterToGo=homeElev-elevData[i+1][k+1]-currentSurfaceWater[i+1][k+1]
						futureSurfaceWater[i+1][k+1] += waterToGo
						realWaterDeficit += waterToGo
						if(printExtras):
							print("Percent Water: ",percentWater)
							print("Tile: ", i+1," ",k+1)
							print("Future water", futureSurfaceWater[i+1][k+1])
					elif(n==4):
						#if(waterToGo>homeElev-elevData[i+1][k]-currentSurfaceWater[i+1][k]): #if it wants to flow more than the actual difference in height just do the difference in height
							#waterToGo=homeElev-elevData[i+1][k]-currentSurfaceWater[i+1][k]
						futureSurfaceWater[i+1][k] += waterToGo
						realWaterDeficit += waterToGo
						if(printExtras):
							print("Percent Water: ",percentWater)
							print("Tile: ", i+1," ",k)
							print("Future water", futureSurfaceWater[i+1][k])
					elif(n==5):
						#if(waterToGo>homeElev-elevData[i+1][k-1]-currentSurfaceWater[i+1][k-1]): #if it wants to flow more than the actual difference in height just do the difference in height
							#waterToGo=homeElev-elevData[i+1][k-1]-currentSurfaceWater[i+1][k-1]
						futureSurfaceWater[i+1][k-1] += waterToGo
						realWaterDeficit += waterToGo
						if(printExtras):
							print("Percent Water: ",percentWater)
							print("Tile: ", i+1," ",k-1)
							print("Future water", futureSurfaceWater[i+1][k-1])
					elif(n==6):
						#if(waterToGo>homeElev-elevData[i][k-1]-currentSurfaceWater[i][k-1]): #if it wants to flow more than the actual difference in height just do the difference in height
							#waterToGo=homeElev-elevData[i][k-1]-currentSurfaceWater[i][k-1]
						futureSurfaceWater[i][k-1] += waterToGo
						realWaterDeficit += waterToGo
						if(printExtras):
							print("Percent Water: ",percentWater)
							print("Tile: ", i," ",k-1)
							print("Future water", futureSurfaceWater[i][k-1])
					elif(n==7):
						#if(waterToGo>homeElev-elevData[i-1][k-1]-currentSurfaceWater[i-1][k-1]): #if it wants to flow more than the actual difference in height just do the difference in height
							#waterToGo=homeElev-elevData[i-1][k-1]-currentSurfaceWater[i-1][k-1]
						futureSurfaceWater[i-1][k-1] += waterToGo
						realWaterDeficit += waterToGo
						if(printExtras):
							print("Percent Water: ",percentWater)
							print("Tile: ", i-1," ",k-1)
							print("Future water", futureSurfaceWater[i-1][k-1])
 
				
				
				if(printExtras):
					print("Percent Water: ",percentWater)
					print("Tile: ", i," ",k)
					print("Future water", futureSurfaceWater[i][k])
				
				n+=1

			#account for water drained

			#Watch for high flow for erosion check for high values of currentSurfaceWater leaving 
			#(but check proportionally to stepSize such that if a certain amount is flowing **per minute** 
			#it will trigger and add it to a list

			
			if(realWaterDeficit/stepSize*tileDim/currentSurfaceWater[i][k]>35*60):  #<-- 35cm/s is cause for concern and given the current surface water is water for a time unit stepSize just moving it from one side to another
				#2cm/min over a tileDim*tileDim through a d*tileDim hole make a velocity of flow
				erosionAreas[i][k] = 1

			#can't apply water deficits yet because what if tile is drained by this then a tile in later on is adjacent to it and flows into that 'empty' tile
			#need to wait to apply deficit same way need to wait to apply future water
			waterDeficitArr[i][k] = realWaterDeficit
			waterMovement[i][k] = totalRate
				
			
			k+=1
		i+=1
			
	
	#Whether data is wanted
	if(printExtras):   
		print("future surface water: ")
		print(futureSurfaceWater)
		print("water deficit:")
		print(waterDeficitArr)
	if(printResults):
		X, Y = np.meshgrid(np.linspace(0, len(elevData[0])-1, len(elevData[0])), np.linspace(0, len(elevData)-1, len(elevData)))

		Z = waterMovement/stepSize

		fig, ax = plt.subplots()

		#Build a 2D mesh with previous values of X and Y and overlay Z values as color
		print("Water Movement:")
		print("Max: ",np.amax(Z)/stepSize," cm/min")
		print("Scale: 0-10cm")
		ax.pcolormesh(X, Y, Z, vmin=0, vmax=5)
		plt.show()
	
	#Now that all water has been moved: Add future to current and set to zero 
	#also remove the water from current that has been moved
	i=0
	while(i<len(futureSurfaceWater)):
		k=0
		while(k<len(futureSurfaceWater[0])):
			#print(currentSurfaceWater[i][k])
			#print(futureSurfaceWater[i][k])
			currentSurfaceWater[i][k] -= waterDeficitArr[i][k]
			currentSurfaceWater[i][k] += futureSurfaceWater[i][k]
			#print("Current Surface water at: ",i, " ", k, " ",currentSurfaceWater[i][k])
			futureSurfaceWater[i][k] = 0
			waterDeficitArr[i][k] = 0
			k+=1
		i+=1

	if(printResults):	
		#build X and Y such that it goes from 0 to len(elevData)-1 in len(elevData) number of bits
		X, Y = np.meshgrid(np.linspace(0, len(elevData[0])-1, len(elevData[0])), np.linspace(0, len(elevData)-1, len(elevData)))

		#build z values from X and Y arrays to a corresponding array for graphing and visualizing
		Z = currentSurfaceWater

		fig, ax = plt.subplots()

		#Build a 2D mesh with previous values of X and Y and overlay Z values as color
		print("Current Surface Water:")
		print("Max: ",np.amax(Z)," cm")
		print("Scale: 0-10cm")
		ax.pcolormesh(X, Y, Z, vmin=0, vmax=10)
		plt.show()

		#reassign z values with data type float to make it cooperate with Z of type 2D float array
		Z = np.array(elevData, dtype = float)

		fig, ax = plt.subplots()

		#Build a 2D mesh with previous values of X and Y and overlay Z values as color
		print("Elevation Data:")
		print("Scale: ",np.amin(Z),"-",np.amax(Z))
		ax.pcolormesh(X, Y, Z, vmin=np.amin(Z), vmax=np.amax(Z))

		plt.show()

		#reassign z values with data type float to make it cooperate with Z of type 2D float array
		Z = np.array(absorptionData, dtype = float)

		fig, ax = plt.subplots()

		#Build a 2D mesh with previous values of X and Y and overlay Z values as color
		#print("Max: ",np.amax(Z))
		print("Absorption Capacity:")
		print("Scale: ",(np.amin(Z)),"-",np.amax(Z)," cm/min")
		ax.pcolormesh(X, Y, Z, vmin=np.amin(Z), vmax=np.amax(Z))

		plt.show()

		n=0
		while(n<8):
			Z = flowsTemp[n]
	
			fig, ax = plt.subplots()
	
			#Build a 2D mesh with previous values of X and Y and overlay Z values as color
			#print("Max: ",np.amax(Z))
			print("flows ",n,": ")
			print("Scale: ",(np.amin(Z)),"-",np.amax(Z)," cm/min")
			ax.pcolormesh(X, Y, Z, vmin=np.amin(Z), vmax=np.amax(Z))
	
			plt.show()
			n+=1
	
	t+=stepSize

	#Check for problem Areas throughout flow model (not only at end)
	i=0
	while(i<len(elevData)):
		k=0
		while(k<len(elevData[0])):
			if(currentSurfaceWater[i][k]>5): #5cm of standing water is ~2in which is enough to damage structures
				#essentially boolean value
				problemAreas[i][k] = 1
			k+=1
		i+=1


	#Update Rainfall according to careful exponential decay function: https://www.desmos.com/calculator/wm53ucf5aq
	i=0
	while(i<len(elevData)):
		k=0
		while(k<len(elevData[0])):
			initial = absorptionDataStatic[i][k]
			absorptionData[i][k] = initial/5*(4*math.exp(-t/4)+1)
			k+=1
		i+=1

	#print more
	percentDone = t/timescale

	#if running for efficiency (not printing a ton of debug measures) then print updates of progress
	if(not(printResults) and not(printExtras) and percentDone<=1):
		#os.system('cls')

		if(percentDone*100>percentGoal): #such that it only prints out 100 updates at most (just cause its annoying scrolling through 520 percent done updates)

			#construct a string with - proportional to percent finished
			result = "("
			while(percentDone > 1/40):
				result+="-"
				percentDone -= 1/20 #twenty because twenty digits in bar
	
			#fill rest with 0s 
			while(len(result)<=20):
				result+="0"
	
			result += ")"
			print(result)
			if(printFiles):
				string1 = "currentWater%d.csv" % (int(percentGoal))
				np.savetxt(string1, currentSurfaceWater, delimiter=",")

	percentDone = t/timescale
	if(percentDone*100>percentGoal):
		#just always print percent done
		if(percentDone<=1):
			end = time.time()
			print("Percent Done: ", math.trunc(percentDone*1000)/10,"% -- ", (end - start), "seconds")
		percentGoal+=1
	#string1 = "currentWater%f.csv" % (t)
	#np.savetxt(string1, currentSurfaceWater, delimiter=",")

#algorithm done --> process and print results
endF = time.time()
print("Total Elapsed Time: ", math.trunc((endF-startF)*1000)/1000)

print("Current Surface Water: ",currentSurfaceWater)

#rebuild X and Y
X, Y = np.meshgrid(np.linspace(0, len(elevData[0])-1, len(elevData[0])), np.linspace(0, len(elevData)-1, len(elevData)))

#build z values
Z = currentSurfaceWater

fig, ax = plt.subplots()

#Build a 2D mesh with previous values of X and Y and overlay Z values as color
print("Final Surface Water:")
print("Max: ",np.amax(Z), "cm")
print("Scale: 0-10 cm")
ax.pcolormesh(X, Y, Z, vmin=0, vmax=10)

plt.show()

fig, ax = plt.subplots()
print("Final Surface Water:")
print("Max: ",np.amax(Z), "cm")
print("Scale: 0-",np.amax(Z)," cm")
ax.pcolormesh(X, Y, Z, vmin=0, vmax=np.amax(Z))

plt.show()

#Print Problem Areas:
Z = np.array(problemAreas, dtype = float) 

#dtype float is essentially just to make it cooperate with being assigned to the same Z that other float arrays have been assigned to earlier in the program 

fig, ax = plt.subplots()

print("Problem Areas:")
ax.pcolormesh(X, Y, Z, vmin=0, vmax=1)

plt.show()

#Print Erosion Problem Areas:
Z = np.array(erosionAreas, dtype = float) 

fig, ax = plt.subplots()

print("Erosion Problem Areas:")
ax.pcolormesh(X, Y, Z, vmin=0, vmax=1)

plt.show()

#print real elevations
i=0
while(i<len(elevData)):
	k=0
	while(k<len(elevData[0])):
		Z[i][k]=elevData[i][k]+currentSurfaceWater[i][k]/100
		k+=1
	i+=1

fig, ax = plt.subplots()

print("Real Elevation (water included):")
print("Scale: ",np.amin(Z),"-",np.amax(Z)," m")
ax.pcolormesh(X, Y, Z, vmin=np.amin(Z), vmax=np.amax(Z))

plt.show()

#Print Elevations for Refrence:
Z = np.array(elevData, dtype = float) 

fig, ax = plt.subplots()

print("Elevation:")
print("Scale: ",np.amin(Z),"-",np.amax(Z)," m")
ax.pcolormesh(X, Y, Z, vmin=np.amin(Z), vmax=np.amax(Z))

plt.show()




#Print flows for north:
Z = np.array(flows[0], dtype = float) 

fig, ax = plt.subplots()

#the way the graphing bit works it takes the array and puts 0,0 in the bottom left whereas I visualize my model as having 0,0 in the top left
#So I say this is northward flow but the graph it shows makes it south (I will interchange North and South from now on) 
print("South Flow:")
print("Scale: ",np.amin(Z),"-",np.amax(Z)," cm/step")
ax.pcolormesh(X, Y, Z, vmin=np.amin(Z), vmax=np.amax(Z))

plt.show()

#Print flows for northeast:
Z = np.array(flows[1], dtype = float) 

fig, ax = plt.subplots()

print("Southeast Flow:")
print("Scale: ",np.amin(Z),"-",np.amax(Z)," cm/step")
ax.pcolormesh(X, Y, Z, vmin=np.amin(Z), vmax=np.amax(Z))

plt.show()

#Print flows for east:
Z = np.array(flows[2], dtype = float) 

fig, ax = plt.subplots()

print("East Flow:")
print("Scale: ",np.amin(Z),"-",np.amax(Z)," cm/step")
ax.pcolormesh(X, Y, Z, vmin=np.amin(Z), vmax=np.amax(Z))

plt.show()

#Print flows for southeast:
Z = np.array(flows[3], dtype = float) 

fig, ax = plt.subplots()

print("Northeast Flow:")
print("Scale: ",np.amin(Z),"-",np.amax(Z)," cm/step")
ax.pcolormesh(X, Y, Z, vmin=np.amin(Z), vmax=np.amax(Z))

plt.show()

#Print flows for south:
Z = np.array(flows[4], dtype = float)

fig, ax = plt.subplots()

print("North Flow:")
print("Scale: ",np.amin(Z),"-",np.amax(Z)," cm/step")
ax.pcolormesh(X, Y, Z, vmin=np.amin(Z), vmax=np.amax(Z))

plt.show()

#Print flows for southwest:
Z = np.array(flows[5], dtype = float)

fig, ax = plt.subplots()

print("Northwest Flow:")
print("Scale: ",np.amin(Z),"-",np.amax(Z)," cm/step")
ax.pcolormesh(X, Y, Z, vmin=np.amin(Z), vmax=np.amax(Z))

plt.show()

#Print flows for west:
Z = np.array(flows[6], dtype = float)

fig, ax = plt.subplots()

print("West Flow:")
print("Scale: ", np.amin(Z), "-", np.amax(Z), " cm/step")
ax.pcolormesh(X, Y, Z, vmin=np.amin(Z), vmax=np.amax(Z))

plt.show()

#Print flows for northwest:
Z = np.array(flows[7], dtype = float)

fig, ax = plt.subplots()

print("Southwest Flow:")
print("Scale: ", np.amin(Z), "-", np.amax(Z), " cm/step")
ax.pcolormesh(X, Y, Z, vmin=np.amin(Z), vmax=np.amax(Z))

plt.show()

np.savetxt("FinalWater.csv", currentSurfaceWater, delimiter=",")
#get output --> tie back to problem or how to get there from what you have if there isn't
#enough time for implimentation
#			 inspiration from Utah State University Hydrology Program and Methodology
#			 absorption differentiation