diff --git a/GeoModelCore/GeoModelKernel/src/GeoXF.cxx b/GeoModelCore/GeoModelKernel/src/GeoXF.cxx index 23cac3f8806844d230ec637264e52b1bfcb85f0b..43dad283bb0f7818f3cb022fd60cffd5d690dffb 100755 --- a/GeoModelCore/GeoModelKernel/src/GeoXF.cxx +++ b/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