sweep

diff --git a/src/spirals.cpp b/src/spirals.cpp
index 375ef3d75f86ae954270388a7c47c8a08b895b4f..b4801326cf9bbc4e67ba2d11759ed64e83aab1ee 100644 (file)
--- a/src/spirals.cpp
+++ b/src/spirals.cpp
@@ -153,9 +153,7 @@ long double fp3 (long double x0, long double y0, long double z0, long double p1,
 }
 
 void linear_kernel_search (long double &x0, long double &y0, long double &z0, long double &p1, long double &p2, long double &p3, std::vector< RVec > x, long double epsi) {
-    long double FX = 0, FX0 = 0, FY0 = 0, FZ0 = 0, FP1 = 0, FP2 = 0, FP3 = 0;
-    long double L0;
-    //int count = 0;
+    long double Ftemp = 0, FX = 0, FX0 = 0, FY0 = 0, FZ0 = 0, FP1 = 0, FP2 = 0, FP3 = 0, L0 = 0;
     while (true) {
         FX = Fx(x0, y0, z0, p1, p2, p3, x);
         FX0 = fx0(x0, y0, z0, p1, p2, p3, x);
@@ -166,28 +164,28 @@ void linear_kernel_search (long double &x0, long double &y0, long double &z0, lo
         FP3 = fp3(x0, y0, z0, p1, p2, p3, x);
 
         L0 = 1;
-        while (Fx(x0 - L0 * FX0, y0 - L0 * FY0, z0 - L0 * FZ0, p1 - L0 * FP1, p2 - L0 * FP2, p3 - L0 * FP3, x) - FX > 0) {
-            L0 /= 2;
-            if ((x0 - L0 * FX0 < epsi) && (y0 - L0 * FY0 < epsi) && (z0 - L0 * FZ0 < epsi) && (p1 - L0 * FP1 < epsi) && (p2 - L0 * FP2 < epsi) && (p3 - L0 * FP3 < epsi)) {
-                L0 = 0;
+        while (true) {
+            Ftemp = Fx(x0 - L0 * FX0, y0 - L0 * FY0, z0 - L0 * FZ0, p1 - L0 * FP1, p2 - L0 * FP2, p3 - L0 * FP3, x);
+            if (Ftemp - FX > 0) {
+                L0 /= 2;
+            } else {
+                x0 -= L0 * FX0;
+                y0 -= L0 * FY0;
+                z0 -= L0 * FZ0;
+                p1 -= L0 * FP1;
+                p2 -= L0 * FP2;
+                p3 -= L0 * FP3;
+                if ((L0 * FX0 < pow(epsi, 2)) && (L0 * FY0 < pow(epsi, 2)) && (L0 * FZ0 < pow(epsi, 2)) && (L0 * FP1 < pow(epsi, 2)) && (L0 * FP2 < pow(epsi, 2)) && (L0 * FP3 < pow(epsi, 2))) {
+                    L0 = 0;
+                }
+                break;
             }
         }
 
-        //std::cout << FX - Fx(x0 - L0 * FX0, y0 - L0 * FY0, z0 - L0 * FZ0, p1 - L0 * FP1, p2 - L0 * FP2, p3 - L0 * FP3, x) << "\n";
-        if (FX - Fx(x0 - L0 * FX0, y0 - L0 * FY0, z0 - L0 * FZ0, p1 - L0 * FP1, p2 - L0 * FP2, p3 - L0 * FP3, x) > epsi) {
-            x0 -= L0 * FX0;
-            y0 -= L0 * FY0;
-            z0 -= L0 * FZ0;
-            p1 -= L0 * FP1;
-            p2 -= L0 * FP2;
-            p3 -= L0 * FP3;
-        } else {
+        if (L0 = 0) {
             break;
         }
-        //count++;
     }
-    //std::cout << "\n";
-    //std::cout << count << "\n";
 }
 
 RVec kernel_pro (double x0, double y0, double z0, double p1, double p2, double p3, RVec x) {
@@ -200,8 +198,12 @@ RVec kernel_pro (double x0, double y0, double z0, double p1, double p2, double p
 }
 
 double left_right_turn (RVec a, RVec b, RVec c) {
-    return a[0] * b[1] * c[2] + a[1] * b[2] * c[0] + b[0] * c[1] * a[2] -
-                (a[2] * b[1] * c[0] + a[0] * b[2] * c[1] + a[1] * b[0] * c[2]);
+    return  a[0] * b[1] * c[2] +
+            a[1] * b[2] * c[0] +
+            b[0] * c[1] * a[2] -
+            a[2] * b[1] * c[0] -
+            a[0] * b[2] * c[1] -
+            a[1] * b[0] * c[2];
 }
 
 /*! \brief
@@ -237,7 +239,7 @@ class Spirals : public TrajectoryAnalysisModule
 
         SelectionList                           sel_;
         int                                     frames      = 0;
-        double                                  epsi        = 0.00001;
+        double                                  epsi        = 0.001;
 
         std::vector< std::vector< RVec > >                  monomers;
         std::vector< kernel_maxima >                        kernel;
@@ -270,8 +272,6 @@ Spirals::initOptions(IOptionsContainer          *options,
                            .description("Position to calculate distances for"));
 }
 
-// -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/test.ndx' -on '/home/toluk/Develop/samples/reca_rd/core.ndx'
-
 void
 Spirals::initAnalysis(const TrajectoryAnalysisSettings &settings,
                                const TopologyInformation         & /*top*/)
@@ -323,11 +323,8 @@ Spirals::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
     t5 = arrow[1];
     t6 = arrow[2];
 
-    //linear_kernel_search(mid[0], mid[1], mid[2], arrow[0], arrow[1], arrow[2], temp, epsi); //изменить формат функции для изменения значений переменных
     linear_kernel_search(t1, t2, t3, t4, t5, t6, temp, epsi);
 
-    //std::cout << t1 << " " << t2 << " " << t3 << " " << t4 << " " << t5 << " " << t6 << "\n";
-
     mid[0] = t1;
     mid[1] = t2;
     mid[2] = t3;
@@ -379,8 +376,6 @@ Spirals::finishAnalysis(int /*nframes*/)
 {
     FILE *file;
 
-
-
     file = std::fopen("linear_kernel.txt", "w+");
     for (int i = 0; i < kernel.size(); i++) {
         for (int j = 0; j < monomers[i].size(); j++) {
@@ -394,8 +389,6 @@ Spirals::finishAnalysis(int /*nframes*/)
     }
     std::fclose(file);
 
-
-
     file = std::fopen("circles_points.txt", "w+");
     for (int i = 0; i < circles.size(); i++) {
         for (int j = 0; j < circles[i].size(); j++) {