Unit test comparing old and new implementation of the GeoXF Pow function

GeoModelCore/GeoModelKernel/src/GeoXF.cxx

@@ -232,6 +232,14 @@ __attribute__ ((flatten))
     // Get the translation part and the rotation part:
     //
     GeoTrf::RotationMatrix3D rotate = m_xf.rotation ();
+      
+    //NEW - start
+      Eigen::Matrix3d linear = m_xf.linear ();
+      Eigen::EigenSolver<GeoTrf::RotationMatrix3D> solver(linear);
+      Eigen::MatrixXcd D = solver.eigenvalues().asDiagonal();
+      Eigen::MatrixXcd V = solver.eigenvectors();
+    //NEW - end
+      
     GeoTrf::Vector3D translate = m_xf.translation ();
     Eigen::AngleAxis<double> aa(rotate);
     //
@@ -244,10 +252,67 @@ __attribute__ ((flatten))
     translate *= nTimes;
     double& delta = aa.angle();
     delta *= nTimes;
+      
+      //NEW - start
+      Eigen::Matrix3cd DPowN=Eigen::Matrix3cd::Zero();
+      for (unsigned int i=0;i<3;i++) DPowN(i,i)=pow(D(i,i),nTimes);
+      //NEW - end
+      
     //
     // Now compose these and return a result:
     //
+      
+      //NEW - start
+      GeoTrf::Transform3D tRPowN=GeoTrf::Transform3D::Identity();
+      tRPowN.linear()=(V*DPowN*V.inverse()).real();
+      //return GeoTrf::Translation3D (translate) * tRPowN;
+      //NEW - end
+      
+      GeoTrf::Transform3D newMatrix=GeoTrf::Translation3D (translate) * tRPowN;
+      
+      GeoTrf::Transform3D oldMatrix=GeoTrf::Translation3D (translate) * GeoTrf::Transform3D(aa);
+      std::cout.precision(17);
+
+//      if(
+//         std::abs(newMatrix(0,0)-oldMatrix(0,0))>1.e-15||
+//         std::abs(newMatrix(0,1)-oldMatrix(0,1))>1.e-15||
+//         std::abs(newMatrix(0,2)-oldMatrix(0,2))>1.e-15||
+//         std::abs(newMatrix(1,0)-oldMatrix(1,0))>1.e-15||
+//         std::abs(newMatrix(1,1)-oldMatrix(1,1))>1.e-15||
+//         std::abs(newMatrix(1,2)-oldMatrix(1,2))>1.e-15||
+//         std::abs(newMatrix(2,0)-oldMatrix(2,0))>1.e-15||
+//         std::abs(newMatrix(2,1)-oldMatrix(2,1))>1.e-15||
+//         std::abs(newMatrix(2,2)-oldMatrix(2,2))>1.e-15
+//         )
+          if(
+             std::abs(newMatrix(0,0)-oldMatrix(0,0))>0||
+             std::abs(newMatrix(0,1)-oldMatrix(0,1))>0||
+             std::abs(newMatrix(0,2)-oldMatrix(0,2))>0||
+             std::abs(newMatrix(1,0)-oldMatrix(1,0))>0||
+             std::abs(newMatrix(1,1)-oldMatrix(1,1))>0||
+             std::abs(newMatrix(1,2)-oldMatrix(1,2))>0||
+             std::abs(newMatrix(2,0)-oldMatrix(2,0))>0||
+             std::abs(newMatrix(2,1)-oldMatrix(2,1))>0||
+             std::abs(newMatrix(2,2)-oldMatrix(2,2))>0
+             )
+      {
+         
+          std::cout<<"Old matrix: \n"<<oldMatrix.matrix()<<std::endl;
+          std::cout<<"New matrix: \n"<<newMatrix.matrix()<<std::endl;
+          std::cout<<"Matrices Diff:\n"<<newMatrix.matrix()-oldMatrix.matrix()<<std::endl;
+          std::cout<<"------\n"<<std::endl;
+          //std::cout<<"Old matrix rotation: \n"<<oldMatrix.rotation()<<std::endl;
+          //std::cout<<"New matrix rotation: \n"<<newMatrix.rotation()<<std::endl;
+          //std::cout<<"Old matrix translation: \n"<<oldMatrix.translation()<<std::endl;
+          //std::cout<<"New matrix translation: \n"<<newMatrix.translation()<<std::endl;
+          exit(-1);
+          
+      }
+
+    //OLD VALUE
     return GeoTrf::Translation3D (translate) * GeoTrf::Transform3D(aa);
+    
+    
   }
 
   GeoTrf::Transform3D Pow::operator        () (const GeoGenfun::Argument & argument) const