diff --git a/GeoLib/AnalyticalGeometry-impl.h b/GeoLib/AnalyticalGeometry-impl.h
index e2ad36a5ec901f02afb885e1277b52a1a492584b..27adefe1b327b39ae930cc20ceb0480a021317a8 100644
--- a/GeoLib/AnalyticalGeometry-impl.h
+++ b/GeoLib/AnalyticalGeometry-impl.h
@@ -78,5 +78,63 @@ void compute3DRotationMatrixToX(MathLib::Vector3  const& v,
     }
 }
 
+template <class T_MATRIX>
+void computeRotationMatrixToXY(MathLib::Vector3 const& n,
+		T_MATRIX & rot_mat)
+{
+	// check if normal points already in the right direction
+	if (sqrt(n[0]*n[0]+n[1]*n[1]) == 0) {
+		rot_mat(0,0) = 1.0;
+		rot_mat(0,1) = 0.0;
+		rot_mat(0,2) = 0.0;
+		rot_mat(1,0) = 0.0;
+		rot_mat(1,1) = 1.0;
+		rot_mat(1,2) = 0.0;
+		rot_mat(2,0) = 0.0;
+		rot_mat(2,1) = 0.0;
+		rot_mat(2,2) = 1.0;
+		return;
+	}
+
+	// sqrt (n_1^2 + n_3^2)
+	double const h0(sqrt(n[0]*n[0]+n[2]*n[2]));
+
+	// In case the x and z components of the normal are both zero the rotation
+	// to the x-z-plane is not required, i.e. only the rotation in the z-axis is
+	// required. The angle is either pi/2 or 3/2*pi. Thus the components of
+	// rot_mat are as follows.
+	if (h0 < std::numeric_limits<double>::epsilon()) {
+		rot_mat(0,0) = 1.0;
+		rot_mat(0,1) = 0.0;
+		rot_mat(0,2) = 0.0;
+		rot_mat(1,0) = 0.0;
+		rot_mat(1,1) = 0.0;
+		if (n[1] > 0)
+			rot_mat(1,2) = -1.0;
+		else
+			rot_mat(1,2) = 1.0;
+		rot_mat(2,0) = 0.0;
+		if (n[1] > 0)
+			rot_mat(2,1) = 1.0;
+		else
+			rot_mat(2,1) = -1.0;
+		rot_mat(2,2) = 0.0;
+		return;
+	}
+
+	double h1(1 / n.getLength());
+
+	// general case: calculate entries of rotation matrix
+	rot_mat(0, 0) = n[2] / h0;
+	rot_mat(0, 1) = 0;
+	rot_mat(0, 2) = - n[0] / h0;
+	rot_mat(1, 0) = - n[1]*n[0]/h0 * h1;
+	rot_mat(1, 1) = h0 * h1;
+	rot_mat(1, 2) = - n[1]*n[2]/h0 * h1;
+	rot_mat(2, 0) = n[0] * h1;
+	rot_mat(2, 1) = n[1] * h1;
+	rot_mat(2, 2) = n[2] * h1;
+}
+
 } // end namespace GeoLib
 
diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index bd9900cdd2e7c9983d0219a653172b873d471d71..67b5f1d6e0f4f1c3e246110d56efd4fe258203dd 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -348,63 +348,6 @@ double calcTetrahedronVolume(MathLib::Point3d const& x1,
 	return std::abs(GeoLib::scalarTriple(ac, ad, ab)) / 6.0;
 }
 
-void computeRotationMatrixToXY(MathLib::Vector3 const& n,
-	MathLib::DenseMatrix<double> & rot_mat)
-{
-	// check if normal points already in the right direction
-	if (sqrt(n[0]*n[0]+n[1]*n[1]) == 0) {
-		rot_mat(0,0) = 1.0;
-		rot_mat(0,1) = 0.0;
-		rot_mat(0,2) = 0.0;
-		rot_mat(1,0) = 0.0;
-		rot_mat(1,1) = 1.0;
-		rot_mat(1,2) = 0.0;
-		rot_mat(2,0) = 0.0;
-		rot_mat(2,1) = 0.0;
-		rot_mat(2,2) = 1.0;
-		return;
-	}
-
-	// sqrt (n_1^2 + n_3^2)
-	double const h0(sqrt(n[0]*n[0]+n[2]*n[2]));
-
-	// In case the x and z components of the normal are both zero the rotation
-	// to the x-z-plane is not required, i.e. only the rotation in the z-axis is
-	// required. The angle is either pi/2 or 3/2*pi. Thus the components of
-	// rot_mat are as follows.
-	if (h0 < std::numeric_limits<double>::epsilon()) {
-		rot_mat(0,0) = 1.0;
-		rot_mat(0,1) = 0.0;
-		rot_mat(0,2) = 0.0;
-		rot_mat(1,0) = 0.0;
-		rot_mat(1,1) = 0.0;
-		if (n[1] > 0)
-			rot_mat(1,2) = -1.0;
-		else
-			rot_mat(1,2) = 1.0;
-		rot_mat(2,0) = 0.0;
-		if (n[1] > 0)
-			rot_mat(2,1) = 1.0;
-		else
-			rot_mat(2,1) = -1.0;
-		rot_mat(2,2) = 0.0;
-		return;
-	}
-
-	double h1(1 / n.getLength());
-
-	// general case: calculate entries of rotation matrix
-	rot_mat(0, 0) = n[2] / h0;
-	rot_mat(0, 1) = 0;
-	rot_mat(0, 2) = - n[0] / h0;
-	rot_mat(1, 0) = - n[1]*n[0]/h0 * h1;
-	rot_mat(1, 1) = h0 * h1;
-	rot_mat(1, 2) = - n[1]*n[2]/h0 * h1;
-	rot_mat(2, 0) = n[0] * h1;
-	rot_mat(2, 1) = n[1] * h1;
-	rot_mat(2, 2) = n[2] * h1;
-}
-
 void computeRotationMatrixToXZ(MathLib::Vector3 const& plane_normal, MathLib::DenseMatrix<double> & rot_mat)
 {
 	// *** some frequently used terms ***
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index 3485bf32617038662998f8f97fb15d98f508b09c..aea8be332bbd01a25bf2c57af53b642a34db1210 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -86,8 +86,9 @@ void compute3DRotationMatrixToX(MathLib::Vector3  const& v, T_MATRIX & rot_mat);
  * @param plane_normal the (3d) vector that is rotated parallel to the \f$z\f$ axis
  * @param rot_mat 3x3 rotation matrix
  */
+template <class T_MATRIX>
 void computeRotationMatrixToXY(MathLib::Vector3 const& plane_normal,
-                               MathLib::DenseMatrix<double> & rot_mat);
+                               T_MATRIX & rot_mat);
 
 /**
  * Method computes the rotation matrix that rotates the given vector parallel to the \f$y\f$ axis.