Changeset 20895b for molecuilder/src/helpers.cpp

Timestamp:

Nov 4, 2009, 5:34:05 PM (16 years ago)

Author:

Frederik Heber <heber@…>

Children:

5f1d021

Parents:

c1b76e

Message:

Extension to the periodic boundary case for analysis_correlation.cpp

• Distance() -> PeriodicDistance(): PairCorrelation(), CorrelationToPoint(), CorrelationToSurface()
• CorrelationToSurface() - goes through periodic cells given by a ranges variable and puts into map as well.
• const: ReturnFullMatrixforSymmetric() and InverseMatrix() now have const * const parameter
• InverseMatrix() shifted from class Vector into helpers.cpp
• adapted due to InverseMatrix(): molecule::CenterInBox(), molecule::BoundInBox(), molecule::TranslatePeriodically()
• new function: RDET3(), RDET2() (former #define's)
• FIX: molecule::CorrectBondDegree() would loop eternally, if bonds cannot be fully corrected.
• FIX for NULL return: Tesselation::FindClosestTrianglesToPoint(), Tesselation::GetCircleOfConnectedPoints()
• FIX: Vector::DistanceToPlane() - did not give correct sign before according to PlaneNormal orientation

other stuff:

• less output: BondedParticle::CorrectBondDegree(), molecule::CyclicStructureAnalysis()

File:

• molecuilder/src/helpers.cpp (modified) (2 diffs)

molecuilder/src/helpers.cpp

@@ -133 +133 @@
  * \return allocated NDIM*NDIM array with the symmetric matrix
  */
-double * ReturnFullMatrixforSymmetric(double *symm)
+double * ReturnFullMatrixforSymmetric(const double * const symm)
 {
   double *matrix = Malloc<double>(NDIM * NDIM, "molecule::ReturnFullMatrixforSymmetric: *matrix");
@@ -148 +148 @@
 };
 
+/** Calculate the inverse of a 3x3 matrix.
+ * \param *matrix NDIM_NDIM array
+ */
+double * InverseMatrix( const double * const A)
+{
+  double *B = Malloc<double>(NDIM * NDIM, "Vector::InverseMatrix: *B");
+  double detA = RDET3(A);
+  double detAReci;
+
+  for (int i=0;i<NDIM*NDIM;++i)
+    B[i] = 0.;
+  // calculate the inverse B
+  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
+    detAReci = 1./detA;
+    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
+    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
+    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
+    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
+    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
+    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
+    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
+    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
+    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
+  }
+  return B;
+};
+
+/** hard-coded determinant of a 3x3 matrix.
+ * \param a[9] matrix
+ * \return \f$det(a)\f$
+ */
+double RDET3(const double a[NDIM*NDIM])
+{
+  return ((a)[0]*(a)[4]*(a)[8] + (a)[3]*(a)[7]*(a)[2] + (a)[6]*(a)[1]*(a)[5] - (a)[2]*(a)[4]*(a)[6] - (a)[5]*(a)[7]*(a)[0] - (a)[8]*(a)[1]*(a)[3]);
+};
+
+/** hard-coded determinant of a 2x2 matrix.
+ * \param a[4] matrix
+ * \return \f$det(a)\f$
+ */
+double RDET2(const double a[4])
+{
+  return ((a[0])*(a[3])-(a[1])*(a[2]));
+};
+
+/** hard-coded determinant of a 2x2 matrix.
+ * \param a0 (0,0) entry of matrix
+ * \param a1 (0,1) entry of matrix
+ * \param a2 (1,0) entry of matrix
+ * \param a3 (1,1) entry of matrix
+ * \return \f$det(a)\f$
+ */
+double RDET2(const double a0, const double a1, const double a2, const double a3)
+{
+  return ((a0)*(a3)-(a1)*(a2));
+};
+
 /** Comparison function for GSL heapsort on distances in two molecules.
  * \param *a