From 1781f2737d865fcf547edd82597f038f417d6a96 Mon Sep 17 00:00:00 2001
From: "W. Trevor King" <wking@drexel.edu>
Date: Wed, 16 Jul 2008 23:50:01 +0000
Subject: [PATCH] Added Fmax, xmin, and xmax options to k_model_utils

git-svn-id: svn://abax.physics.drexel.edu/sawsim/trunk@9 865a22a6-13cc-4084-8aa6-876098d8aa20
---
 sawsim.nw | 45 +++++++++++++++++++++++++++++++--------------
 1 file changed, 31 insertions(+), 14 deletions(-)

diff --git a/sawsim.nw b/sawsim.nw
index 5a67a2c..3b78bcc 100644
--- a/sawsim.nw
+++ b/sawsim.nw
@@ -4307,10 +4307,11 @@ int main(int argc, char **argv)
   unsigned int flags;
   int num_param_args; /* for INIT_MODEL() */
   char **param_args;  /* for INIT_MODEL() */
+  double Fmax;
   double special_xmin;
   double special_xmax;
   get_options(argc, argv, &env, NUM_K_MODELS, k_models, &mode, &model,
-              &special_xmin, &special_xmax, &flags);
+              &Fmax, &special_xmin, &special_xmax, &flags);
   if (flags & VFLAG) {
     printf("#initializing model %s with parameters %s\n", model->name, model->params);
   }
@@ -4321,28 +4322,39 @@ int main(int argc, char **argv)
         printf("No k function for model %s\n", model->name);
 	exit(0);
       }
+      if (Fmax <= 0) {
+        printf("Fmax = %g <= 0\n", Fmax);
+	exit(1);
+      }
       {
-        double dF=1e-11, F=0, k=1.0;
+        double F, k=1.0;
+	int i, N=200;
         printf("#F (N)\tk (%% pop. per s)\n");
-        while (k > 0 && k < 1e5) {
-          k = (*model->k)(F, &env, params);
+        for (i=0; i<=N; i++) {
+	  F = Fmax*i/(double)N;
+	  k = (*model->k)(F, &env, params);
+	  if (k < 0) break;
           printf("%g\t%g\n", F, k);
-          F += dF;
         }
       }
       break;
     case M_SPECIAL :
       if (model->k == NULL || model->k == &bell_k) {
         printf("No special mode for model %s\n", model->name);
-	exit(0);
+	exit(1);
+      }
+      if (special_xmax <= special_xmin) {
+        printf("Xmax = %g <= %g = Xmin\n", special_xmax, special_xmin);
+	exit(1);
       }
       {
-        double dx=(special_xmax-special_xmin)/500, x=special_xmin, E;
+        double Xrng=(special_xmax-special_xmin), x, E;
+	int i, N=500;
         printf("#x (m)\tE (J)\n");
-        while (x <= special_xmax) {
+	for (i=0; i<=N; i++) {
+	  x = special_xmin + Xrng*i/(double)N;
           E = multi_model_E(model->k, params, &env, x);
           printf("%g\t%g\n", x, E);
-          x += dx;
         }
       }
       break;
@@ -4428,6 +4440,7 @@ void help(char *prog_name,
   printf("  ...\n");
   printf("-m\tChange output to standard mode\n");
   printf("-M\tChange output to special mode\n");
+  printf("-F\tSet the maximum F value for the standard mode k(F) output\n");
   printf("-x\tSet the minimum x value for the special mode E(x) output\n");
   printf("-X\tSet the maximum x value for the special mode E(x) output\n");
   printf("-V\tChange output to verbose mode\n");
@@ -4447,10 +4460,11 @@ void help(char *prog_name,
 void get_options(int argc, char **argv, environment_t *env,
 	         int n_k_models, k_model_getopt_t *k_models,
                  enum mode_t *mode, k_model_getopt_t **model,
-		 double *special_xmin, double *special_xmax, unsigned int *flags)
+		 double *Fmax, double *special_xmin, double *special_xmax,
+		 unsigned int *flags)
 {
   char *prog_name = NULL;
-  char c, options[] = "T:C:k:K:mMx:X:Vh";
+  char c, options[] = "T:C:k:K:mMF:x:X:Vh";
   int k_model=0;
   extern char *optarg;
   extern int optind, optopt, opterr;
@@ -4464,7 +4478,9 @@ void get_options(int argc, char **argv, environment_t *env,
   *mode = M_K_OF_F;
   *flags = 0;
   *model = k_models;
+  *Fmax = 1e-9;
   *special_xmax = 1e-8;
+  *special_xmin = 0.0;
 
   while ((c=getopt(argc, argv, options)) != -1) {
     switch(c) {
@@ -4477,11 +4493,12 @@ void get_options(int argc, char **argv, environment_t *env,
     case 'K':
       k_models[k_model].params = optarg;
       break;
-    case 'm': *mode = M_K_OF_F;   break;
-    case 'M': *mode = M_SPECIAL;  break;
+    case 'm': *mode = M_K_OF_F;             break;
+    case 'M': *mode = M_SPECIAL;            break;
+    case 'F': *Fmax = atof(optarg);         break;
     case 'x': *special_xmin = atof(optarg); break;
     case 'X': *special_xmax = atof(optarg); break;
-    case 'V': *flags |= VFLAG;    break;
+    case 'V': *flags |= VFLAG;              break;
     case '?':
       fprintf(stderr, "unrecognized option '%c'\n", optopt);
       /* fall through to default case */
-- 
2.26.2