Partially done shock calculation

flood.py

@@ -3,32 +3,6 @@ from matplotlib import pyplot as plt
 from scipy import integrate as sp_int
 import numpy as np
 
-# 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 plot_expr(f, x_sym, lim, N):
     F = lambdify(x_sym, f)
     x = []
@@ -86,6 +60,36 @@ def derive_wetted_perim(w, y, A_sym, debug=False):
 
         return l_of_A
 
+def shock_step(y, t, vs, vc, dvc_dc, dc_dx0):
+    x = y[0]
+    c_L = y[1]
+    c_R = y[2]
+    
+    dx_dt = vs(c_L, c_R)
+
+    vc_L = vc(c_L)
+    vc_R = vc(c_R)
+    x0L = x - t*vc_L
+    x0R = x - t*vc_R
+    dcL_dx0 = dc_dx0(x0L)
+    dcR_dx0 = dc_dx0(x0R)
+
+    dcL_dt = dcL_dx0 * (vc_L - dx_dt) / (1 - t*dvc_dc(c_L)*dcL_dx0)
+    dcR_dt = dcR_dx0 * (vc_R - dx_dt) / (1 - t*dvc_dc(c_R)*dcR_dx0)
+
+    return [dx_dt, dcL_dt, dcR_dt]
+    
+def solve_shock(Q, c_init, x0, t0):
+    pass 
+    
+def solve_characteristics(Q, c_init, delta_x, delta_t, T):
+    ch = [{"pts":[],"x":c} for c in c_init]
+
+    for i in range(int(T/delta_t)):
+        t = delta_t * i
+        for char in ch:
+            char["pts"].append()
+    
 def godunov_solve(Q, c_init, c_inflow, delta, t):
     gstep = lambda c, t: godunov_step(Q, delta, c, c_inflow)
     return sp_int.odeint(gstep, c_init, t)