Update to the interactive model, pretty sweet
From GeoMod
Here is the update. I have rewritten the interactive code. I give credit to Dr. Urbano for the slider code I clipped from his coriolis code. I edited it somewhat to show the current value of the slider position once the slider is moved. I gave up on the plot display for now. Apparently trying to annalyze 30,000 points at once really bogs down the processor and makes on the fly adjustments difficult. So, instead I opted for a readout of flow and Alkalinity for the 4 sample point reaches. This is much faster and creates a real time model with all adjustable elements.
IN the Calculation phase, I set Nitrogen, Phosphate, an pH at a fixed amount because the possible range is so small that it would have a negligable effect on the outcome. So, only Carbonate and organic carbon are variable factors for Alkalinity. Still a bug or two in the flow calculation, but I should have it soon enough.
heres the code:
from visual import *
from visual.controls import *
from visual.graph import *
from math import *
### Code by Daniel A. Neilans
### University of Memphis
### Dec. 2005
butrad = 1.0
# 4 sets for surface runoff
oc1zero = 20
oc1line = cylinder(pos=(-83,oc1zero,0), axis = (0,5,0), radius = 0.5)
oc1but = sphere(pos=(-83,oc1zero,0), radius = butrad, color=(0,1,0))
oc1stop = sphere(pos=(-86,oc1zero,0), radius = butrad, color=(1,0,0))
oc1label = label(pos=oc1line.pos+vector(0,5,0), text="Organic Carbon 1" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
flow1zero = 10
flow1line = cylinder(pos=(-73,flow1zero,0), axis = (0,5,0), radius = 0.5)
flow1but = sphere(pos=(-73,flow1zero,0), radius = butrad, color=(0,1,0))
flow1stop = sphere(pos=(-76,flow1zero,0), radius = butrad, color=(1,0,0))
flow1label = label(pos=flow1line.pos+vector(0,5,0), text="Overland Flow Rate 1" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
oc2zero = 0
oc2line = cylinder(pos=(-63,oc2zero,0), axis = (0,5,0), radius = 0.5)
oc2but = sphere(pos=(-63,oc2zero,0), radius = butrad, color=(0,1,0))
oc2stop = sphere(pos=(-66,oc2zero,0), radius = butrad, color=(1,0,0))
oc2label = label(pos=oc2line.pos+vector(0,5,0), text="Organic Carbon 2" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
flow2zero = -10
flow2line = cylinder(pos=(-53,flow2zero,0), axis = (0,5,0), radius = 0.5)
flow2but = sphere(pos=(-53,flow2zero,0), radius = butrad, color=(0,1,0))
flow2stop = sphere(pos=(-56,flow2zero,0), radius = butrad, color=(1,0,0))
flow2label = label(pos=flow2line.pos+vector(0,5,0), text="Overland Flow Rate 2" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
oc3zero = 20
oc3line = cylinder(pos=(-43,oc3zero,0), axis = (0,5,0), radius = 0.5)
oc3but = sphere(pos=(-43,oc3zero,0), radius = butrad, color=(0,1,0))
oc3stop = sphere(pos=(-46,oc3zero,0), radius = butrad, color=(1,0,0))
oc3label = label(pos=oc3line.pos+vector(0,5,0), text="Organic Carbon 3" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
flow3zero = 10
flow3line = cylinder(pos=(-33,flow3zero,0), axis = (0,5,0), radius = 0.5)
flow3but = sphere(pos=(-33,flow3zero,0), radius = butrad, color=(0,1,0))
flow3stop = sphere(pos=(-36,flow3zero,0), radius = butrad, color=(1,0,0))
flow3label = label(pos=flow3line.pos+vector(0,5,0), text="Overland Flow Rate 3" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
oc4zero = 0
oc4line = cylinder(pos=(-23,oc4zero,0), axis = (0,5,0), radius = 0.5)
oc4but = sphere(pos=(-23,oc4zero,0), radius = butrad, color=(0,1,0))
oc4stop = sphere(pos=(-26,oc4zero,0), radius = butrad, color=(1,0,0))
oc4label = label(pos=oc4line.pos+vector(0,5,0), text="Organic Carbon 4" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
flow4zero = -10
flow4line = cylinder(pos=(-13,flow4zero,0), axis = (0,5,0), radius = 0.5)
flow4but = sphere(pos=(-13,flow4zero,0), radius = butrad, color=(0,1,0))
flow4stop = sphere(pos=(-16,flow4zero,0), radius = butrad, color=(1,0,0))
flow4label = label(pos=flow4line.pos+vector(0,5,0), text="Overland Flow Rate 4" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
# 4 sets for ground water
carb1zero = 20
carb1line = cylinder(pos=(-3,carb1zero,0), axis = (0,5,0), radius = 0.5)
carb1but = sphere(pos=(-3,carb1zero,0), radius = butrad, color=(0,1,0))
carb1stop = sphere(pos=(-6,carb1zero,0), radius = butrad, color=(1,0,0))
carb1label = label(pos=carb1line.pos+vector(0,5,0), text="Carbonate Concentration 1" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
aq1zero = 10
aq1line = cylinder(pos=(3,aq1zero,0), axis = (0,5,0), radius = 0.5)
aq1but = sphere(pos=(3,aq1zero,0), radius = butrad, color=(0,1,0))
aq1stop = sphere(pos=(6,aq1zero,0), radius = butrad, color=(1,0,0))
aq1label = label(pos=aq1line.pos+vector(0,5,0), text="Aquifer Discharge 1" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
carb2zero = 0
carb2line = cylinder(pos=(13,carb2zero,0), axis = (0,5,0), radius = 0.5)
carb2but = sphere(pos=(13,carb2zero,0), radius = butrad, color=(0,1,0))
carb2stop = sphere(pos=(16,carb2zero,0), radius = butrad, color=(1,0,0))
carb2label = label(pos=carb2line.pos+vector(0,5,0), text="Carbonate Concentration 2" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
aq2zero = -10
aq2line = cylinder(pos=(23,aq2zero,0), axis = (0,5,0), radius = 0.5)
aq2but = sphere(pos=(23,aq2zero,0), radius = butrad, color=(0,1,0))
aq2stop = sphere(pos=(26,aq2zero,0), radius = butrad, color=(1,0,0))
aq2label = label(pos=aq2line.pos+vector(0,5,0), text="Aquifer Discharge 2" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
carb3zero = 20
carb3line = cylinder(pos=(33,carb3zero,0), axis = (0,5,0), radius = 0.5)
carb3but = sphere(pos=(33,carb3zero,0), radius = butrad, color=(0,1,0))
carb3stop = sphere(pos=(36,carb3zero,0), radius = butrad, color=(1,0,0))
carb3label = label(pos=carb3line.pos+vector(0,5,0), text="Carbonate Concentration 3" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
aq3zero = 10
aq3line = cylinder(pos=(43,aq3zero,0), axis = (0,5,0), radius = 0.5)
aq3but = sphere(pos=(43,aq3zero,0), radius = butrad, color=(0,1,0))
aq3stop = sphere(pos=(46,aq3zero,0), radius = butrad, color=(1,0,0))
aq3label = label(pos=aq3line.pos+vector(0,5,0), text="Aquifer Discharge 3" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
carb4zero = 0
carb4line = cylinder(pos=(53,carb4zero,0), axis = (0,5,0), radius = 0.5)
carb4but = sphere(pos=(53,carb4zero,0), radius = butrad, color=(0,1,0))
carb4stop = sphere(pos=(56,carb4zero,0), radius = butrad, color=(1,0,0))
carb4label = label(pos=carb4line.pos+vector(0,5,0), text="Carbonate Concentration 4" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
aq4zero = -10
aq4line = cylinder(pos=(63,aq4zero,0), axis = (0,5,0), radius = 0.5)
aq4but = sphere(pos=(63,aq4zero,0), radius = butrad, color=(0,1,0))
aq4stop = sphere(pos=(66,aq4zero,0), radius = butrad, color=(1,0,0))
aq4label = label(pos=aq4line.pos+vector(0,5,0), text="Aquifer Discharge 4" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
### set up display labels
alk1label = label(pos=flow2line.pos+vector(0,-15,0), text="Alkalinity 1" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
river1label = label(pos=flow2line.pos+vector(0,-25,0), text="River Flow 1" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
alk2label = label(pos=flow4line.pos+vector(0,-15,0), text="Alkalinity 2" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
river2label = label(pos=flow4line.pos+vector(0,-25,0), text="River Flow 2" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
alk3label = label(pos=aq2line.pos+vector(0,-15,0), text="Alkalinity 3" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
river3label = label(pos=aq2line.pos+vector(0,-25,0), text="River Flow 3" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
alk4label = label(pos=aq4line.pos+vector(0,-15,0), text="Alkalinity 4" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
river4label = label(pos=aq4line.pos+vector(0,-25,0), text="River Flow 4" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
alk1val = label(pos=flow2line.pos+vector(0,-35,0), text="Alkalinity 1" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
river1val = label(pos=flow2line.pos+vector(0,-45,0), text="River Flow 1" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
alk2val = label(pos=flow4line.pos+vector(0,-35,0), text="Alkalinity 2" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
river2val = label(pos=flow4line.pos+vector(0,-45,0), text="River Flow 2" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
alk3val = label(pos=aq2line.pos+vector(0,-35,0), text="Alkalinity 3" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
river3val = label(pos=aq2line.pos+vector(0,-45,0), text="River Flow 3" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
alk4val = label(pos=aq4line.pos+vector(0,-35,0), text="Alkalinity 4" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
river4val = label(pos=aq4line.pos+vector(0,-45,0), text="River Flow 4" ,height=10,border=2, xoffset=6, yoffset=5, space = 2)
### set the defaults for sliders
pick = None
framerate = 360.0
new_pos= vector(0,0,0)
oc1 = 0
flow1 = 0
oc2 = 0
flow2 = 0
oc3 = 0
flow3 = 0
oc4 = 0
flow4 = 0
carb1 = 0
aq1 = 0
carb2 = 0
aq2 = 0
carb3 = 0
aq3 = 0
carb4 = 0
aq4 = 0
alkval = 0
alkval = 0
alkval = 0
alkval = 0
river1val = 0
river2val = 0
river3val = 0
river4val = 0
##info for graphing
pos1 = array([35.311139,-89.639073],Float)#hwy 70
pos2 = array([35.306869,-89.668077],Float)# Hwy 205
pos3 = array([35.281040,-89.766085],Float)# New Brunswick
pos4 = array([35.275161,-89.887326],Float)#204 Singleton Pkwy
pos5 = array([35.263172,-89.929469],Float)#Raliegh Millington
pos6 = array([35.259577,-89.993882],Float)#Hwy 51
pos7 = array([35.254232,-90.026146],Float)#N. Watkins
wellpos1 = array((35.244580,-89.956463))
wellpos2 = array((35.037847,-89.855829))
wellpos3 = array((35.020714,-90.121397))
#flow for pos location in months(Jan,Feb,mar...Dec)
#flow1 = array([498,673,659,574,399,274,203,163,171,179,349,633],Float)
#flow2 = 0
#flow3 = array([1368,1591,1308,1108,699,358,262,198,228,141,444,831],Float)
#flow4 = 0
#flow5 = array([2.03,1.19,2.99,2.97,3.32,0.99,0.92,0.38,0.55,0.24,1.43,2.83],Float)
#flow6 = array([0,0,0,0,0,0,353,213,164,0,0,0],Float)
d1_2 = int(sqrt(pow((pos2[0]-pos1[0]),2)+pow((pos2[1]-pos1[1]),2))*91021.855)
d2_3 = int(sqrt(pow((pos3[0]-pos2[0]),2)+pow((pos3[1]-pos2[1]),2))*91021.855)
d3_4 = int(sqrt(pow((pos4[0]-pos3[0]),2)+pow((pos4[1]-pos3[1]),2))*91021.855)
d4_5 = int(sqrt(pow((pos5[0]-pos4[0]),2)+pow((pos5[1]-pos4[1]),2))*91021.855)
d5_6 = int(sqrt(pow((pos6[0]-pos5[0]),2)+pow((pos6[1]-pos5[1]),2))*91021.855)
d6_7 = int(sqrt(pow((pos7[0]-pos6[0]),2)+pow((pos7[1]-pos6[1]),2))*91021.855)
# deg * 91021.855 = meters
p1 = (0)
p2 = d1_2
p3 = p2 + d2_3
p4 = p3 + d3_4
p5 = p4 + d4_5
p6 = p5 + d5_6
p7 = p6 + d6_7
scale = 1000
p1 = int(p1/scale)
p2 = int(p2/scale)
p3 = int(p3/scale)
p4 = int(p3/scale)
p5 = int(p5/scale)
p6 = int(p6/scale)
p7 = int(p7/scale)
#print d1_2,d2_3,d3_4,d4_5,d5_6,d6_7
sec1 = arrayrange(0, d1_2, 1)
sec2 = arrayrange(d1_2, d2_3, 1)
sec3 = arrayrange(d2_3, d3_4, 1)
sec4 = arrayrange(d3_4, d4_5, 1)
sec5 = arrayrange(d4_5, d5_6, 1)
sec6 = arrayrange(d5_6, d6_7, 1)
plotvals = ones((1,(p7)),Float)
flow_vol = ones((1,(p7)),Float) * 150
#print plotvals
#print plotpos
#print flow_vol
#flow_profile = display(xtitle='Distance(m)', ytitle='flow rate(ft3/s)', title='Flow Profile Plot' )
#functa = gdots(color=color.blue)
#functb = gcurve(color=color.red)
while 1:
rate(framerate)
if scene.mouse.events:
m1 = scene.mouse.getevent() # obtain drag or drop event
if m1.drag and (m1.pick == oc1but or m1.pick == flow1but or m1.pick == oc2but or m1.pick == flow2but
or m1.pick == oc3but or m1.pick == flow3but or m1.pick == oc4but or m1.pick == flow4but
or m1.pick == carb1but or m1.pick == aq1but or m1.pick == carb2but or m1.pick == aq2but
or m1.pick == carb3but or m1.pick == aq3but or m1.pick == carb4but or m1.pick == aq4but):
drag_pos = m1.pickpos
pick = m1.pick
scene.cursor.visible = 0 # make cursor invisible
elif m1.drop:
pick = None # end dragging
scene.cursor.visible = 1 # cursor visible
if pick:
new_pos = scene.mouse.project(normal=(0,0,1))
if new_pos != drag_pos:
pick.pos += new_pos - drag_pos
drag_pos = new_pos
#print new_pos.y
# initialize overland controls
if pick == oc1but:
oc1label.text = "%0.2f" % new_pos.y
pick.pos.x = oc1line.pos.x
pick.pos.z = oc1line.pos.z
if new_pos.y > 5+oc1zero:
pick.pos.y = 5+oc1zero
elif new_pos.y < oc1zero:
pick.pos.y = oc1zero
oc1 = new_pos.y
if pick == flow1but:
flow1label.text = "%0.2f" % new_pos.y
pick.pos.x = flow1line.pos.x
pick.pos.z = flow1line.pos.z
if new_pos.y > 5+flow1zero:
pick.pos.y = 5+flow1zero
elif new_pos.y < flow1zero:
pick.pos.y = flow1zero
flow1 = new_pos.y
if pick == oc2but:
oc2label.text = "%0.2f" % new_pos.y
pick.pos.x = oc2line.pos.x
pick.pos.z = oc2line.pos.z
if new_pos.y > 5+oc2zero:
pick.pos.y = 5+oc2zero
elif new_pos.y < oc2zero:
pick.pos.y = oc2zero
oc2 = new_pos.y
if pick == flow2but:
flow2label.text = "%0.2f" % new_pos.y
pick.pos.x = flow2line.pos.x
pick.pos.z = flow2line.pos.z
if new_pos.y > 5+flow2zero:
pick.pos.y = 5+flow2zero
elif new_pos.y < flow2zero:
pick.pos.y = flow2zero
flow2 = new_pos.y
if pick == oc3but:
oc3label.text = "%0.2f" % new_pos.y
pick.pos.x = oc3line.pos.x
pick.pos.z = oc3line.pos.z
if new_pos.y > 5+oc3zero:
pick.pos.y = 5+oc3zero
elif new_pos.y < oc3zero:
pick.pos.y = oc3zero
oc3 = new_pos.y
if pick == flow3but:
flow3label.text = "%0.2f" % new_pos.y
pick.pos.x = flow3line.pos.x
pick.pos.z = flow3line.pos.z
if new_pos.y > 5+flow3zero:
pick.pos.y = 5+flow3zero
elif new_pos.y < flow3zero:
pick.pos.y = flow3zero
flow3 = new_pos.y
if pick == oc4but:
oc4label.text = "%0.2f" % new_pos.y
pick.pos.x = oc4line.pos.x
pick.pos.z = oc4line.pos.z
if new_pos.y > 5+oc4zero:
pick.pos.y = 5+oc4zero
elif new_pos.y < oc4zero:
pick.pos.y = oc4zero
oc4 = new_pos.y
if pick == flow4but:
flow4label.text = "%0.2f" % new_pos.y
pick.pos.x = flow4line.pos.x
pick.pos.z = flow4line.pos.z
if new_pos.y > 5+flow4zero:
pick.pos.y = 5+flow4zero
elif new_pos.y < flow4zero:
pick.pos.y = flow4zero
flow4 = new_pos.y
#initialize aquifer controls
if pick == carb1but:
carb1label.text = "%0.2f" % new_pos.y
pick.pos.x = carb1line.pos.x
pick.pos.z = carb1line.pos.z
if new_pos.y > 5+carb1zero:
pick.pos.y = 5+carb1zero
elif new_pos.y < carb1zero:
pick.pos.y = carb1zero
carb1 = new_pos.y
if pick == aq1but:
aq1label.text = "%0.2f" % new_pos.y
pick.pos.x = aq1line.pos.x
pick.pos.z = aq1line.pos.z
if new_pos.y > 5+aq1zero:
pick.pos.y = 5+aq1zero
elif new_pos.y < aq1zero:
pick.pos.y = aq1zero
aq1 = new_pos.y
if pick == carb2but:
carb2label.text = "%0.2f" % new_pos.y
pick.pos.x = carb2line.pos.x
pick.pos.z = carb2line.pos.z
if new_pos.y > 5+carb2zero:
pick.pos.y = 5+carb2zero
elif new_pos.y < carb2zero:
pick.pos.y = carb2zero
carb2 = new_pos.y
if pick == aq2but:
aq2label.text = "%0.2f" % new_pos.y
pick.pos.x = aq2line.pos.x
pick.pos.z = aq2line.pos.z
if new_pos.y > 5+aq2zero:
pick.pos.y = 5+aq2zero
elif new_pos.y < aq2zero:
pick.pos.y = aq2zero
aq2 = new_pos.y
if pick == carb3but:
carb3label.text = "%0.2f" % new_pos.y
pick.pos.x = carb3line.pos.x
pick.pos.z = carb3line.pos.z
if new_pos.y > 5+carb3zero:
pick.pos.y = 5+carb3zero
elif new_pos.y < carb3zero:
pick.pos.y = carb3zero
carb3 = new_pos.y
if pick == aq3but:
aq3label.text = "%0.2f" % new_pos.y
pick.pos.x = aq3line.pos.x
pick.pos.z = aq3line.pos.z
if new_pos.y > 5+aq3zero:
pick.pos.y = 5+aq3zero
elif new_pos.y < aq3zero:
pick.pos.y = aq3zero
aq3 = new_pos.y
if pick == carb4but:
carb4label.text = "%0.2f" % new_pos.y
pick.pos.x = carb4line.pos.x
pick.pos.z = carb4line.pos.z
if new_pos.y > 5+carb4zero:
pick.pos.y = 5+carb4zero
elif new_pos.y < carb4zero:
pick.pos.y = carb4zero
carb4 = new_pos.y
if pick == aq4but:
aq4label.text = "%0.2f" % new_pos.y
pick.pos.x = aq4line.pos.x
pick.pos.z = aq4line.pos.z
if new_pos.y > 5+aq4zero:
pick.pos.y = 5+aq4zero
elif new_pos.y < aq4zero:
pick.pos.y = aq4zero
aq4 = new_pos.y
### analyze all points for flow and Alkalinity
for i in range((p1+1),p2):
plotvals[0,i] = flow1*(oc1 - .15 - .50 - .85) + aq1*(carb1 - .80) + flow_vol[0,i]*plotvals[0,i]
flow_vol[0,i] = flow_vol[0,i-1] + flow1 + aq1
for i in range(p2+1,p3):
plotvals[0,i] = flow2*(oc2 - .15 - .50 - .85) + aq2*(carb2 - .80) + flow_vol[0,i]*plotvals[0,i]
flow_vol[0,i] = flow_vol[0,i-1] + flow2 + aq2
for i in range(p3+1,p4):
plotvals[0,i] = flow3*(oc3 - .15 - .50 - .85) + aq3*(carb3 - .80) + flow_vol[0,i]*plotvals[0,i]
flow_vol[0,i] = flow_vol[0,i-1] + flow3 + aq3
for i in range(p4+1,p7):
plotvals[0,i] = flow4*(oc4 - .15 - .50 - .85) + aq4*(carb4 - .80) + flow_vol[0,i]*plotvals[0,i]
flow_vol[0,i] = flow_vol[0,i-1] + flow4 + aq4
print plotvals[0,p1], plotvals[0,p1+10]
### display the calculated values
alk1val.text = "%0.2f" % plotvals[0,p2-1]
alk2val.text = "%0.2f" % plotvals[0,p3-1]
alk3val.text = "%0.2f" % plotvals[0,p4-1]
alk4val.text = "%0.2f" % plotvals[0,p6-1]
plotvals = ones((1,(p7)),Float)
flow_vol = ones((1,(p7)),Float) * 150
'''
river1val.text = "%0.3f" % flow_vol[0,p2]
river2val.text = "%0.3f" % flow_vol[0,p3]
river3val.text = "%0.3f" % flow_vol[0,p4]
river4val.text = "%0.3f" % flow_vol[0,p6]
'''
#for i in range(p1, p7):
#lin = curve(color=(0,1,0))
#lin.append(pos=(i,flow_vol[0,i],0))
#for i in range(p1,p7):
#functa.plot(pos=(i,flow_vol[0,i]))
#functb.plot(pos=(i,plotvals[0,i]))
