Fix omega -> \omega in manual's timescale discussion.
[sawsim.git] / src / sawsim.nw
index c9b9a3436757e4f96c765b38dc3f755c393113b5..5f89121e6f23443520fdf9bdeabd19a3c47f85c2 100644 (file)
@@ -1,3 +1,23 @@
+% sawsim - simulating single molecule mechanical unfolding.
+% Copyright (C) 2008-2010  William Trevor King
+%
+% This program is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program.  If not, see <http://www.gnu.org/licenses/>.
+%
+% The author may be contacted at <wking@drexel.edu> on the Internet, or
+% write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
+% Philadelphia PA 19104, USA.
+
 %%%%%%%%%%%%%%%%%%%%%%%%% Start LaTeX boilerplate %%%%%%%%%%%%%%%%%%%%%%%%%%%
 % we have our own preamble,
 % use `noweave -delay` to not print noweb's automatic one
@@ -107,7 +127,29 @@ to share the burden and increase the transparency of data analysis.
 
 \subsection{Related work}
 
-TODO References
+Sawim has been published\citep{king10}!  It is, as far as I know, the
+only open-source Monte Carlo force spectroscopy simulator, but similar
+closed source simulators have been around since the seminal paper by
+\citet{rief97a}.  In chronological order, major contributions have
+come from
+
+\begin{packed_item}
+  \item \citet{rief97a} \citeyear{rief97a}:
+  \item \citet{rief97b} \citeyear{rief97b}:
+  \item \citet{kellermayer97} \citeyear{kellermayer97}:
+  \item \citet{rief98} \citeyear{rief98}:
+    first paper to focus mainly on the simulation
+  \item \citet{oberhauser98} \citeyear{oberhauser98}:
+  \item \citet{carrion-vazquez99a} \citeyear{carrion-vazquez99a}:
+  \item \citet{best02} \citeyear{best02}:
+    pointed out large fit valley
+  \item \citet{zinober02} \citeyear{zinober02}:
+    introduced the scaffold effect
+  \item \citet{jollymore09} \citeyear{jollymore09}:
+  \item \citet{king10} \citeyear{king10}:
+    introduced sawsim
+\end{packed_item}
+
 
 \subsection{About this document}
 
@@ -175,8 +217,8 @@ handle domains with undefined stiffness.
 \subsection{License}
 
 <<license>>=
- sawsim - program for simulating single molecule mechanical unfolding.
- Copyright (C) 2008-2009, William Trevor King
+ sawsim - simulating single molecule mechanical unfolding.
+ Copyright (C) 2008-2010, William Trevor King
 
  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
@@ -189,9 +231,7 @@ handle domains with undefined stiffness.
  See the GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
- along with this program; if not, write to the Free Software
- Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
- 02111-1307, USA.
+ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
  The author may be contacted at <wking@drexel.edu> on the Internet, or
  write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
@@ -645,7 +685,7 @@ along the chain.  The quality of this assumption depends on your
 particular chain.  For example, a damped spring thermalizes on a
 timescale of order $\tau = 1/\gamma$ where the damping ratio $\gamma
 \equiv \omega_0(-\zeta \pm \sqrt{\zeta^2-1}$, the natural angular
-frequency $omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
+frequency $\omega_0 \equiv \sqrt{k/m}$, $\zeta \equiv b/2\sqrt{km}$,
 and $b$ sets the drag $F_b = -bv$.  For our cantilevers $\tau$ is
 on the order of milliseconds, which is longer than a timestep.
 However, the approximation is still reasonable, because there is
@@ -719,11 +759,12 @@ allows interchangeable models, and we are currently focusing on the
 simulation itself, so we remove the actual model implementation to the
 appendices.
 
-The tension models are defined in Appendix \ref{sec.tension_models}
+The tension models are defined in Appendix \ref{sec.tension_models},
 and unfolding models are defined in Appendix \ref{sec.k_models}.
 
 <<model globals>>=
 #define k_B   1.3806503e-23  /* J/K */
+
 @
 
 
@@ -1241,6 +1282,7 @@ We include [[math.h]], so don't forget to link to the libm with `-lm'.
 #include "tension_model.h"
 #include "parse.h"
 #include "accel_k.h"
+
 @
 
 <<definitions>>=
@@ -1325,6 +1367,7 @@ We also define a macro for our [[check]] unit testing
                 -(received-expected)/expected, max_err, #received,           \
                  expected, received);                                        \
   } while(0)
+
 @
 
 <<happens>>=
@@ -2031,6 +2074,7 @@ void free_parsed_list(int num, char **string_array)
 
 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
 <<check-parse.c>>=
+<<license comment>>
 <<parse check includes>>
 <<string comparison definition and globals>>
 <<check relative error>>
@@ -2055,10 +2099,10 @@ Here we check to make sure the various functions work as expected, using \citeta
 <<parse suite definition>>=
 Suite *test_suite (void)
 {
-  Suite *s = suite_create ("k model");
+  Suite *s = suite_create ("parse");
   <<parse list string test case defs>>
 
-  <<parse list string test case add>>
+  <<parse list string test case adds>>
   return s;
 }
 @
@@ -2132,7 +2176,7 @@ END_TEST
 <<parse list string test case defs>>=
 TCase *tc_parse_list_string = tcase_create("parse list string");
 @
-<<parse list string test case add>>=
+<<parse list string test case adds>>=
 //tcase_add_test(tc_parse_list_string, test_next_delim_index);
 tcase_add_test(tc_parse_list_string, test_parse_list_null);
 tcase_add_test(tc_parse_list_string, test_parse_list_single_simple);
@@ -2166,7 +2210,6 @@ Here we check to make sure the various functions work as expected, using \citeta
 @
 
 <<test suite>>=
-<<F tests>>
 <<determine dt tests>>
 <<happens tests>>
 <<does domain transition tests>>
@@ -2178,17 +2221,15 @@ Here we check to make sure the various functions work as expected, using \citeta
 Suite *test_suite (void)
 {
   Suite *s = suite_create ("sawsim");
-  <<F test case defs>>
   <<determine dt test case defs>>
   <<happens test case defs>>
   <<does domain transition test case defs>>
   <<randomly unfold domains test case defs>>
 
-  <<F test case add>>
-  <<determine dt test case add>>
-  <<happens test case add>>
-  <<does domain transition test case add>>
-  <<randomly unfold domains test case add>>
+  <<determine dt test case adds>>
+  <<happens test case adds>>
+  <<does domain transition test case adds>>
+  <<randomly unfold domains test case adds>>
 
 /*
   tcase_add_checked_fixture(tc_strip_address,
@@ -2212,54 +2253,6 @@ int main(void)
 }
 @
 
-\subsection{F tests}
-<<F tests>>=
-<<worm-like chain tests>>
-@
-<<F test case defs>>=
-<<worm-like chain test case def>>
-@
-<<F test case add>>=
-<<worm-like chain test case add>>
-@
-
-<<worm-like chain tests>>=
-<<worm-like chain function>>
-
-START_TEST(test_wlc_at_zero)
-{
-  double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
-  fail_unless(abs(wlc(x, T, p, L)) < lim, \
-              "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
-              x, T, p, L, abs(wlc(x, T, p, L)), lim);
-}
-END_TEST
-
-START_TEST(test_wlc_at_half)
-{
-  double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
-  /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
-   * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
-   *                = 0.25 * 3 = 3/4
-   * linear term = x/L = 0.5
-   * nonlinear + linear = 0.75 + 0.5 = 1.25
-   * wlc = 10e21*1.25 = 12.5e21
-   */
-  fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
-              "wlc(%g, %g, %g, %g) = %g != %g",
-               x, T, p, L, wlc(x, T, p, L), 12.5e21);
-}
-END_TEST
-@
-<<worm-like chain test case def>>=
-TCase *tc_wlc = tcase_create("WLC");
-@
-
-<<worm-like chain test case add>>=
-tcase_add_test(tc_wlc, test_wlc_at_zero);
-tcase_add_test(tc_wlc, test_wlc_at_half);
-suite_add_tcase(s, tc_wlc);
-@
 
 \subsection{Model tests}
 
@@ -2300,7 +2293,7 @@ END_TEST
 <<determine dt test case defs>>=
 TCase *tc_determine_dt = tcase_create("Determine dt");
 @
-<<determine dt test case add>>=
+<<determine dt test case adds>>=
 tcase_add_test(tc_determine_dt, test_determine_dt_linear_k);
 suite_add_tcase(s, tc_determine_dt);
 @
@@ -2309,21 +2302,21 @@ suite_add_tcase(s, tc_determine_dt);
 @
 <<happens test case defs>>=
 @
-<<happens test case add>>=
+<<happens test case adds>>=
 @
 
 <<does domain transition tests>>=
 @
 <<does domain transition test case defs>>=
 @
-<<does domain transition test case add>>=
+<<does domain transition test case adds>>=
 @
 
 <<randomly unfold domains tests>>=
 @
 <<randomly unfold domains test case defs>>=
 @
-<<randomly unfold domains test case add>>=
+<<randomly unfold domains test case adds>>=
 @
 
 
@@ -2447,9 +2440,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
@@ -2790,6 +2785,7 @@ double length_scale = 1e-10; /* HACK */
 
 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
 <<check-tension-balance.c>>=
+<<license comment>>
 <<tension balance check includes>>
 <<tension balance test suite>>
 <<main check program>>
@@ -3296,6 +3292,7 @@ The user must free the data pointed to by the node on their own.
 
 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
 <<check-list.c>>=
+<<license comment>>
 <<list check includes>>
 <<list test suite>>
 <<main check program>>
@@ -3593,6 +3590,7 @@ The [[while]] loop around [[wait]] protects [[wait]] from interrupting signals.
 
 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
 <<check-string-eval.c>>=
+<<license comment>>
 <<string eval check includes>>
 <<string comparison definition and globals>>
 <<string eval test suite>>
@@ -3620,7 +3618,7 @@ Suite *test_suite (void)
   Suite *s = suite_create ("string eval");
   <<string eval test case defs>>
 
-  <<string eval test case add>>
+  <<string eval test case adds>>
   return s;
 }
 @
@@ -3638,7 +3636,7 @@ START_TEST(test_invalid_command)
   FILE *in, *out;
   char input[80], output[80]={};
   string_eval_setup(&in, &out);
-  sprintf(input, "print ABCDefg\n");
+  sprintf(input, "print undefined_name__traceback_should_be_printed_to_stderr\n");
   string_eval(in, out, input, 80, output);
   string_eval_teardown(&in, &out);
 }
@@ -3685,7 +3683,7 @@ END_TEST
 <<string eval test case defs>>=
 TCase *tc_string_eval = tcase_create("string_eval");
 @
-<<string eval test case add>>=
+<<string eval test case adds>>=
 tcase_add_test(tc_string_eval, test_setup_teardown);
 tcase_add_test_raise_signal(tc_string_eval, test_invalid_command, SIGKILL);
 tcase_add_test(tc_string_eval, test_simple_eval);
@@ -3708,6 +3706,7 @@ The number of evaluation calls could be drastically reduced, however, by impleme
 <<accel-k.h>>=
 #ifndef ACCEL_K_H
 #define ACCEL_K_H
+<<license comment>>
 double accel_k(k_func_t *k, double F, environment_t *env, void *params);
 void free_accels();
 #endif /* ACCEL_K_H */
@@ -3715,11 +3714,12 @@ void free_accels();
 
 <<accel k module makefile lines>>=
 NW_SAWSIM_MODS += accel_k
-#CHECK_BINS += check_accel_k
-check_accel_k_MODS =
+CHECK_BINS += check_accel_k
+check_accel_k_MODS = interp k_model parse string_eval tavl
 @
 
 <<accel-k.c>>=
+<<license comment>>
 #include <assert.h>  /* assert()                */
 #include <stdlib.h>  /* realloc(), free(), NULL */
 #include "global.h"  /* environment_t           */
@@ -3795,6 +3795,77 @@ double accel_k(k_func_t *k, double F, environment_t *env, void *params)
 }
 @
 
+\subsection{Unit tests}
+
+Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
+<<check-accel-k.c>>=
+<<license comment>>
+<<accel k check includes>>
+<<check relative error>>
+<<model globals>>
+<<accel k test suite>>
+<<main check program>>
+@
+
+<<accel k check includes>>=
+#include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
+#include <stdio.h>  /* sprintf()                  */
+#include <assert.h> /* assert()                   */
+#include <math.h>   /* exp()                      */
+#include <errno.h>  /* errno, ERANGE              */
+<<check includes>>
+#include "global.h"
+#include "k_model.h"
+#include "accel_k.h"
+
+@
+
+<<accel k test suite>>=
+<<accel k suite tests>>
+<<accel k suite definition>>
+@
+
+<<accel k suite definition>>=
+Suite *test_suite (void)
+{
+  Suite *s = suite_create ("accelerated k model");
+  <<accel k test case defs>>
+
+  <<accel k test case adds>>
+  return s;
+}
+@
+
+<<accel k test case defs>>=
+TCase *tc_accel_k = tcase_create("accelerated k model");
+@
+
+<<accel k test case adds>>=
+tcase_add_test(tc_accel_k, test_accel_k_accuracy);
+suite_add_tcase(s, tc_accel_k);
+@
+
+<<accel k suite tests>>=
+START_TEST(test_accel_k_accuracy)
+{
+  double k, k_accel, F, Fm=0.0, FM=300e-12, dF=0.5e-12;
+  k_func_t *k_func = bell_k;
+  char *param_strings[] = {"3.3e-4", "0.25e-9"};
+  void *params = string_create_bell_param_t(param_strings);
+  environment_t env = {};
+  env.T = 300.0;
+
+  for(F=Fm; F <= FM; F+=dF) {
+    k = k_func(F, &env, params);
+    k_accel = accel_k(k_func, F, &env, params);
+    CHECK_ERR(1e-8, k, k_accel);
+  }
+  destroy_bell_param_t(params);
+}
+END_TEST
+
+@
+
 \section{Tension models}
 \label{sec.tension_models}
 
@@ -3821,8 +3892,8 @@ The interface is defined in [[tension_model.h]], the implementation in [[tension
 
 <<tension model module makefile lines>>=
 NW_SAWSIM_MODS += tension_model
-#CHECK_BINS += check_tension_model
-check_tension_model_MODS =
+CHECK_BINS += check_tension_model
+check_tension_model_MODS = list tension_balance
 @
 <<tension model utils makefile lines>>=
 TENSION_MODEL_MODS = tension_model parse list tension_balance
@@ -3832,7 +3903,7 @@ TENSION_MODEL_SRC = $(BUILD_DIR)/tension_model_utils.c $(BUILD_DIR)/global.h \
 $(BUILD_DIR)/tension_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD_DIR)
        notangle -Rtension-model-utils.c $< > $@
 $(BIN_DIR)/tension_model_utils : $(TENSION_MODEL_SRC) | $(BIN_DIR)
-       gcc -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+       $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(TENSION_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
 clean_tension_model_utils :
        rm -f $(BUILD_DIR)/tension_model_utils.c $(BIN_DIR)/tension_model_utils
 @ The pipe symbol [[|]] marks the prerequisits following it (in this
@@ -4361,14 +4432,13 @@ $$
 $$
 where $L=Nl$ is the total length of the chain, and $\mathcal{L}(\alpha) \equiv \coth{\alpha} - \frac{1}{\alpha}$ is the Langevin function\citep{hatfield99}.
 <<freely-jointed chain function>>=
-double langevin(double x, void *params)
+static double langevin(double x, void *params)
 {
   if (x==0) return 0;
   return 1.0/tanh(x) - 1/x;
 }
-<<one dimensional bracket declaration>>
-<<one dimensional solve declaration>>
-double inv_langevin(double x)
+
+static double inv_langevin(double x)
 {
   double lb=0.0, ub=1.0, ret;
   oneD_bracket(&langevin, NULL, x, 1.0, &lb, &ub);
@@ -4376,13 +4446,14 @@ double inv_langevin(double x)
   return ret;
 }
 
-double fjc(double x, double T, double l, int N)
+static double fjc(double x, double T, double l, int N)
 {
   double L = l*(double)N;
   if (x == 0) return 0;
   //printf("x %g\tL%g\tx/L %g\n", x, L, x/L);
   return k_B*T/l * inv_langevin(x/L);
 }
+
 @
 
 <<freely-jointed chain handler declaration>>=
@@ -4547,6 +4618,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);
 @
@@ -4855,6 +4947,239 @@ void get_options(int argc, char **argv, environment_t *env,
 }
 @
 
+\subsection{Tension model unit tests}
+
+Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
+<<check-tension-model.c>>=
+<<license comment>>
+<<tension model check includes>>
+<<check relative error>>
+<<model globals>>
+<<tension model test suite>>
+<<main check program>>
+@
+
+<<tension model check includes>>=
+#include <stdlib.h> /* EXIT_SUCCESS, EXIT_FAILURE */
+#include <stdio.h>  /* sprintf()                  */
+#include <assert.h> /* assert()                   */
+#include <math.h>   /* exp()                      */
+#include <errno.h>  /* errno, ERANGE              */
+<<check includes>>
+#include "global.h"
+#include "list.h"
+#include "tension_balance.h" /* oneD_*() */
+#include "tension_model.h"
+@
+
+<<tension model test suite>>=
+<<safe strto*>>
+<<const tension tests>>
+<<hooke tests>>
+<<worm-like chain tests>>
+<<freely-jointed chain tests>>
+<<piston tests>>
+<<tension model suite definition>>
+@
+
+<<tension model suite definition>>=
+Suite *test_suite (void)
+{
+  Suite *s = suite_create ("tension model");
+  <<const tension test case defs>>
+  <<hooke test case defs>>
+  <<worm-like chain test case defs>>
+  <<freely-jointed chain test case defs>>
+  <<piston test case defs>>
+
+  <<const tension test case adds>>
+  <<hooke test case adds>>
+  <<worm-like chain test case adds>>
+  <<freely-jointed chain test case adds>>
+  <<piston test case adds>>
+  return s;
+}
+@
+
+\subsubsection{Constant}
+
+<<const tension test case defs>>=
+TCase *tc_const_tension = tcase_create("Constant tension");
+@
+
+<<const tension test case adds>>=
+tcase_add_test(tc_const_tension, test_const_tension_create_destroy);
+suite_add_tcase(s, tc_const_tension);
+@
+
+<<const tension tests>>=
+START_TEST(test_const_tension_create_destroy)
+{
+  char *tension[] = {"1","2.2", "3"};
+  char *params[] = {tension[0]};
+  void *p = NULL;
+  int i;
+  for( i=0; i < sizeof(tension)/sizeof(char *); i++) {
+    params[0] = tension[i];
+    p = string_create_const_tension_param_t(params);
+    destroy_const_tension_param_t(p);
+  }
+}
+END_TEST
+
+@ TODO: In order to test the constant tension handler itself, we'd
+have to construct a group.
+
+
+\subsubsection{Hooke}
+
+<<hooke test case defs>>=
+TCase *tc_hooke = tcase_create("Hooke");
+@
+
+<<hooke test case adds>>=
+tcase_add_test(tc_hooke, test_hooke_create_destroy);
+suite_add_tcase(s, tc_hooke);
+
+@
+
+<<hooke tests>>=
+START_TEST(test_hooke_create_destroy)
+{
+  char *k[] = {"1","2.2", "3"};
+  char *params[] = {k[0]};
+  void *p = NULL;
+  int i;
+  for( i=0; i < sizeof(k)/sizeof(char *); i++) {
+    params[0] = k[i];
+    p = string_create_hooke_param_t(params);
+    destroy_hooke_param_t(p);
+  }
+}
+END_TEST
+
+@ TODO: In order to test the Hooke tension handler itself, we'd
+have to construct a group.
+
+
+\subsubsection{Worm-like chain}
+
+<<worm-like chain test case defs>>=
+TCase *tc_wlc = tcase_create("WLC");
+@
+
+<<worm-like chain test case adds>>=
+tcase_add_test(tc_wlc, test_wlc_at_zero);
+tcase_add_test(tc_wlc, test_wlc_at_half);
+suite_add_tcase(s, tc_wlc);
+
+@
+
+<<worm-like chain tests>>=
+<<worm-like chain function>>
+START_TEST(test_wlc_at_zero)
+{
+  double T=1.0, L=1.0, p=0.1, x=0.0, lim=1e-30;
+  fail_unless(abs(wlc(x, T, p, L)) < lim, \
+              "abs(wlc(x=%g,T=%g,p=%g,L=%g)) = %g >= %g",
+              x, T, p, L, abs(wlc(x, T, p, L)), lim);
+}
+END_TEST
+
+START_TEST(test_wlc_at_half)
+{
+  double T=1.0, L=1.0, p=0.1*k_B, x=0.5;
+  /* prefactor = k_B T / p = k_B 1.0 / (k_B*0.1) = 10.0 J/nm = 10.0e21 pN
+   * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
+   *                = 0.25 * 3 = 3/4
+   * linear term = x/L = 0.5
+   * nonlinear + linear = 0.75 + 0.5 = 1.25
+   * wlc = 10e21*1.25 = 12.5e21
+   */
+  fail_unless(wlc(x, T, p, L)-12.5e21 < 1e16,
+              "wlc(%g, %g, %g, %g) = %g != %g",
+               x, T, p, L, wlc(x, T, p, L), 12.5e21);
+}
+END_TEST
+
+@
+
+
+\subsubsection{Freely-jointed chain}
+
+<<freely-jointed chain test case defs>>=
+TCase *tc_fjc = tcase_create("FJC");
+@
+
+<<freely-jointed chain test case adds>>=
+tcase_add_test(tc_fjc, test_fjc_at_zero);
+tcase_add_test(tc_fjc, test_fjc_at_half);
+suite_add_tcase(s, tc_fjc);
+
+@
+
+<<freely-jointed chain tests>>=
+<<freely-jointed chain function>>
+START_TEST(test_fjc_at_zero)
+{
+  int N=10;
+  double T=1.0, l=1.0, p=0.1, x=0.0, lim=1e-30;
+  fail_unless(abs(fjc(x, T, l, N)) < lim, \
+              "abs(fjc(x=%g,T=%g,l=%g,N=%d)) = %g >= %g",
+              x, T, l, N, abs(fjc(x, T, l, N)), lim);
+}
+END_TEST
+
+START_TEST(test_fjc_at_half)
+{
+  int N = 10;
+  double T=1.0/k_B, l=1.0, x=5.0;
+  /* prefactor = k_B T / l = k_B (1.0/k_B) / 1.0 = 1.0 J/nm = 1.0e21 pN
+   * nonlinear term = 0.25 (1/(1-x/L)^2-1) = 0.25(1/.25 - 1)
+   *                = 0.25 * 3 = 3/4
+   * linear term = x/L = 0.5
+   * nonlinear + linear = 0.75 + 0.5 = 1.25
+   * fjc = 10e21*1.25 = 12.5e21
+   */
+  fail_unless(fjc(x, T, l, N)-12.5e21 < 1e16,
+              "fjc(%g, %g, %g, %d) = %g != %g",
+               x, T, l, N, fjc(x, T, l, N), 12.5e21);
+}
+END_TEST
+
+@
+
+
+\subsubsection{Piston}
+
+<<piston test case defs>>=
+TCase *tc_piston = tcase_create("Piston");
+@
+
+<<piston test case adds>>=
+tcase_add_test(tc_piston, test_piston_create_destroy);
+suite_add_tcase(s, tc_piston);
+
+@
+
+<<piston tests>>=
+START_TEST(test_piston_create_destroy)
+{
+  char *L[] = {"1","2.2", "3"};
+  char *params[] = {L[0]};
+  void *p = NULL;
+  int i;
+  for( i=0; i < sizeof(L)/sizeof(char *); i++) {
+    params[0] = L[i];
+    p = string_create_piston_tension_param_t(params);
+    destroy_piston_tension_param_t(p);
+  }
+}
+END_TEST
+
+@ TODO: In order to test the piston tension handler itself, we'd
+have to construct a group.
+
 
 \section{Unfolding rate models}
 \label{sec.k_models}
@@ -4908,7 +5233,7 @@ K_MODEL_SRC = $(BUILD_DIR)/k_model_utils.c $(BUILD_DIR)/global.h \
 $(BUILD_DIR)/k_model_utils.c : $(SRC_DIR)/sawsim.nw | $(BUILD)
        notangle -Rk-model-utils.c $< > $@
 $(BIN_DIR)/k_model_utils : $(K_MODEL_SRC) | $(BIN_DIR)
-       gcc -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+       $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(K_MODEL_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
 clean_k_model_utils :
        rm -f $(BUILD_DIR)/k_model_utils.c $(BIN_DIR)/k_model_utils
 @
@@ -5107,7 +5432,7 @@ char bell_param_string[]="3.3e-4,0.25e-9";
 <<bell k model getopt>>=
 {"bell", "thermalized, two-state transitions", &bell_k, num_bell_params, bell_param_descriptions, (char *)bell_param_string, &string_create_bell_param_t, &destroy_bell_param_t}
 @
-Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
+Initialized with Titin I27 parameters\citep{carrion-vazquez99a}.
 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
 
 
@@ -5208,7 +5533,7 @@ char kbell_param_string[]="3.3e-4,0.25e-9";
 <<kbell k model getopt>>=
 {"kbell", "thermalized, two-state transitions, corrected for effective linker stiffness", &kbell_k, num_kbell_params, kbell_param_descriptions, (char *)kbell_param_string, &string_create_kbell_param_t, &destroy_kbell_param_t}
 @
-Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
+Initialized with Titin I27 parameters\citep{carrion-vazquez99a}.
 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
 
 
@@ -5724,7 +6049,7 @@ char kramers_integ_param_string[]="0.2e-12,{(2e-9,0),(2.1e-9,-3.7e-20),(2.2e-9,-
 <<kramers integ k model getopt>>=
 {"kramers_integ", "thermalized, diffusion-limited transitions", &kramers_integ_k, num_kramers_integ_params, kramers_integ_param_descriptions, (char *)kramers_integ_param_string, &string_create_kramers_integ_param_t, &destroy_kramers_integ_param_t}
 @
-Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
+Initialized with Titin I27 parameters\citep{carrion-vazquez99a}.
 % Titin I27 parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696
 
 \subsection{Unfolding model implementation}
@@ -5783,6 +6108,7 @@ Initialized with Titin I27 parameters\citep{carrion-vazquez99b}.
 
 Here we check to make sure the various functions work as expected, using \citetalias{sw:check}.
 <<check-k-model.c>>=
+<<license comment>>
 <<k model check includes>>
 <<check relative error>>
 <<model globals>>
@@ -5952,10 +6278,10 @@ END_TEST
 <<kramers k tests>>=
 @
 
-<<kramers k test case def>>=
+<<kramers k test case defs>>=
 @
 
-<<kramers k test case add>>=
+<<kramers k test case adds>>=
 @
 
 <<k model function tests>>=
@@ -6128,7 +6454,7 @@ void help(char *prog_name,
          k_models[k_model].params);
   printf("There are two output modes.  In standard mode, k(F) is printed\n");
   printf("For example:\n");
-  printf("  #Force (N)\tk (% pop. per s)\n");
+  printf("  #Force (N)\tk (%% pop. per s)\n");
   printf("  123.456\t7.89\n");
   printf("  ...\n");
   printf("In special mode, the output depends on the model.\n");
@@ -6227,6 +6553,8 @@ $(DOC_DIR)/sawsim.pdf : $(BUILD_DIR)/sawsim.tex $(BUILD_DIR)/sawsim.bib \
        $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
        $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
        $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
+       $(SHELL) -e -c "cd $(BUILD_DIR) && bibtex sawsim"
+       $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
        $(SHELL) -e -c "cd $(BUILD_DIR) && pdflatex sawsim"
        mv $(BUILD_DIR)/sawsim.pdf $@
 @
@@ -6263,6 +6591,11 @@ Extract the makefile with
 The sed is needed because notangle replaces tabs with eight spaces,
 but [[make]] needs tabs.
 <<makefile>>=
+# Customize compilation
+CC = gcc
+CFLAGS =
+LDFLAGS =
+
 # Decide what directories to put things in
 SRC_DIR = src
 BUILD_DIR = build
@@ -6273,6 +6606,7 @@ DOC_DIR = doc
 
 # Modules (X.c, X.h) defined in the noweb file
 NW_SAWSIM_MODS =
+# Check binaries (check_*) defined in the noweb file
 CHECK_BINS =
 # Modules defined outside the noweb file
 FREE_SAWSIM_MODS = interp tavl
@@ -6298,7 +6632,11 @@ BINS = sawsim tension_model_utils k_model_utils sawsim_profile
 DOCS = sawsim.pdf
 
 # Define the major targets
-all : ./Makefile $(BINS:%=$(BIN_DIR)/%) $(DOCS:%=$(DOC_DIR)/%) ;
+all : ./Makefile all_bin all_doc
+
+.PHONY: all_bin all_doc
+all_bin : $(BINS:%=$(BIN_DIR)/%)
+all_doc : $(DOCS:%=$(DOC_DIR)/%)
 
 view : $(DOC_DIR)/sawsim.pdf
        xpdf $< &
@@ -6316,15 +6654,15 @@ clean : $(CHECK_BINS:%=clean_%) $(SAWSIM_MODS:%=clean_%) \
                $(BUILD_DIR)/interp.c $(BUILD_DIR)/interp.h \
                $(BUILD_DIR)/tavl.c $(BUILD_DIR)/tavl.h \
                $(BUILD_DIR)/global.h ./gmon.out
-       $(SHELL) -e -c "rmdir $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR)"
+       $(SHELL) -e -c "rm -rf $(BUILD_DIR) $(DOC_DIR)"
 
 # Various builds of sawsim
 $(BIN_DIR)/sawsim : $(SAWSIM_SRC) | $(BIN_DIR)
-       gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+       $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
 $(BIN_DIR)/sawsim_static : $(SAWSIM_SRC) | $(BIN_DIR)
-       gcc -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+       $(CC) $(CFLAGS) $(LDFLAGS) -g -static -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
 $(BIN_DIR)/sawsim_profile : $(SAWSIM_SRC) | $(BIN_DIR)
-       gcc -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
+       $(CC) $(CFLAGS) $(LDFLAGS) -g -pg -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(LIBS:%=-l%)
 
 # Create the directories
 $(BUILD_DIR) $(BIN_DIR) $(DOC_DIR) :
@@ -6365,7 +6703,7 @@ $(BUILD_DIR)/check_sawsim.c : $(SRC_DIR)/sawsim.nw
 $(BIN_DIR)/check_sawsim : $(BUILD_DIR)/check_sawsim.c $(BUILD_DIR)/global.h \
                $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) \
                $(SAWSIM_MODS:%=$(BUILD_DIR)/%.h) | $(BIN_DIR)
-       gcc -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
+       $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(SAWSIM_MODS:%=$(BUILD_DIR)/%.c) $(CHECK_LIBS:%=-l%)
 clean_check_sawsim :
        rm -f $(BIN_DIR)/check_sawsim $(BUILD_DIR)/check_sawsim.c
 # ... and the modules
@@ -6376,10 +6714,10 @@ $(CHECK_BINS:%=$(BIN_DIR)/%) : $$(subst $(BIN_DIR),$(BUILD_DIR),$$@).c \
                $$(patsubst %,$(BUILD_DIR)/%.c,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
                $$(patsubst %,$(BUILD_DIR)/%.h,$$($$(subst $(BIN_DIR)/,,$$@_MODS)))\
                $$(BUILD_DIR)/global.h | $(BIN_DIR)
-       gcc -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
+       $(CC) $(CFLAGS) $(LDFLAGS) -g -o $@ $< $(@:$(BIN_DIR)/check_%=$(BUILD_DIR)/%.c) \
                $(patsubst %,$(BUILD_DIR)/%.c,$($(subst $(BIN_DIR)/,,$@_MODS)))\
                $(CHECK_LIBS:%=-l%)
-# todo: clean up the required modules to
+# TODO: clean up the required modules too
 $(CHECK_BINS:%=clean_%) :
        rm -f $(@:clean_%=$(BIN_DIR)/%) $(@:clean_%=$(BUILD_DIR)/%.c)