move functions to top of file.

1 files changed, 89 insertions, 86 deletions

diff --git a/makeBoundary.py b/makeBoundary.py
index 8dc7ad0..639160f 100755
--- a/makeBoundary.py
+++ b/makeBoundary.py
@@ -16,6 +16,95 @@ def read_definition(filename):
                 ddict[items[0]] = ast.literal_eval(items[1])
     return ddict
 
+def conveyance(numH, n_co, xregion, zregion, zmin, zmax):
+    p_i = [] # wetted perimeter
+    A_i = [] # area
+    r_h = [] # hydraulic radius
+    h_i = [] # list of heights
+    K_i = [] # conveyance
+    Q_i = [] # discharge
+    x_sub = [[] for i in range(numH)] # list of x values in subregion
+    z_sub = [[] for i in range(numH)] # list of z values in subregion
+    for i in range(numH):
+        h_i.append(zmin + (i+1)*(zmax-zmin)/numH)
+        #print(zregion[zregion < h_i[i]])
+        booleanArray = zregion < h_i[i]
+        #print(booleanArray[i])
+        x_sub[i] += list(xregion[booleanArray])
+        z_sub[i] += list(zregion[booleanArray])
+        for interval in range(len(xregion)-1):
+            if booleanArray[interval+1] != booleanArray[interval]:
+                x_extra = xregion[interval] + (h_i[i]-zregion[interval])*(xregion[interval+1]-xregion[interval])/(zregion[interval+1]-zregion[interval])
+                bisect.insort(x_sub[i], x_extra) # add intercept value
+                ind_x = x_sub[i].index(x_extra)
+                z_sub[i].insert(ind_x, h_i[i])   # add height value
+        #print(z_sub[i])
+
+        dp = 0
+        dA = 0
+        eps = 1e-06
+        for j in range(len(x_sub[i])-1):
+            if abs(z_sub[i][j+1] - h_i[i]) > eps or abs(z_sub[i][j] - h_i[i]) > eps:
+                dp += np.hypot(x_sub[i][j+1] - x_sub[i][j],
+                               abs(z_sub[i][j+1] - z_sub[i][j]))
+            #print(dp)
+            # calculate area using trapezium rule
+            dA += (h_i[i] - (z_sub[i][j+1] + z_sub[i][j])/2)*(x_sub[i][j+1] - x_sub[i][j])
+            #print('Area =', dA)
+
+        p_i.append(dp)
+        A_i.append(dA)
+
+        r_h.append(A_i[i]/p_i[i])  # ratio of area and wetted perimeter
+        #print('hydraulic radius =', r_h[i])
+        K_i.append(A_i[i]*(1/n_co)*r_h[i]**(2/3)) # conveyance
+        Q_i.append(K_i[i]*slope**0.5)             # discharge
+
+    return p_i, A_i, r_h, h_i, K_i, Q_i
+
+def plot_region(xdata, labelx,
+                ydata1, labely1,
+                ydata2, labely2,
+                ydata3, labely3, titlep):
+    plt.xlabel(labelx)
+    plt.ylabel(labely1)
+    plt.title(titlep)
+    plt.plot(xdata, ydata1)
+    plt.show()
+
+    plt.xlabel(labelx)
+    plt.ylabel(labely2)
+    plt.title(titlep)
+    plt.plot(xdata, ydata2)
+    plt.show()
+
+    plt.xlabel(labelx)
+    plt.ylabel(labely3)
+    plt.title(titlep)
+    plt.plot(xdata, ydata3)
+    plt.show()
+
+def save_bc():
+    with open(outputFilename, 'w') as f:
+        f.write('{:6} {:>2} {:>2} {:>10}\n'.format('#x', 'c', 'q', 'h'))
+        for ind_z, (xitem, zitem) in enumerate(zip(xin, zin)):
+            panel_x = bisect.bisect(markers, xitem)
+            if csa[ind_z] == 0:
+                # wall boundary condition
+                f.write('{:6.2f} {:>2}\n'.format(xitem, 2))
+            elif panel_x == ind_p:
+                # imposed discharge within part-filled panel
+                f.write('{:6.2f} {:>2} {:10.6f} {:9.6f}\n'.format(
+                    xitem, 5,
+                    -csa[ind_z]*panel_target_flow/csa_p[panel_x],
+                    h_extra-zitem))
+            else:
+                # imposed discharge within filled panels
+                f.write('{:6.2f} {:>2} {:10.6f} {:9.6f}\n'.format(
+                    xitem, 5,
+                    -csa[ind_z]*Q_i[panel_x][-1]/csa_p[panel_x],
+                    zmax-zitem))
+
 # read boundary definition file.
 definition_dict = read_definition('boundaryDefinition.txt')
 #for dd in definition_dict:
@@ -126,75 +215,9 @@ print(zmin[1])
 # channel overtopping height (minimum of left bank and right bank heights):
 zmax = min(zregion[panel[0]][0], zregion[panel[0]][-1]) - ztol
 
-def conveyance(numH, n_co, xregion, zregion, zmin, zmax):
-    p_i = [] # wetted perimeter
-    A_i = [] # area
-    r_h = [] # hydraulic radius
-    h_i = [] # list of heights
-    K_i = [] # conveyance
-    Q_i = [] # discharge
-    x_sub = [[] for i in range(numH)] # list of x values in subregion
-    z_sub = [[] for i in range(numH)] # list of z values in subregion
-    for i in range(numH):
-        h_i.append(zmin + (i+1)*(zmax-zmin)/numH)
-        #print(zregion[zregion < h_i[i]])
-        booleanArray = zregion < h_i[i]
-        #print(booleanArray[i])
-        x_sub[i] += list(xregion[booleanArray])
-        z_sub[i] += list(zregion[booleanArray])
-        for interval in range(len(xregion)-1):
-            if booleanArray[interval+1] != booleanArray[interval]:
-                x_extra = xregion[interval] + (h_i[i]-zregion[interval])*(xregion[interval+1]-xregion[interval])/(zregion[interval+1]-zregion[interval])
-                bisect.insort(x_sub[i], x_extra) # add intercept value
-                ind_x = x_sub[i].index(x_extra)
-                z_sub[i].insert(ind_x, h_i[i])   # add height value
-        #print(z_sub[i])
-
-        dp = 0
-        dA = 0
-        eps = 1e-06
-        for j in range(len(x_sub[i])-1):
-            if abs(z_sub[i][j+1] - h_i[i]) > eps or abs(z_sub[i][j] - h_i[i]) > eps:
-                dp += np.hypot(x_sub[i][j+1] - x_sub[i][j],
-                               abs(z_sub[i][j+1] - z_sub[i][j]))
-            #print(dp)
-            # calculate area using trapezium rule
-            dA += (h_i[i] - (z_sub[i][j+1] + z_sub[i][j])/2)*(x_sub[i][j+1] - x_sub[i][j])
-            #print('Area =', dA)
-
-        p_i.append(dp)
-        A_i.append(dA)
-
-        r_h.append(A_i[i]/p_i[i])  # ratio of area and wetted perimeter
-        #print('hydraulic radius =', r_h[i])
-        K_i.append(A_i[i]*(1/n_co)*r_h[i]**(2/3)) # conveyance
-        Q_i.append(K_i[i]*slope**0.5)             # discharge
-
-    return p_i, A_i, r_h, h_i, K_i, Q_i
 
 #print(h_i)
 
-def plot_region(xdata, labelx,
-                ydata1, labely1,
-                ydata2, labely2,
-                ydata3, labely3, titlep):
-    plt.xlabel(labelx)
-    plt.ylabel(labely1)
-    plt.title(titlep)
-    plt.plot(xdata, ydata1)
-    plt.show()
-
-    plt.xlabel(labelx)
-    plt.ylabel(labely2)
-    plt.title(titlep)
-    plt.plot(xdata, ydata2)
-    plt.show()
-
-    plt.xlabel(labelx)
-    plt.ylabel(labely3)
-    plt.title(titlep)
-    plt.plot(xdata, ydata3)
-    plt.show()
 
 p_i = [[] for _ in range(num_panels)]
 A_i = [[] for _ in range(num_panels)]
@@ -281,25 +304,5 @@ for i, p in enumerate(panel):
 #print('csa_p[0] = {} csa_p[1] = {}'.format(csa_p[0], csa_p[1]))
 #print('A_i_west = {} A_i = {} A_i_east = {}'.format(A_extra, A_i[-1], A_i_east[-1]))
 
-def save_bc():
-    with open(outputFilename, 'w') as f:
-        f.write('{:6} {:>2} {:>2} {:>10}\n'.format('#x', 'c', 'q', 'h'))
-        for ind_z, (xitem, zitem) in enumerate(zip(xin, zin)):
-            panel_x = bisect.bisect(markers, xitem)
-            if csa[ind_z] == 0:
-                # wall boundary condition
-                f.write('{:6.2f} {:>2}\n'.format(xitem, 2))
-            elif panel_x == ind_p:
-                # imposed discharge within part-filled panel
-                f.write('{:6.2f} {:>2} {:10.6f} {:9.6f}\n'.format(
-                    xitem, 5,
-                    -csa[ind_z]*panel_target_flow/csa_p[panel_x],
-                    h_extra-zitem))
-            else:
-                # imposed discharge within filled panels
-                f.write('{:6.2f} {:>2} {:10.6f} {:9.6f}\n'.format(
-                    xitem, 5,
-                    -csa[ind_z]*Q_i[panel_x][-1]/csa_p[panel_x],
-                    zmax-zitem))
 
 save_bc()