thrown code into outdoor functions to minimize code in main functions

and futher implement window scanning

diff --git a/src/spirals.cpp b/src/spirals.cpp
index 89b06aaac26eb9a4cc19162f8bb84b45f3ca3aeb..a69a0bb5533efedb4f316230b8e292fecc92987c 100644 (file)
--- a/src/spirals.cpp
+++ b/src/spirals.cpp
@@ -198,6 +198,49 @@ RVec kernel_pro (double x0, double y0, double z0, double p1, double p2, double p
     return pro;
 }
 
+void make_kernel (std::vector< RVec > temp, double epsi, std::vector< kernel_maxima > &kernel) {
+    RVec mid, arrow;
+
+    mid[0] = 0;
+    mid[1] = 0;
+    mid[2] = 0;
+
+    arrow[0] = 0;
+    arrow[1] = 0;
+    arrow[2] = 0;
+
+    for (int i = 0; i < temp.size(); i++) {
+        rvec_inc(temp[i], mid);
+    }
+    mid[0] /= temp.size();
+    mid[1] /= temp.size();
+    mid[2] /= temp.size();
+    rvec_sub(temp.back(), temp.front(), arrow);
+
+    long double t1, t2, t3, t4, t5, t6;
+    t1 = mid[0];
+    t2 = mid[1];
+    t3 = mid[2];
+    t4 = arrow[0];
+    t5 = arrow[1];
+    t6 = arrow[2];
+
+    linear_kernel_search(t1, t2, t3, t4, t5, t6, temp, epsi);
+
+    mid[0] = t1;
+    mid[1] = t2;
+    mid[2] = t3;
+    arrow[0] = t4;
+    arrow[1] = t5;
+    arrow[2] = t6;
+
+    kernel.back().x = mid;
+    kernel.back().p = arrow;
+    for (int i = 0; i < temp.size(); i++) {
+        kernel.back().krnl.push_back(kernel_pro(mid[0], mid[1], mid[2], arrow[0], arrow[1], arrow[2], temp[i]));
+    }
+}
+
 double left_right_turn (RVec a, RVec b, RVec c) {
     return  a[0] * b[1] * c[2] +
             a[1] * b[2] * c[0] +
@@ -207,6 +250,35 @@ double left_right_turn (RVec a, RVec b, RVec c) {
             a[1] * b[0] * c[2];
 }
 
+void make_circles (std::vector< std::vector< std::vector< int > > > &circles, std::vector< RVec > temp, std::vector< kernel_maxima > kernel)  {
+    bool st1 = true, st2 = false;
+    double turn = -1, tempt;
+    RVec a, b, c;
+    rvec_sub(temp[0], kernel.back().krnl.front(), a);
+    rvec_sub(kernel.back().krnl.front(), kernel.back().krnl.back(), b);
+    for (int i = 1; i < temp.size(); i++) {
+        rvec_sub(temp[i], kernel.back().krnl[i], c);
+        tempt = left_right_turn(a, b, c);
+        if (tempt > turn) {
+            st1 = true;
+        } else {
+            st1 = false;
+            st2 = true;
+        }
+        turn = tempt;
+        if (st1 && !st2 || !st1 && st2) {
+            if (circles.back().size() == 0) {
+                circles.back().resize(1);
+            }
+            circles.back().back().push_back(i);
+        } else {
+            circles.back().resize(circles.back().size() + 1);
+            circles.back().back().push_back(i);
+            st2 = false;
+        }
+    }
+}
+
 /*! \brief
  * Template class to serve as a basis for user analysis tools.
  */
@@ -241,6 +313,7 @@ class Spirals : public TrajectoryAnalysisModule
         SelectionList                           sel_;
         int                                     frames      = 0;
         double                                  epsi        = 0.00001;
+        int                                     window      = 6;
 
         std::vector< std::vector< RVec > >                  monomers;
         std::vector< kernel_maxima >                        kernel;
@@ -265,9 +338,13 @@ Spirals::initOptions(IOptionsContainer          *options,
     // Add the descriptive text (program help text) to the options
     settings->setHelpText(desc);
     // Add option for output file name
-    options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
+    /*options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
                             .store(&fnNdx_).defaultBasename("rcore")
-                            .description("Index file from the rcore"));
+                            .description("Index file from the rcore"));*/
+    // Add option for tau constant
+    options->addOption(gmx::IntegerOption("window_length")
+                            .store(&window)
+                            .description("window length for local parameters"));
     // Add option for selection list
     options->addOption(SelectionOption("select").storeVector(&sel_)
                            .required().dynamicMask().multiValue()
@@ -300,76 +377,12 @@ Spirals::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
         monomers.back().push_back(temp[i]);
     }
 
-    RVec mid, arrow;
-
-    mid[0] = 0;
-    mid[1] = 0;
-    mid[2] = 0;
-
-    arrow[0] = 0;
-    arrow[1] = 0;
-    arrow[2] = 0;
-
-    for (int i = 0; i < sel.size(); i++) {
-        rvec_inc(temp[i], mid);
-    }
-    mid[0] /= sel.size();
-    mid[1] /= sel.size();
-    mid[2] /= sel.size();
-    rvec_sub(temp.back(), temp.front(), arrow);
-
-    long double t1, t2, t3, t4, t5, t6;
-    t1 = mid[0];
-    t2 = mid[1];
-    t3 = mid[2];
-    t4 = arrow[0];
-    t5 = arrow[1];
-    t6 = arrow[2];
-
-    linear_kernel_search(t1, t2, t3, t4, t5, t6, temp, epsi);
-
-    mid[0] = t1;
-    mid[1] = t2;
-    mid[2] = t3;
-    arrow[0] = t4;
-    arrow[1] = t5;
-    arrow[2] = t6;
-
     kernel.resize(kernel.size() + 1);
-    kernel.back().x = mid;
-    kernel.back().p = arrow;
-    for (int i = 0; i < sel.size(); i++) {
-        kernel.back().krnl.push_back(kernel_pro(mid[0], mid[1], mid[2], arrow[0], arrow[1], arrow[2], temp[i]));
-    }
+    make_kernel(temp, epsi, kernel);
 
     circles.resize(circles.size() + 1);
+    make_circles(circles, temp, kernel);
 
-    bool st1 = true, st2 = false;
-    double turn = -1, tempt;
-    RVec a, b, c;
-    rvec_sub(temp[0], kernel.back().krnl.front(), a);
-    rvec_sub(kernel.back().krnl.front(), kernel.back().krnl.back(), b);
-    for (int i = 1; i < sel.size(); i++) {
-        rvec_sub(temp[i], kernel.back().krnl[i], c);
-        tempt = left_right_turn(a, b, c);
-        if (tempt > turn) {
-            st1 = true;
-        } else {
-            st1 = false;
-            st2 = true;
-        }
-        turn = tempt;
-        if (st1 && !st2 || !st1 && st2) {
-            if (circles.back().size() == 0) {
-                circles.back().resize(1);
-            }
-            circles.back().back().push_back(i);
-        } else {
-            circles.back().resize(circles.back().size() + 1);
-            circles.back().back().push_back(i);
-            st2 = false;
-        }
-    }
 
     /*groups.resize(groups.size() + 1);
     int itemp = circles.back().size() / 2;
@@ -423,7 +436,7 @@ Spirals::finishAnalysis(int /*nframes*/)
 
     file = std::fopen("steps_Rspiral_Nmonomers.txt", "w+");
     for (int i = 0; i < circles.size(); i++) {
-        std::fprintf(file, "Frame # %3.2d\n", i);
+        std::fprintf(file, "Frame # %6.2d\n", i);
         std::fprintf(file, "Spiral steps:\n");
         for (int j = 0; j < circles[i].size(); j++) {
             std::fprintf(file, "%3.2f ", std::sqrt(distance2(kernel[i].krnl[circles[i][j].front() - 1], kernel[i].krnl[circles[i][j].back() - 1])));