added missing header file

diff --git a/include/newfit.h b/include/newfit.h
new file mode 100644 (file)
index 0000000..0cd4ee0
--- /dev/null
+++ b/include/newfit.h
@@ -0,0 +1,173 @@
+#ifndef NEWFIT_H
+#define NEWFIT_H
+
+#include "gromacs/math/vectypes.h"
+#include "vector"
+
+using gmx::RVec;
+
+namespace NewFit {
+
+    void implement_fit(std::vector< RVec > &if_frame, RVec if_step, double if_A, double if_B, double if_C) {
+        for (int i = 0; i < if_frame.size(); i++) {
+            if_frame[i] += if_step;
+        }
+        for (int i = 0; i < if_frame.size(); i++) {
+            if_step[0] = if_frame[i][0];
+            if_step[1] = if_frame[i][1];
+            if_step[2] = if_frame[i][2];
+            if_frame[i][1] = if_step[1] * cos(if_A) - if_step[2] * sin(if_A);
+            if_frame[i][2] = if_step[1] * sin(if_A) + if_step[2] * cos(if_A);
+        }
+        for (int i = 0; i < if_frame.size(); i++) {
+            if_step[0] = if_frame[i][0];
+            if_step[1] = if_frame[i][1];
+            if_step[2] = if_frame[i][2];
+            if_frame[i][0] = if_step[0] * cos(if_B) + if_step[2] * sin(if_B);
+            if_frame[i][2] = - if_step[0] * sin(if_B) + if_step[2] * cos(if_B);
+        }
+        for (int i = 0; i < if_frame.size(); i++) {
+            if_step[0] = if_frame[i][0];
+            if_step[1] = if_frame[i][1];
+            if_step[2] = if_frame[i][2];
+            if_frame[i][0] = if_step[0] * cos(if_C) - if_step[1] * sin(if_C);
+            if_frame[i][1] = if_step[0] * sin(if_C) + if_step[1] * cos(if_C);
+        }
+    }
+
+    double F(double aix, double aiy, double aiz, double bix, double biy, double biz, double cA, double sA, double cB, double sB, double cC, double sC) {
+        return  gmx::square(aix + biy*(cA*sC - cC*sA*sB) - biz*(sA*sC + cA*cC*sB) - cB*cC*bix) +
+                gmx::square(aiy - biy*(cA*cC + sA*sB*sC) + biz*(cC*sA - cA*sB*sC) - cB*sC*bix) +
+                gmx::square(aiz + sB*bix - cA*cB*biz - cB*sA*biy) ;
+    }
+
+    double Fx(double temp1, double temp2, double temp3, double cB, double sB, double cC, double sC) {
+        return 2*sB*temp3 -
+               2*cB*cC*temp1 -
+               2*cB*sC*temp2 ;
+    }
+
+    double Fy(double temp1, double temp2, double temp3, double cA, double sA, double cB, double sB, double cC, double sC) {
+        return 2*(cA*sC - cC*sA*sB)*temp1 -
+               2*(cA*cC + sA*sB*sC)*temp2 -
+               2*cB*sA*temp3 ;
+    }
+
+    double Fz(double temp1, double temp2, double temp3, double cA, double sA, double cB, double sB, double cC, double sC) {
+        return 2*(cC*sA - cA*sB*sC)*temp2 -
+               2*(sA*sC + cA*cC*sB)*temp1 -
+               2*cA*cB*temp3 ;
+    }
+
+    double FA(double temp1, double temp2, double temp3, double biy_pl_y, double biz_pl_z, double cA, double sA, double cB, double sB, double cC, double sC) {
+        return 2*(biy_pl_y*(cC*sA - cA*sB*sC) + biz_pl_z*(cA*cC + sA*sB*sC))*temp2 -
+               2*(cA*cB*biy_pl_y - cB*sA*biz_pl_z)*temp3 -
+               2*(biy_pl_y*(sA*sC + cA*cC*sB) + biz_pl_z*(cA*sC - cC*sA*sB))*temp1 ;
+    }
+
+    double FB(double temp1, double temp2, double temp3, double bix_pl_x, double biy_pl_y, double biz_pl_z, double cA, double sA, double cB, double sB, double cC, double sC) {
+        return 2*(cB*bix_pl_x + sA*sB*biy_pl_y + cA*sB*biz_pl_z)*temp3 -
+               2*(cA*cB*cC*biz_pl_z - cC*sB*bix_pl_x + cB*cC*sA*biy_pl_y)*temp1 -
+               2*(cA*cB*sC*biz_pl_z - sB*sC*bix_pl_x + cB*sA*sC*biy_pl_y)*temp2 ;
+    }
+
+    double FC(double temp1, double temp2, double bix_pl_x, double biy_pl_y, double biz_pl_z, double cA, double sA, double cB, double sB, double cC, double sC) {
+        return 2*(biy_pl_y*(cA*cC + sA*sB*sC) - biz_pl_z*(cC*sA - cA*sB*sC) + cB*sC*bix_pl_x)*temp1 -
+               2*(biz_pl_z*(sA*sC + cA*cC*sB) - biy_pl_y*(cA*sC - cC*sA*sB) + cB*cC*bix_pl_x)*temp2 ;
+    }
+
+    double dist(std::vector< RVec > alfa, std::vector< RVec > beta, std::vector< std::pair< int, int > > cros) {
+        double distance = 0;
+        for (int i = 0; i < cros.size(); i++) {
+            distance += F(alfa[cros[i].first][0], alfa[cros[i].first][1], alfa[cros[i].first][2], beta[cros[i].second][0], beta[cros[i].second][1], beta[cros[i].second][2], 1, 0, 1, 0, 1, 0);
+        }
+        return distance;
+    }
+
+    void new_fit(std::vector< RVec > alfa, std::vector< RVec > &beta, std::vector< std::pair< int, int > > cros, double prec) {
+        double bx, by, bz;
+        RVec step;
+        double x1 = 0, y1 = 0, z1 = 0, A1 = 0, B1 = 0, C1 = 0;
+        double x2 = 0, y2 = 0, z2 = 0, A2 = 0, B2 = 0, C2 = 0;
+        double cA, sA, cB, sB, cC, sC, temp1, temp2, temp3;
+        double H1 = 0.00001, H2 = 0.00001, F1 = 0, F2 = 0;
+        do {
+            if (F1 > F2 || F1 == 0) {
+                F1 = F2;
+            }
+            F2 = 0;
+            for (int i = 0; i < cros.size(); i++) {
+                bx = beta[cros[i].second][0] + x1;
+                by = beta[cros[i].second][1] + y1;
+                bz = beta[cros[i].second][2] + z1;
+                cA = cos(A1);
+                sA = sin(A1);
+                cB = cos(B1);
+                sB = sin(B1);
+                cC = cos(C1);
+                sC = sin(C1);
+                temp1 = alfa[cros[i].first][0] + by*(cA*sC - cC*sA*sB) - bz*(sA*sC + cA*cC*sB) - cB*cC*bx;
+                temp2 = alfa[cros[i].first][1] - by*(cA*cC + sA*sB*sC) + bz*(cC*sA - cA*sB*sC) - cB*sC*bx;
+                temp3 = alfa[cros[i].first][2] + sB*bx - cA*cB*bz - cB*sA*by;
+                x2 += Fx(temp1, temp2, temp3,                     cB, sB, cC, sC);
+                y2 += Fy(temp1, temp2, temp3,             cA, sA, cB, sB, cC, sC);
+                z2 += Fz(temp1, temp2, temp3,             cA, sA, cB, sB, cC, sC);
+                A2 += FA(temp1, temp2, temp3,     by, bz, cA, sA, cB, sB, cC, sC);
+                B2 += FB(temp1, temp2, temp3, bx, by, bz, cA, sA, cB, sB, cC, sC);
+                C2 += FC(temp1, temp2,        bx, by, bz, cA, sA, cB, sB, cC, sC);
+            }
+            if (std::max(std::fabs(H1*x2), std::max(std::fabs(H1*y2), std::fabs(H1*z2))) > 1) {
+                H1 = 1 / std::max(std::fabs(H1*x2), std::max(std::fabs(H1*y2), std::fabs(H1*z2)));
+            }
+            x2 = x1 - H1 * x2;
+            y2 = y1 - H1 * y2;
+            z2 = z1 - H1 * z2;
+            if (std::max(std::fabs(H2*A2), std::max(std::fabs(H2*B2), std::fabs(H2*C2))) > M_PI/60) {
+                H2 = M_PI/60 / std::max(std::fabs(H2*A2), std::max(std::fabs(H2*B2), std::fabs(H2*C2)));
+            }
+            A2 = A1 - H2 * A2;
+            if (std::fabs(A2) > M_PI) {
+                A2 = 0;
+            }
+            B2 = B1 - H2 * B2;
+            if (std::fabs(B2) > M_PI) {
+                B2 = 0;
+            }
+            C2 = C1 - H2 * C2;
+            if (std::fabs(C2) > M_PI) {
+                C2 = 0;
+            }
+            for (int i = 0; i < cros.size(); i++) {
+                cA = cos(A2);
+                sA = sin(A2);
+                cB = cos(B2);
+                sB = sin(B2);
+                cC = cos(C2);
+                sC = sin(C2);
+                F2 += F(alfa[cros[i].first][0], alfa[cros[i].first][1], alfa[cros[i].first][2], beta[cros[i].second][0] + x2, beta[cros[i].second][1] + y2, beta[cros[i].second][2] + z2, cA, sA, cB, sB, cC, sC);
+            }
+            if (F2 > F1 && F1 != 0) {
+                H1 /= 10;
+                H2 /= 10;
+            } else {
+                x1 = x2;
+                y1 = y2;
+                z1 = z2;
+                A1 = A2;
+                B1 = B2;
+                C1 = C2;
+                x2 = 0;
+                y2 = 0;
+                z2 = 0;
+                A2 = 0;
+                B2 = 0;
+                C2 = 0;
+            }
+        } while (std::fabs(F1 - F2) > prec/*0.00001*/);
+        step[0] = x1;
+        step[1] = y1;
+        step[2] = z1;
+        implement_fit(beta, step, A1, B1, C1);
+    }
+}
+#endif // NEWFIT_H