Added piston tension model.

Also typographic cleanup 'assert (' -> 'assert('.

Currently broken: stiffness calculation with the tension model.

diff --git a/src/sawsim.nw b/src/sawsim.nw
index 9917130a961264dfc716da6aa1e0afc4847d6e75..c247d7c49f48dca8b1abe9a619c55a80f2f0c546 100644 (file)
--- a/src/sawsim.nw
+++ b/src/sawsim.nw
@@ -264,6 +264,7 @@ following models have been implemented:
  \item  Hookean spring (Appendix \ref{sec.hooke}),
  \item  Worm-like chain (WLC, Appendix \ref{sec.wlc}), and
  \item  Freely-jointed chain (FJC, Appendix \ref{sec.fjc}).
+ \item  Piston (Appendix \ref{sec.piston_tension}),
 \end{packed_item}
 
 The tension model and parameters used for a particular domain depend
@@ -418,7 +419,7 @@ you accidentally assign incompatible domains to the same group you
 should get an appropriate error message.
 
 <<find tension definitions>>=
-#define NUM_TENSION_MODELS 5
+#define NUM_TENSION_MODELS 6
 
 @
 
@@ -806,7 +807,8 @@ tension_model_getopt_t tension_models[NUM_TENSION_MODELS] = {
   <<constant tension model getopt>>,
   <<hooke tension model getopt>>,
   <<worm-like chain tension model getopt>>,
-  <<freely-jointed chain tension model getopt>>
+  <<freely-jointed chain tension model getopt>>,
+  <<piston tension model getopt>>
 };
 @
 
@@ -947,7 +949,7 @@ void get_options(int argc, char **argv,
   extern char *optarg;
   extern int optind, optopt, opterr;
 
-  assert (argc > 0);
+  assert(argc > 0);
 
 #ifdef DEBUG
   for (i=0; i<n_tension_models; i++) {
@@ -2517,7 +2519,7 @@ double group_stiffness(one_dim_fn_t *tension_handler, void *handler_data, double
   F = (*tension_handler)(x, handler_data);
   Fl = (*tension_handler)(xl, handler_data);
   stiffness = (F-Fl)/dx;
-  assert(stiffness > 0);
+  assert(stiffness >= 0);
   return stiffness;
 }
 @
@@ -3695,6 +3697,7 @@ The interface is defined in [[tension_model.h]], the implementation in [[tension
 <<hooke tension model declarations>>
 <<worm-like chain tension model declarations>>
 <<freely-jointed chain tension model declarations>>
+<<piston tension model declarations>>
 <<find tension definitions>>
 <<tension model global declarations>>
 #endif /* TENSION_MODEL_H */
@@ -3764,7 +3767,7 @@ double const_tension_handler(double x, void *pdata)
   list_t *list = data->group_tension_params;
   double F;
 
-  assert (x >= 0.0);
+  assert(x >= 0.0);
   list = head(list);
   assert(list != NULL); /* empty active group?! */
   F = ((const_tension_param_t *)list->d)->F;
@@ -3829,7 +3832,7 @@ char const_tension_param_string[] = "0";
 @
 
 <<constant tension model getopt>>=
-{&const_tension_handler, NULL, "const", "an infinitely stretchable domain with constant tension", 1, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
+{&const_tension_handler, NULL, "const", "an infinitely stretchable domain with constant tension", num_const_tension_params, const_tension_param_descriptions, const_tension_param_string, &string_create_const_tension_param_t, &destroy_const_tension_param_t}
 @
 
 \subsection{Hooke}
@@ -3896,7 +3899,7 @@ double hooke_handler(double x, void *pdata)
 {
   hooke_param_t param = {0};
 
-  assert (x >= 0.0);
+  assert(x >= 0.0);
   hooke_param_grouper((tension_handler_data_t *)pdata, &param);
   return param.k*x;
 }
@@ -3909,7 +3912,7 @@ double inverse_hooke_handler(double F, void *pdata)
 {
   hooke_param_t param = {0};
 
-  assert (F >= 0.0);
+  assert(F >= 0.0);
   hooke_param_grouper((tension_handler_data_t *)pdata, &param);
   return F/param.k;
 }
@@ -4132,6 +4135,7 @@ double wlc_handler(double x, void *pdata)
   tension_handler_data_t *data = (tension_handler_data_t *)pdata;
   wlc_param_t param = {0};
 
+  assert(x >= 0.0);
   wlc_param_grouper(data, &param);
   return wlc(x, data->env->T, param.p, param.L);
 }
@@ -4214,7 +4218,7 @@ char wlc_param_string[]="0.39e-9,28.4e-9";
 @
 
 <<worm-like chain tension model getopt>>=
-{wlc_handler, inverse_wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t}
+{&wlc_handler, &inverse_wlc_handler, "wlc", "a worm-like chain (F=kbT/p{.25[(1-x/L)^-2 -1] + x/L})", num_wlc_params, wlc_param_descriptions, (char *)wlc_param_string, &string_create_wlc_param_t, &destroy_wlc_param_t}
 @
 Titin parameters from Carrion-Vazquez 1999, PNAS USA 96, 3696.
 
@@ -4348,6 +4352,145 @@ char fjc_param_string[]="0.5e-9,200";
 {fjc_handler, NULL, "fjc", "a freely-jointed chain [F=kbT/l L^-1(x/L)]", num_fjc_params, fjc_param_descriptions, (char *)fjc_param_string, &string_create_fjc_param_t, &destroy_fjc_param_t}
 @
 
+\subsection{Piston}
+\label{sec.piston_tension}
+
+An infinitely stretchable domain with no tension for extensions less
+than a particular threshold $L$, and infinite tension for $x>L$.  The
+tension at $x=L$ is undefined, but may be any positive value.  The
+model is called the ``piston'' model because it resembles a
+frictionless piston sliding in a rigid cylinder of length $L$.
+
+Because the tension at $x=L$ is undefined, the user must make sure
+domains with this tension model are never the initial active group,
+because it would break [[tension_balance()]] in [[find_tension()]]
+(see Section \ref{sec.tension_balance}).
+
+<<piston tension functions>>=
+<<piston tension parameter grouper>>
+<<piston tension handler>>
+<<inverse piston tension handler>>
+<<piston tension parameter create/destroy functions>>
+@
+
+<<piston tension model declarations>>=
+<<piston tension handler declaration>>
+<<inverse piston tension handler declaration>>
+<<piston tension parameter create/destroy function declarations>>
+<<piston tension model global declarations>>
+
+@
+
+First we define a function that computes the effective piston parameters
+(stored in a single [[piston_tension_param_t]]) for the entire group.
+<<piston tension parameter grouper>>=
+static void piston_tension_param_grouper(tension_handler_data_t *data,
+                                         piston_tension_param_t *param) {
+  list_t *list = data->group_tension_params;
+  double L=0;
+
+  list = head(list);
+  assert(list != NULL); /* empty active group?! */
+  while (list != NULL) {
+    L += ((piston_tension_param_t *)list->d)->L;
+    list = list->next;
+  }
+  param->L = L;
+}
+
+<<piston tension handler declaration>>=
+extern double piston_tension_handler(double x, void *pdata);
+@
+<<piston tension handler>>=
+double piston_tension_handler(double x, void *pdata)
+{
+  tension_handler_data_t *data = (tension_handler_data_t *)pdata;
+  piston_tension_param_t param = {0};
+  double F;
+
+  piston_tension_param_grouper(data, &param);
+  if (x <= param.L) F = 0;
+  else F = HUGE_VAL;
+#ifdef DEBUG
+  fprintf(stderr, "%s: x %g, L %g, tension %g N\n", __FUNCTION__, x, param.L, F);
+#endif
+  return F;
+}
+
+@
+
+<<inverse piston tension handler declaration>>=
+extern double inverse_piston_tension_handler(double x, void *pdata);
+@
+<<inverse piston tension handler>>=
+double inverse_piston_tension_handler(double F, void *pdata)
+{
+  tension_handler_data_t *data = (tension_handler_data_t *)pdata;
+  piston_tension_param_t param = {0};
+
+  assert(F >= 0.0);
+  piston_tension_param_grouper(data, &param);
+  if (F == 0.0) return 0;
+  return param.L;
+}
+
+@
+
+<<piston tension structure>>=
+typedef struct piston_tension_param_struct {
+  double L; /* length of piston in m */
+} piston_tension_param_t;
+
+@
+
+<<piston tension parameter create/destroy function declarations>>=
+extern void *string_create_piston_tension_param_t(char **param_strings);
+extern void destroy_piston_tension_param_t(void *p);
+
+@
+
+<<piston tension parameter create/destroy functions>>=
+piston_tension_param_t *create_piston_tension_param_t(double L)
+{
+  piston_tension_param_t *ret
+    = (piston_tension_param_t *) malloc(sizeof(piston_tension_param_t));
+  assert(ret != NULL);
+  ret->L = L;
+  return ret;
+}
+
+void *string_create_piston_tension_param_t(char **param_strings)
+{
+  return create_piston_tension_param_t(
+      safe_strtod(param_strings[0],__FUNCTION__));
+}
+
+void destroy_piston_tension_param_t(void *p)
+{
+  if (p)
+    free(p);
+}
+
+@
+
+<<piston tension model global declarations>>=
+extern int num_piston_tension_params;
+extern char *piston_tension_param_descriptions[];
+extern char piston_tension_param_string[];
+
+@
+
+<<piston tension model globals>>=
+int num_piston_tension_params = 1;
+char *piston_tension_param_descriptions[] = {"length L, m"};
+char piston_tension_param_string[] = "0";
+
+@
+
+<<piston tension model getopt>>=
+{&piston_tension_handler, &inverse_piston_tension_handler, "piston", "a domain that slides without tension for x<L, but has infinite tension for x>L.", num_piston_tension_params, piston_tension_param_descriptions, piston_tension_param_string, &string_create_piston_tension_param_t, &destroy_piston_tension_param_t}
+@
+
 \subsection{Tension model implementation}
 
 <<tension-model.c>>=
@@ -4376,6 +4519,7 @@ char fjc_param_string[]="0.5e-9,200";
 <<hooke structure>>
 <<worm-like chain structure>>
 <<freely-jointed chain structure>>
+<<piston tension structure>>
 @
 
 <<tension model globals>>=
@@ -4384,6 +4528,7 @@ char fjc_param_string[]="0.5e-9,200";
 <<hooke tension model globals>>
 <<worm-like chain tension model globals>>
 <<freely-jointed chain tension model globals>>
+<<piston tension model globals>>
 @
 
 <<tension model functions>>=
@@ -4392,9 +4537,9 @@ char fjc_param_string[]="0.5e-9,200";
 <<hooke functions>>
 <<worm-like chain functions>>
 <<freely-jointed chain functions>>
+<<piston tension functions>>
 @
 
-
 \subsection{Utilities}
 
 The tension models can get complicated, and users may want to reassure
@@ -4545,7 +4690,7 @@ void get_options(int argc, char **argv, environment_t *env,
   extern char *optarg;
   extern int optind, optopt, opterr;
 
-  assert (argc > 0);
+  assert(argc > 0);
 
   /* setup defaults */
 
@@ -5894,7 +6039,7 @@ void get_options(int argc, char **argv, environment_t *env,
   extern char *optarg;
   extern int optind, optopt, opterr;
 
-  assert (argc > 0);
+  assert(argc > 0);
 
   /* setup defaults */