- merged with FitNew

- switched full corr output place in code

diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp
index f930fb4d8fbaea087dedbd0965e8c2e0df290f37..6fa1cbec19f4b8ec1fe3d6c84d59ab3ba01c5d9a 100644 (file)
--- a/src/spacetimecorr.cpp
+++ b/src/spacetimecorr.cpp
@@ -49,12 +49,255 @@
 #include <gromacs/utility/smalloc.h>
 #include <gromacs/math/do_fit.h>
 
 #include <gromacs/utility/smalloc.h>
 #include <gromacs/math/do_fit.h>
 

-#include "fittingn.cpp"
+#include <vector>
+#include <math.h>
+
+//#include "fittingn.cpp"

 
 using namespace gmx;
 
 using gmx::RVec;
 
 
 using namespace gmx;
 
 using gmx::RVec;
 

+double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
+    return  sqrt(   pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2) +
+                    pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2) +
+                    pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2)  );
+}
+
+void searchF0xyzabc (double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
+    double t01, t02, t03, t04, t05, t06, t07;
+    t01 = pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2);
+    t02 = pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2);
+    t03 = pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2);
+    t04 = sqrt(t01 + t02 + t03);
+    t05 = (aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x);
+    t06 = (aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y);
+    t07 = (aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x);
+
+    //#pragma omp parallel sections
+    //{
+        //#pragma omp section
+            F += t04;
+        //#pragma omp section
+            Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
+        //#pragma omp section
+            Fy += -(2 * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) * t07 - 2 * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) * t05 + 2 * cos(B) * sin(A) * t06) / (2 * t04);
+        //#pragma omp section
+            Fz += -(2 * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) * t05 - 2 * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) * t07 + 2 * cos(A) * cos(B) * t06) / (2 * t04);
+        //#pragma omp section
+            Fa += -(2 * (cos(A) * cos(B) * biy_plus_y - cos(B) * sin(A) * biz_plus_z) * t06 -
+                2 * (biy_plus_y * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + biz_plus_z * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C))) * t07 +
+                2 * (biy_plus_y * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) + biz_plus_z * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B))) * t05) / (2 * t04);
+        //#pragma omp section
+            Fb += -(2 * (cos(A) * cos(B) * sin(C) * biz_plus_z - sin(B) * sin(C) * bix_plus_x + cos(B) * sin(A) * sin(C) * biy_plus_y) * t07 +
+                2 * (cos(A) * cos(B) * cos(C) * biz_plus_z - cos(C) * sin(B) * bix_plus_x + cos(B) * cos(C) * sin(A) * biy_plus_y) * t05 -
+                2 * (cos(B) * bix_plus_x + sin(A) * sin(B) * biy_plus_y + cos(A) * sin(B) * biz_plus_z) * t06) / (2 * t04);
+        //#pragma omp section
+            Fc += (2 * (biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * bix_plus_x) * t05 -
+                2 * (biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * bix_plus_x) * t07) / (2 * t04);
+    //}
+}
+
+void searchL (double &fLd, double &fLdd, double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C, double L, double fx, double fy, double fz, double fa, double fb, double fc) {
+    double AmLfa, BmLfb, CmLfc, BIXpXmLfx, BIYpYmLfy, BIZpZmLfz, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17, t18, t19, t20;
+    AmLfa = A - L * fa;
+    BmLfb = B - L * fb;
+    CmLfc = C - L * fc;
+    BIXpXmLfx = bix + x - L * fx;
+    BIYpYmLfy = biy + y - L * fy;
+    BIZpZmLfz = biz + z - L * fz;
+    t1 = (cos(AmLfa) * cos(CmLfc) + sin(AmLfa) * sin(BmLfb) * sin(CmLfc));
+    t2 = (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc));
+    t3 = (cos(AmLfa) * sin(CmLfc) - cos(CmLfc) * sin(AmLfa) * sin(BmLfb));
+    t4 = (sin(AmLfa) * sin(CmLfc) + cos(AmLfa) * cos(CmLfc) * sin(BmLfb));
+    t5 = pow(aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy, 2);
+    t6 = pow(aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx, 2);
+    t7 = pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2);
+
+    t8 = fc * sin(AmLfa) * sin(CmLfc);
+    t9 = fa * cos(AmLfa) * cos(CmLfc);
+    t10 = fb * cos(AmLfa) * cos(BmLfb) * sin(CmLfc);
+    t11 = fc * cos(AmLfa) * cos(CmLfc) * sin(BmLfb);
+    t12 = fa * sin(AmLfa) * sin(BmLfb) * sin(CmLfc);
+
+    t13 = fa * cos(AmLfa) * sin(BmLfb) * sin(CmLfc);
+    t14 = fc * cos(AmLfa) * sin(CmLfc);
+    t15 = fa * cos(CmLfc) * sin(AmLfa);
+    t16 = fb * cos(BmLfb) * sin(AmLfa) * sin(CmLfc);
+    t17 = fc * cos(CmLfc) * sin(AmLfa) * sin(BmLfb);
+
+    t18 = fx * cos(BmLfb) * sin(CmLfc);
+    t19 = fc * cos(BmLfb) * cos(CmLfc);
+    t20 = fb * sin(BmLfb) * sin(CmLfc);
+
+    fLd +=  (2 * (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
+                (   BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
+                2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
+                (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) -
+                    fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy + fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
+                2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
+                (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
+                    BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
+                    fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * t3 +
+                    fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) /
+            (2 * sqrt(  t5 + t6 + pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2)));
+
+
+    fLdd += (   (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
+                (2 * fb * fx * cos(BmLfb) - pow(fb, 2) * sin(BmLfb) * BIXpXmLfx - 2 * fa * fy * cos(AmLfa) * cos(BmLfb) + 2 * fa * fz * cos(BmLfb) * sin(AmLfa) + 2 * fb * fz * cos(AmLfa) * sin(BmLfb) +
+                    2 * fb * fy * sin(AmLfa) * sin(BmLfb) + pow(fa, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz + pow(fb, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz +
+                    pow(fa, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + pow(fb, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + 2 * fa * fb * cos(AmLfa) * sin(BmLfb) * BIYpYmLfy -
+                    2 * fa * fb * sin(AmLfa) * sin(BmLfb) * BIZpZmLfz) +
+                (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
+                    fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
+                    BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
+                    fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * (cos(AmLfa) * sin(CmLfc) -
+                    cos(CmLfc) * sin(AmLfa) * sin(BmLfb)) + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) -
+                    fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
+                (2 * BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) +
+                    fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) -
+                    fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) + 2 * BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) +
+                    t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
+                    2 * fy * t3 + 2 * fz * t4 +
+                    2 * fx * cos(BmLfb) * cos(CmLfc) - 2 * fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - 2 * fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
+                (2 * aix + 2 * t3 * BIYpYmLfy - 2 * (sin(AmLfa) * sin(CmLfc) +
+                    cos(AmLfa) * cos(CmLfc) * sin(BmLfb)) * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
+                (BIZpZmLfz * (pow(fa, 2) * sin(AmLfa) * sin(CmLfc) +
+                    pow(fc, 2) * sin(AmLfa) * sin(CmLfc) - 2 * fa * fc * cos(AmLfa) * cos(CmLfc) + pow(fa, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) +
+                    pow(fb, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) + pow(fc, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) + 2 * fa * fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) +
+                    2 * fb * fc * cos(AmLfa) * cos(BmLfb) * sin(CmLfc) - 2 * fa * fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) + BIYpYmLfy * (pow(fa, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) -
+                    pow(fc, 2) * cos(AmLfa) * sin(CmLfc) - 2 * fa * fc * cos(CmLfc) * sin(AmLfa) - pow(fa, 2) * cos(AmLfa) * sin(CmLfc) + pow(fb, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) +
+                    pow(fc, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) - 2 * fa * fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) + 2 * fa * fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) +
+                    2 * fb * fc * cos(BmLfb) * sin(AmLfa) * sin(CmLfc)) - 2 * fz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
+                    fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
+                    2 * fy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) -
+                    fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) + 2 * fb * fx * cos(CmLfc) * sin(BmLfb) + 2 * fc * t18 + pow(fb, 2) * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx +
+                    pow(fc, 2) * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx - 2 * fb * fc * sin(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
+                (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) *
+                (2 * BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + 2 * fy * t1 - 2 * fz * t2 + 2 * BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + 2 * t18 + 2 * t19 * BIXpXmLfx - 2 * t20 * BIXpXmLfx) +
+                (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
+                    fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) *
+                (2 * fx * sin(BmLfb) - 2 * fy * cos(BmLfb) * sin(AmLfa) + 2 * fb * cos(BmLfb) * BIXpXmLfx - 2 * fz * cos(AmLfa) * cos(BmLfb) - 2 * fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
+                    2 * fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + 2 * fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + 2 * fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
+                (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
+                (BIZpZmLfz * (pow(fa, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) -
+                    pow(fc, 2) * cos(CmLfc) * sin(AmLfa) - 2 * fa * t14 - pow(fa, 2) * cos(CmLfc) * sin(AmLfa) + pow(fb, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) +
+                    pow(fc, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) - 2 * fb * fc * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) + 2 * fa * t16 +
+                    2 * fa * t17) + BIYpYmLfy * (pow(fa, 2) * cos(AmLfa) * cos(CmLfc) + pow(fc, 2) * cos(AmLfa) * cos(CmLfc) - 2 * fa * t8 +
+                    pow(fa, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) + pow(fb, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) + pow(fc, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) -
+                    2 * fa * t10 - 2 * fa * t11 - 2 * fb * t19 * sin(AmLfa)) - 2 * fz * (t8 - t9 + t10 + t11 - t12) - 2 * fy * (t13 - t14 - t15 + t16 + t17) -
+                    2 * fc * fx * cos(BmLfb) * cos(CmLfc) + 2 * fb * fx * sin(BmLfb) * sin(CmLfc) +
+                    pow(fb, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + pow(fc, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + 2 * fb * fc * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx)) /
+            (2 * sqrt(  t5 + t6 + t7)) -
+            (   (2 *    (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
+                        (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
+                        2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
+                        (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
+                            fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
+                        2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
+                        (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
+                        BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
+                        fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) *
+                (   (2 * aix + 2 * t3 * BIYpYmLfy - 2 * t4 * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
+                    (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
+                        BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
+                        fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
+                    (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * t2 * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
+                    (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * (cos(AmLfa) * cos(CmLfc) + sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) +
+                        t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) - (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
+                    (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
+                        fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy))) /
+            (4 * pow    (   t5 + t6 + t7, 3 / 2));
+
+}
+
+
+void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
+    RVec temp;
+    //#pragma omp parallel
+    //{
+        //#pragma omp for schedule(dynamic)
+        for (int i = 0; i < b.size(); i++) {
+            temp = b[i];
+            b[i][0] = (temp[2] + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - (temp[1] + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * (temp[0] + x);
+            b[i][1] = (temp[1] + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - (temp[2] + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * (temp[0] + x);
+            b[i][2] = cos(A) * cos(B) * (temp[2] + z) - sin(B) * (temp[0] + x) + cos(B) * sin(A) * (temp[1] + y);
+        }
+    //}
+}
+
+void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< int, int > > pairs) {
+    midA[0] = 0;
+    midA[1] = 0;
+    midA[2] = 0;
+
+    midB[0] = 0;
+    midB[1] = 0;
+    midB[2] = 0;
+
+    for (int i = 0; i < pairs.size(); i++) {
+        rvec_inc(midA, a[pairs[i].first]);
+        rvec_inc(midB, b[pairs[i].second]);
+    }
+    midA[0] /= pairs.size();
+    midA[1] /= pairs.size();
+    midA[2] /= pairs.size();
+
+    midB[0] /= pairs.size();
+    midB[1] /= pairs.size();
+    midB[2] /= pairs.size();
+}
+
+void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< int, int > > pairs) {
+    double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
+    RVec ma, mb;
+    CalcMid(a, b, ma, mb, pairs);
+    rvec_dec(ma, mb);
+    double x = ma[0], y = ma[1], z = ma[2], A = 0, B = 0, C = 0;
+    for (int i = 0; i < pairs.size(); i++) {
+        f1 += F(a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + x, b[pairs[i].second][1] + y, b[pairs[i].second][2] + z, A, B, C);
+    }
+    if (f1 == 0) {
+        return;
+    } else {
+        double fLd, fLdd;
+        //l = 0;
+        int count = 0;
+            while (f1 - f2 < 0 || f1 - f2 > 0.00001) {
+                f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
+                for (int i = 0; i < pairs.size(); i++) {
+                    searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc, a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + x, b[pairs[i].second][1] + y, b[pairs[i].second][2] + z, A, B, C);
+                }
+                /*do {
+                    fLd = 0;
+                    fLdd = 0;
+                    for (int i = 0; i < pairs.size(); i++) {
+                        searchL(fLd, fLdd, a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0], b[pairs[i].second][1], b[pairs[i].second][2], x, y, z, A, B, C, l, fx, fy, fz, fa, fb, fc);
+                    }
+                    l -= fLd / fLdd;
+                } while (std::abs(fLd) > 1);*/
+                //std::cout << " Ltempo= " << Ltempo << " iter= " << Lco << "\n";
+                while (true) {
+                    f2 = 0;
+                    for (int i = 0; i < pairs.size(); i++) {
+                        f2 += F(a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + (x - l * fx), b[pairs[i].second][1] + (y - l * fy), b[pairs[i].second][2] + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
+                    }
+                    count++;
+                    if (f2 < f1) {
+                        x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
+                        fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
+                        break;
+                    } else {
+                        l *= 0.85;
+                    }
+                }
+                count++;
+            }
+        //std::cout << "iteration count: " << count << "\n";
+        ApplyFit(b, x, y, z, A, B, C);
+    }
+}
+

 void make_pdb_trajectory(std::vector< std::vector< RVec > > trajectory, const char* file_name)
 {
     FILE *file;
 void make_pdb_trajectory(std::vector< std::vector< RVec > > trajectory, const char* file_name)
 {
     FILE *file;


@@ -68,7 +311,7 @@ void make_pdb_trajectory(std::vector< std::vector< RVec > > trajectory, const ch
     }
 }
 
     }
 }
 

-void make_correlation_file(std::vector< std::vector< std::vector< double > > > correlations, const char* file_name, int start)
+void make_correlation_file(std::vector< std::vector< std::vector< float > > > correlations, const char* file_name, int start)

 {
     FILE *file;
     file = std::fopen(file_name, "w+");
 {
     FILE *file;
     file = std::fopen(file_name, "w+");


@@ -89,7 +332,7 @@ void make_correlation_file(std::vector< std::vector< std::vector< double > > > c
     std::fclose(file);
 }
 
     std::fclose(file);
 }
 

-void make_rout_file2(double crl_border, std::vector< int > indx, std::vector< std::vector< std::pair< int, int > > > rout, const char* file_name)
+void make_rout_file2(float crl_border, std::vector< int > indx, std::vector< std::vector< std::pair< int, int > > > rout, const char* file_name)

 {
     FILE *file;
     file = std::fopen(file_name, "w+");
 {
     FILE *file;
     file = std::fopen(file_name, "w+");


@@ -103,7 +346,7 @@ void make_rout_file2(double crl_border, std::vector< int > indx, std::vector< st
     std::fclose(file);
 }
 
     std::fclose(file);
 }
 

-void make_best_corrs_graphics(std::vector< std::vector< std::vector< double > > > correlations,
+void make_best_corrs_graphics(std::vector< std::vector< std::vector< float > > > correlations,

                               std::vector< std::vector< std::pair< int, int > > > rout_pairs,
                               std::vector< int > indx,
                               const char* file_name)
                               std::vector< std::vector< std::pair< int, int > > > rout_pairs,
                               std::vector< int > indx,
                               const char* file_name)


@@ -143,7 +386,7 @@ bool isitsubset (std::vector< int > a, std::vector< int > b) {
     }
 }
 
     }
 }
 

-void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame, std::vector< std::vector< std::vector< double > > > &crl, int tauS, int tauE) {
+void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame, std::vector< std::vector< std::vector< float > > > &crl, int tauS, int tauE) {

     crl.resize(tauE + 1);
     for (int i = 0; i < crl.size(); i++) {
         crl[i].resize(traj.front().size());
     crl.resize(tauE + 1);
     for (int i = 0; i < crl.size(); i++) {
         crl[i].resize(traj.front().size());


@@ -155,12 +398,16 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame
     temp_zero[0] = 0;
     temp_zero[1] = 0;
     temp_zero[2] = 0;
     temp_zero[0] = 0;
     temp_zero[1] = 0;
     temp_zero[2] = 0;

-    #pragma omp parallel shared(crl)
+
+    std::vector< float > d;
+    d.resize(traj.front().size(), 0);
+
+    #pragma omp parallel shared(crl, temp_zero, d)

     {
         #pragma omp for schedule(dynamic)
         for (int i = tauS; i <= tauE; i += 1) {
     {
         #pragma omp for schedule(dynamic)
         for (int i = tauS; i <= tauE; i += 1) {

-            real tmp2;
-            rvec *medx, *medy, temp, temp1, temp2;
+            rvec *medx, *medy, temp1, temp2;
+            real tmp;

             snew(medx, traj.front().size());
             snew(medy, traj.front().size());
             for (int i = 0; i < traj.front().size(); i++) {
             snew(medx, traj.front().size());
             snew(medy, traj.front().size());
             for (int i = 0; i < traj.front().size(); i++) {


@@ -169,39 +416,30 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame
             }
             for (int j = 0; j < traj.size() - i - 1; j++) {
                 for (int f = 0; f < traj.front().size(); f++) {
             }
             for (int j = 0; j < traj.size() - i - 1; j++) {
                 for (int f = 0; f < traj.front().size(); f++) {

-                    rvec_sub(traj[b_frame][f], traj[j][f], temp);
-                    rvec_inc(medx[f], temp);
+                    rvec_sub(traj[b_frame][f], traj[j][f], temp1);
+                    rvec_inc(medx[f], temp1);

 
 

-                    rvec_sub(traj[b_frame][f], traj[j + i][f], temp);
-                    rvec_inc(medy[f], temp);
+                    rvec_sub(traj[b_frame][f], traj[j + i][f], temp1);
+                    rvec_inc(medy[f], temp1);

                 }
             }
                 }
             }

-
-            for (int j = 0; j < traj.front().size(); j++) {
-                tmp2 = 1 / (traj.size() - 1 - i);
-
-                /*temp2 = traj.size() - 1;
-                temp2 -= i;
-                temp2 = 1 / temp2;*/
-
-                copy_rvec(medx[j], temp);
-                svmul(tmp2, temp, medx[j]);
-                copy_rvec(medy[j], temp);
-                svmul(tmp2, temp, medy[j]);
-            }
-            std::vector< std::vector< double > > a, b, c;
-            a.resize(traj.front().size());
-            b.resize(traj.front().size());
-            c.resize(traj.front().size());

             for (int j = 0; j < traj.front().size(); j++) {
             for (int j = 0; j < traj.front().size(); j++) {

-                a[j].resize(traj.front().size(), 0);
-                b[j].resize(traj.front().size(), 0);
-                c[j].resize(traj.front().size(), 0);
+                tmp = traj.size() - 1;
+                tmp -= i;
+                tmp = 1 / tmp;
+                //tmp = 1 / (traj.size() - 1 - i);
+                copy_rvec(medx[j], temp1);
+                svmul(tmp, temp1, medx[j]);
+                copy_rvec(medy[j], temp1);
+                svmul(tmp, temp1, medy[j]);

             }
             }

+            std::vector< std::vector< float > > a, b, c;
+            a.resize(traj.front().size(), d);
+            b.resize(traj.front().size(), d);
+            c.resize(traj.front().size(), d);

             for (int j = 0; j < traj.size() - i - 1; j++) {
                 for (int f1 = 0; f1 < traj.front().size(); f1++) {
                     for (int f2 = 0; f2 < traj.front().size(); f2++) {
             for (int j = 0; j < traj.size() - i - 1; j++) {
                 for (int f1 = 0; f1 < traj.front().size(); f1++) {
                     for (int f2 = 0; f2 < traj.front().size(); f2++) {

-

                         rvec_sub(traj[b_frame][f1], traj[j][f1], temp1);
                         rvec_dec(temp1, medx[f1]);
 
                         rvec_sub(traj[b_frame][f1], traj[j][f1], temp1);
                         rvec_dec(temp1, medx[f1]);
 


@@ -221,7 +459,7 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame
             }
             sfree(medx);
             sfree(medy);
             }
             sfree(medx);
             sfree(medy);

-            std::cout << i << "\n";
+            std::cout << i << " corr done\n";

         }
     }
     #pragma omp barrier
         }
     }
     #pragma omp barrier


@@ -292,9 +530,9 @@ void domain_chopping(SelectionList seList, std::vector< int > indx, std::vector<
     }
 }
 
     }
 }
 

-void graph_calculation(std::vector< std::vector< std::pair< double, int > > > &graph, std::vector< std::vector< int > > &s_graph, std::vector< std::vector< std::pair< int, int > > > &s_graph_rbr,
+void graph_calculation(std::vector< std::vector< std::pair< float, int > > > &graph, std::vector< std::vector< int > > &s_graph, std::vector< std::vector< std::pair< int, int > > > &s_graph_rbr,

                       std::vector< std::vector< RVec > > traj, int b_frame,
                       std::vector< std::vector< RVec > > traj, int b_frame,

-                      std::vector< std::vector< std::vector< double > > > crl, double crl_b, double e_rad, int tauE) {
+                      std::vector< std::vector< std::vector< float > > > crl, float crl_b, float e_rad, int tauE) {

     graph.resize(traj.front().size());
     for (int i = 0; i < traj.front().size(); i++) {
         graph[i].resize(traj.front().size(), std::make_pair(0, -1));
     graph.resize(traj.front().size());
     for (int i = 0; i < traj.front().size(); i++) {
         graph[i].resize(traj.front().size(), std::make_pair(0, -1));


@@ -316,6 +554,7 @@ void graph_calculation(std::vector< std::vector< std::pair< double, int > > > &g
             }
         }
     }
             }
         }
     }

+    std::cout << "crl analysed\n";

     std::vector< bool > graph_flags;
     graph_flags.resize(traj.front().size(), true);
     std::vector< int > a;
     std::vector< bool > graph_flags;
     graph_flags.resize(traj.front().size(), true);
     std::vector< int > a;


@@ -361,15 +600,15 @@ void graph_calculation(std::vector< std::vector< std::pair< double, int > > > &g
     }
 }
 
     }
 }
 

-bool myfunction (const std::pair< int, double > i, const std::pair< int, double > j) {
+bool myfunction (const std::pair< int, float > i, const std::pair< int, float > j) {

     return i.second < j.second;
 }
 
 void graph_back_bone_evaluation(std::vector< std::vector < std::pair< int, int > > > &rout_n, int indxSize,
     return i.second < j.second;
 }
 
 void graph_back_bone_evaluation(std::vector< std::vector < std::pair< int, int > > > &rout_n, int indxSize,

-                                std::vector< std::vector< std::pair< double, int > > > graph, std::vector< std::vector< int > > s_graph, std::vector< std::vector< std::pair< int, int > > > s_graph_rbr) {
-    std::vector< double > key;
+                                std::vector< std::vector< std::pair< float, int > > > graph, std::vector< std::vector< int > > s_graph, std::vector< std::vector< std::pair< int, int > > > s_graph_rbr) {
+    std::vector< float > key;

     std::vector< int > path;
     std::vector< int > path;

-    std::vector< std::pair< int, double > > que;
+    std::vector< std::pair< int, float > > que;

     std::vector< std::pair< int, int > > a;
     int v;
     for (int i = 0; i < s_graph.size(); i++) {
     std::vector< std::pair< int, int > > a;
     int v;
     for (int i = 0; i < s_graph.size(); i++) {


@@ -473,9 +712,8 @@ class Domains : public TrajectoryAnalysisModule
         int                                                         frames              = 0;
         int                                                         basic_frame         = 0;
         int                                                         tau                 = 0;
         int                                                         frames              = 0;
         int                                                         basic_frame         = 0;
         int                                                         tau                 = 0;

-        int                                                         graph_tau           = 0;
-        double                                                      crl_border          = 0;
-        double                                                      eff_rad             = 1.5;
+        float                                                      crl_border          = 0;
+        float                                                      eff_rad             = 1.5;

         real                                                        **w_rls;
         std::vector< std::vector< std::pair< int, int > > >         w_rls2;
 
         real                                                        **w_rls;
         std::vector< std::vector< std::pair< int, int > > >         w_rls2;
 


@@ -504,20 +742,16 @@ Domains::initOptions(IOptionsContainer          *options,
     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
                             .store(&fnNdx_).defaultBasename("domains")
                             .description("Index file from the domains"));
     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
                             .store(&fnNdx_).defaultBasename("domains")
                             .description("Index file from the domains"));

-    // Add option for graph_tau constant
-    options->addOption(gmx::IntegerOption("Gtau")
-                            .store(&graph_tau)
-                            .description("number of frames for graph to see for time travel"));

     // Add option for tau constant
     options->addOption(gmx::IntegerOption("tau")
                             .store(&tau)
                             .description("number of frames for time travel"));
     // Add option for crl_border constant
     // Add option for tau constant
     options->addOption(gmx::IntegerOption("tau")
                             .store(&tau)
                             .description("number of frames for time travel"));
     // Add option for crl_border constant

-    options->addOption(DoubleOption("crl")
+    options->addOption(FloatOption("crl")

                             .store(&crl_border)
                             .description("make graph based on corrs > constant"));
     // Add option for eff_rad constant
                             .store(&crl_border)
                             .description("make graph based on corrs > constant"));
     // Add option for eff_rad constant

-    options->addOption(DoubleOption("ef_rad")
+    options->addOption(FloatOption("ef_rad")

                             .store(&eff_rad)
                             .description("effective radius for atoms to evaluate corrs"));
     // Add option for selection list
                             .store(&eff_rad)
                             .description("effective radius for atoms to evaluate corrs"));
     // Add option for selection list


@@ -558,6 +792,8 @@ Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
     frankenstein_trajectory.resize(trajectory.size());
     frankenstein_trajectory.back() = trajectory.back();
     trajectory.resize(trajectory.size() + 1);
     frankenstein_trajectory.resize(trajectory.size());
     frankenstein_trajectory.back() = trajectory.back();
     trajectory.resize(trajectory.size() + 1);

+    trajectory.back().resize(0);
+    trajectory.back().resize(index.size());

 
     w_rls2.resize(domains_local.size() + 1);
     for (int i = 0; i < domains_local.size(); i++) {
 
     w_rls2.resize(domains_local.size() + 1);
     for (int i = 0; i < domains_local.size(); i++) {


@@ -574,27 +810,30 @@ Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
 }
 
 // -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/test6.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelectionList5'  -tau 5 -crl 0.10
 }
 
 // -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/test6.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelectionList5'  -tau 5 -crl 0.10

+// -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/CorrsTestDomainsNfit.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelectionListDomainsNFit' -tau 5000 -crl 0.75 -ef_rad 1
+// -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/TestCa.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelListCa' -tau 100 -crl 0.75 -ef_rad 1

 
 void
 Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
                       TrajectoryAnalysisModuleData *pdata)
 {
 
 void
 Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
                       TrajectoryAnalysisModuleData *pdata)
 {

-    trajectory.back().resize(0);
-    trajectory.back().resize(index.size());

     for (int i = 0; i < index.size(); i++) {
         trajectory.back()[i] = fr.x[index[i]];
     }
     frankenstein_trajectory.resize(frankenstein_trajectory.size() + 1);
     frankenstein_trajectory.back().resize(index.size());
     for (int i = 0; i < index.size(); i++) {
         trajectory.back()[i] = fr.x[index[i]];
     }
     frankenstein_trajectory.resize(frankenstein_trajectory.size() + 1);
     frankenstein_trajectory.back().resize(index.size());

-    for (int i = 0; i < domains_local.size(); i++) {
-        std::vector< RVec > basic, traj;
-        basic = trajectory[basic_frame];
-        traj = trajectory.back();
-        MyFitNew(basic, traj, w_rls2[i]);
-        for (int j = 0; j < domains_local[i].size(); j++) {
-            frankenstein_trajectory.back()[domains_local[i][j]] = traj[domains_local[i][j]];
+    #pragma omp parallel shared(frankenstein_trajectory, trajectory)
+    {
+        #pragma omp for schedule(dynamic)
+        for (int i = 0; i < domains_local.size(); i++) {
+            std::vector< RVec > traj = trajectory.back();
+            MyFitNew(trajectory[basic_frame], traj, w_rls2[i]);
+            for (int j = 0; j < domains_local[i].size(); j++) {
+                frankenstein_trajectory.back()[domains_local[i][j]] = traj[domains_local[i][j]];
+            }

         }
     }
         }
     }

+    #pragma omp barrier

     std::cout << "frame № " << frames + 1 <<" analysed\n";
     frames++;
 }
     std::cout << "frame № " << frames + 1 <<" analysed\n";
     frames++;
 }


@@ -602,7 +841,7 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 void
 Domains::finishAnalysis(int nframes)
 {
 void
 Domains::finishAnalysis(int nframes)
 {

-    std::vector< std::vector< std::vector< double > > > crltns;
+    std::vector< std::vector< std::vector< float > > > crltns;

     int k = 1000, m = 0;
     if (tau > 0) {
         k = tau;
     int k = 1000, m = 0;
     if (tau > 0) {
         k = tau;


@@ -616,6 +855,7 @@ Domains::finishAnalysis(int nframes)
 
     std::cout << "Correlation's evaluation - start\n\n";
     correlation_evaluation(frankenstein_trajectory, basic_frame, crltns, m, k);
 
     std::cout << "Correlation's evaluation - start\n\n";
     correlation_evaluation(frankenstein_trajectory, basic_frame, crltns, m, k);

+    make_correlation_file(crltns, "full crltns file Ca.txt", 0);

     std::cout << "Correlation's evaluation - end\n\n";
 
     /*
     std::cout << "Correlation's evaluation - end\n\n";
 
     /*


@@ -630,7 +870,7 @@ Domains::finishAnalysis(int nframes)
         for (int i2 = 1; i2 < crltns.size(); i2++) {
             for (int i3 = 0; i3 < crltns[i2].size(); i3++) {
                 for (int i4 = 0; i4 < crltns[i2][i3].size(); i4++) {
         for (int i2 = 1; i2 < crltns.size(); i2++) {
             for (int i3 = 0; i3 < crltns[i2].size(); i3++) {
                 for (int i4 = 0; i4 < crltns[i2][i3].size(); i4++) {

-                    if (std::abs(crltns[i2][i3][i4]) >= (double)i1 / 100 && i3 != i4) {
+                    if (std::abs(crltns[i2][i3][i4]) >= (float)i1 / 100 && i3 != i4) {

                         i5++;
                     }
                 }
                         i5++;
                     }
                 }


@@ -644,12 +884,12 @@ Domains::finishAnalysis(int nframes)
 
     std::cout << "graph evaluation: start\n";
 
 
     std::cout << "graph evaluation: start\n";
 

-    std::vector< std::vector< std::pair< double, int > > > graph;
+    std::vector< std::vector< std::pair< float, int > > > graph;

     std::vector< std::vector< int > > sub_graph;
     std::vector< std::vector< std::pair< int, int > > > sub_graph_rbr;
 
 
     std::vector< std::vector< int > > sub_graph;
     std::vector< std::vector< std::pair< int, int > > > sub_graph_rbr;
 
 

-    graph_calculation(graph, sub_graph, sub_graph_rbr, frankenstein_trajectory, basic_frame, crltns, crl_border, eff_rad, graph_tau/*k*/);
+    graph_calculation(graph, sub_graph, sub_graph_rbr, frankenstein_trajectory, basic_frame, crltns, crl_border, eff_rad, k);

 
 
     std::vector< std::vector < std::pair< int, int > > > rout_new;
 
 
     std::vector< std::vector < std::pair< int, int > > > rout_new;


@@ -661,8 +901,11 @@ Domains::finishAnalysis(int nframes)
 
     std::cout << "routs evaluation: end\n";
 
 
     std::cout << "routs evaluation: end\n";
 

-    make_rout_file2(crl_border, index, rout_new, "Routs_DomainsNFit_5000_0.70_1_t.txt");
-    make_best_corrs_graphics(crltns, rout_new, index, "best_graphics_DomainsNFit_5000_0.70_1_t.txt");
+
+    // Routs_DomainsNFit_5000_0.75_1_t.txt
+    make_rout_file2(crl_border, index, rout_new, "Routs_Ca.txt");
+    make_best_corrs_graphics(crltns, rout_new, index, "best_graphics_Ca.txt");
+

     std::cout << "Finish Analysis - end\n\n";
 }
 
     std::cout << "Finish Analysis - end\n\n";
 }