Initial

.gitignore

+*~

flood.py

+from sympy import *
+
+# def wrt_w():
+#     A_sym = Symbol("A")
+#     l = Function("l")(A_sym)
+
+#     w_sym = w = Symbol("w", positive=True)
+#     x = Symbol("x")
+    
+#     a = Symbol("a", positive=True)
+#     large = 100000
+#     h = Piecewise((large*(a-x), x<-a), (0, x<=a), (large*(x-a), True))
+#     w_of_A = a
+    
+#     dh_dx = diff(h, x)
+#     dl_dx = sqrt(dh_dx**2 + 1)
+#     l = integrate(dl_dx, (x, -w/2, w/2))
+    
+#     A = w*h.subs(x, w/2) - integrate(dh_dx, (x, -w/2, w/2))
+    
+#     if w_of_A is None:
+#         w_of_A = solve(A_sym - A, w)[0]
+
+#     w = w_of_A.subs("A", A)
+
+#     print(simplify(l))
+#     print(simplify(l.subs(w_sym, w)))
+
+def main():
+    # Define symbols and functions
+    A_sym = Symbol("A", positive=True, real=True)
+    l = Function("l", positive=True, real=True)(A_sym)
+
+    h_sym = h = Symbol("h", positive=True, real=True)
+    y = Symbol("y", positive=True, real=True)
+    
+    # Width of channel wrt vertical coordinate (y)
+    a = Symbol("a", positive=True, real=True)
+    w = sqrt(y/a)
+
+    # Calculate length by integrating segment lengths between 0 and h.
+    # Add the width at the base to the total length for the case of a flat bottom.
+    dx_dy = diff(w/2, y) # x = +- w/2
+    dl_dy = 2 * sqrt(dx_dy**2 + 1) # 2* for both left and right side of channel
+    l = integrate(dl_dy, (y, 0, h)) + w.subs(y, 0)
+    print("l(h) =", l)
+
+    # Calculate channel area by integrating width wrt vertical coordinate.
+    A = simplify(integrate(w, (y, 0, h)))
+    print("A(h) =", A)
+
+    # Get h in terms of A (there may be multiple solutions)
+    h_of_A_arr = solve(A_sym - A, h)
+
+    # For each solution...
+    for h_of_A in h_of_A_arr:
+        print("================")
+        print("h(A) =")
+        pprint(h_of_A)
+
+        # Substitute the expression for h(A) into l(h) to get l(A)
+        l_of_A = simplify(l.subs(h, h_of_A))
+        print("l(A) =")
+        pprint(l_of_A)
+    
+main()