From c6c03cede40563fe7dbb83ef0dd0e49b0400d74e Mon Sep 17 00:00:00 2001
From: Joe Pater <02joepater06@gmail.com>
Date: Tue, 25 Feb 2025 11:02:27 +0000
Subject: [PATCH] Initial
---
.gitignore | 1 +
flood.py | 66 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 files changed, 67 insertions(+)
create mode 100644 .gitignore
create mode 100644 flood.py
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..b25c15b
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+*~
diff --git a/flood.py b/flood.py
new file mode 100644
index 0000000..b3ea20d
--- /dev/null
+++ b/flood.py
@@ -0,0 +1,66 @@
+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()
--
GitLab