Added inverse piston comments + whitespace in tension_balance()

diff --git a/src/sawsim.nw b/src/sawsim.nw
index c9b9a3436757e4f96c765b38dc3f755c393113b5..4f2e5c37008cd79dab3fba8b31ea53a6ad40d227 100644 (file)
--- a/src/sawsim.nw
+++ b/src/sawsim.nw
@@ -2447,9 +2447,11 @@ double tension_balance(int num_tension_groups,
         lb = x_xo_data.xi[0]- (x-last_x); /* apply all change to x0 */
         ub = x_xo_data.xi[0];
       } else { /* x == last_x */
+#ifdef DEBUG
         fprintf(stderr, "not moving\n");
-        lb= x_xo_data.xi[0];
-        ub= x_xo_data.xi[0];
+#endif
+        lb = x_xo_data.xi[0];
+        ub = x_xo_data.xi[0];
       }
     }
 #ifdef DEBUG
@@ -4547,6 +4549,27 @@ double piston_tension_handler(double x, void *pdata)
 
 @
 
+We abuse [[x_of_xo()]] and [[oneD_solve()]] with this inverse
+definition (see Section \ref{sec.tension_balance}), because the
+$x(F=0)<=L$, but our function returns $x(F=0)=0$.  This is fine when
+the total extension \emph{is} zero, but cheats a bit when there is
+some total extension.  For any non-zero extension, the initial active
+group produces some tension (we hope.  True for all our other tension
+models).  This tension fully extends the piston group (of which there
+should be only one, since all piston states can get lumped into the
+same tension group.).  If the total extension is $\le$ the target
+extension, the full extension is accurate, so the inverse was valid
+after all.  If not, [[oneD_solve()]] will halve the extension of the
+first group, reducing the overall tension.  For total extension less
+than the piston extension, this bisection continues for [[max_steps]],
+leaving $x_0$ reduced by a factor of $2^\text{max\_steps}$, a very
+small number.  Because of this, the returned tension $F_0(x_0)$ is
+very small, even though the total extension found by [[x_of_xo()]]
+is still $>$ the piston length.
+
+While this works, it is clearly not the most elegant or robust
+solution.  TODO: think of (and implement) a better procedure.
+
 <<inverse piston tension handler declaration>>=
 extern double inverse_piston_tension_handler(double x, void *pdata);
 @