introduce function for linear interpolation.

1 files changed, 21 insertions, 8 deletions

diff --git a/makeBoundary b/makeBoundary
index 76eaf21..590bcaf 100755
--- a/makeBoundary
+++ b/makeBoundary
@@ -112,6 +112,12 @@ def save_bc(outputfile):
                     -csa[ind_z]*Q_i[panel_x][-1]/csa_p[panel_x],
                     zmax-zitem))
 
+def interp(extra2, max1, min1, max2, min2):
+    # use similar triangles to perform linear interpolation
+    extra1 = min1 + (max1 - min1)*(extra2 - min2)/(max2 - min2)
+    
+    return extra1
+
 # read command line argument:
 parser = argparse.ArgumentParser(
     description="generate FullSWOF boundary files")
@@ -294,16 +300,23 @@ ind_q = bisect.bisect(Q_i[panel[ind_p]], panel_target_flow)
 print('insertion point =', ind_q)
 
 # find height at target flow by linear interpolation
-h_extra = h_i[panel[ind_p]][ind_q-1] \
-    + (panel_target_flow - Q_i[panel[ind_p]][ind_q-1])\
-    *(h_i[panel[ind_p]][ind_q] - h_i[panel[ind_p]][ind_q-1])\
-    /(Q_i[panel[ind_p]][ind_q] - Q_i[panel[ind_p]][ind_q-1])
+h_extra = interp(
+    panel_target_flow,
+    h_i[panel[ind_p]][ind_q],
+    h_i[panel[ind_p]][ind_q-1],
+    Q_i[panel[ind_p]][ind_q],
+    Q_i[panel[ind_p]][ind_q-1])
+
 print('heights:', h_i[panel[ind_p]][ind_q-1], h_extra, h_i[panel[ind_p]][ind_q])
+
 # find area at target flow by linear interpolation
-A_extra = A_i[panel[ind_p]][ind_q-1] \
-    + (h_extra - h_i[panel[ind_p]][ind_q-1])\
-    *(A_i[panel[ind_p]][ind_q] - A_i[panel[ind_p]][ind_q-1])\
-    /(h_i[panel[ind_p]][ind_q] - h_i[panel[ind_p]][ind_q-1])
+A_extra = interp(
+    h_extra,
+    A_i[panel[ind_p]][ind_q],
+    A_i[panel[ind_p]][ind_q-1],
+    h_i[panel[ind_p]][ind_q],
+    h_i[panel[ind_p]][ind_q-1])
+    
 print('hydraulic radii:', r_h[panel[ind_p]][ind_q-1], r_h[panel[ind_p]][ind_q])
 
 velocity_panel    = panel_target_flow/A_extra