diff --git a/src/phasicFlow/CMakeLists.txt b/src/phasicFlow/CMakeLists.txt
index 4449aa88..58117785 100644
--- a/src/phasicFlow/CMakeLists.txt
+++ b/src/phasicFlow/CMakeLists.txt
@@ -63,6 +63,7 @@ containers/Field/Fields.cpp
 containers/symArrayHD/symArrays.cpp
 containers/List/anyList/anyList.cpp
 
+structuredData/zAxis/zAxis.cpp
 structuredData/box/box.cpp
 structuredData/sphere/sphere.cpp
 structuredData/cylinder/cylinder.cpp
@@ -73,8 +74,6 @@ structuredData/peakableRegion/peakableRegion/peakableRegion.cpp
 structuredData/pointStructure/internalPoints/internalPointsKernels.cpp
 structuredData/pointStructure/internalPoints/internalPoints.cpp
 
-structuredData/zAxis/zAxis.cpp
-
 structuredData/line/line.cpp
 structuredData/infinitePlane/infinitePlane.cpp
 structuredData/plane/plane.cpp
diff --git a/src/phasicFlow/structuredData/cylinder/cylinder.cpp b/src/phasicFlow/structuredData/cylinder/cylinder.cpp
index b8ce3c68..8c706d3e 100644
--- a/src/phasicFlow/structuredData/cylinder/cylinder.cpp
+++ b/src/phasicFlow/structuredData/cylinder/cylinder.cpp
@@ -27,7 +27,6 @@ FUNCTION_H
 bool pFlow::cylinder::calculateParams()
 {
 
-	WARNING<<"Use of cylinder requires modifications to zAxis"<<END_WARNING;
 	auto p1p2 = p2_ - p1_;
 
 	if( p1p2.length() > smallValue )
@@ -45,9 +44,12 @@ bool pFlow::cylinder::calculateParams()
 	realx3 minPinZ(-sqrt(radius2_), -sqrt(radius2_), 0.0);
 	realx3 maxPinZ( sqrt(radius2_),  sqrt(radius2_), sqrt(axisVector2_));
 
-	minPoint_ = zA.transferBackZ(minPinZ);
-	maxPoint_ = zA.transferBackZ(maxPinZ);
+	auto minp = zA.transferBackZ(minPinZ);
+	auto maxp = zA.transferBackZ(maxPinZ);
 
+	minPoint_ = min(minp, maxp);
+	maxPoint_ = max(minp, maxp);
+	
 	return true;
 }
 
diff --git a/src/phasicFlow/structuredData/zAxis/array2D.hpp b/src/phasicFlow/structuredData/zAxis/array2D.hpp
new file mode 100644
index 00000000..f76e6080
--- /dev/null
+++ b/src/phasicFlow/structuredData/zAxis/array2D.hpp
@@ -0,0 +1,188 @@
+
+#ifndef __array2D_hpp__
+#define __array2D_hpp__
+
+#include "iOstream.hpp"
+
+namespace pFlow
+{
+
+template<typename T, size_t nRow, size_t nCol>
+struct array2D
+{
+
+    T   elements_[nRow][nCol];
+
+    inline
+    T& operator()(size_t i, size_t j) noexcept
+    {
+        return elements_[i][j];
+    }
+
+    inline
+    const T& operator()(size_t i, size_t j)const noexcept
+    {
+        return elements_[i][j];
+    }
+
+    constexpr size_t size()const noexcept
+    {
+        return  nRow*nCol;
+    }
+
+
+    constexpr size_t nCols()const noexcept
+    {
+        return nCols;
+    }
+
+    constexpr size_t nRows()const noexcept
+    {
+        return nRows;
+    }
+    
+};
+
+template<typename T, size_t nRow, size_t nCol>
+array2D<T, nRow, nCol> operator+(
+    const array2D<T, nRow, nCol>& arr1, 
+    const array2D<T, nRow, nCol>& arr2)
+{
+    array2D<T, nRow, nCol> res;
+
+    for(size_t i=0; i<nRow; i++)
+        for(size_t j=0; j<nCol; j++)
+            res(i,j)= arr1(i,j)+arr2(i,j);
+    return res;
+}
+
+template<typename T, size_t nRow, size_t nCol>
+array2D<T, nRow, nCol> operator-(
+    const array2D<T, nRow, nCol>& arr1, 
+    const array2D<T, nRow, nCol>& arr2)
+{
+    array2D<T, nRow, nCol> res;
+
+    for(size_t i=0; i<nRow; i++)
+        for(size_t j=0; j<nCol; j++)
+            res(i,j)= arr1(i,j)-arr2(i,j);
+    return res;
+}
+
+template<typename T, size_t nRow, size_t nCol>
+array2D<T, nRow, nCol> operator*(
+    const array2D<T, nRow, nCol>& arr1, 
+    const array2D<T, nRow, nCol>& arr2)
+{
+    array2D<T, nRow, nCol> res;
+
+    for(size_t i=0; i<nRow; i++)
+        for(size_t j=0; j<nCol; j++)
+            res(i,j)= arr1(i,j)*arr2(i,j);
+    return res;
+}
+
+template<typename T, size_t nRow, size_t nCol>
+array2D<T, nRow, nCol> operator*(
+    const T& s, 
+    const array2D<T, nRow, nCol>& arr2)
+{
+    array2D<T, nRow, nCol> res;
+
+    for(size_t i=0; i<nRow; i++)
+        for(size_t j=0; j<nCol; j++)
+            res(i,j)= s*arr2(i,j);
+    return res;
+}
+
+template<typename T, size_t nRow, size_t nCol>
+array2D<T, nRow, nCol> operator/(
+    const array2D<T, nRow, nCol>& arr1, 
+    const array2D<T, nRow, nCol>& arr2)
+{
+    array2D<T, nRow, nCol> res;
+
+    for(size_t i=0; i<nRow; i++)
+        for(size_t j=0; j<nCol; j++)
+            res(i,j)= arr1(i,j)/arr2(i,j);
+    return res;
+}
+
+template<typename T, size_t nRow, size_t nCol>
+array2D<T, nRow, nCol> operator/(
+    const T& val,
+    const array2D<T, nRow, nCol>& arr2)
+{
+    array2D<T, nRow, nCol> res;
+
+    for(size_t i=0; i<nRow; i++)
+        for(size_t j=0; j<nCol; j++)
+            res(i,j)= val/arr2(i,j);
+    return res;
+}
+
+template <typename T, size_t nRow, size_t nInner, size_t nCol >
+array2D<T, nRow, nCol> 
+MatMul(
+    const array2D<T, nRow, nInner>& A, 
+    const array2D<T, nInner, nCol>& B)
+{
+    array2D<T, nRow, nCol> C;
+	
+    for (size_t row = 0; row < nRow; row++)
+	{
+		for (size_t col = 0; col < nCol; col++)
+		{
+			T sum = 0;
+			for (size_t inner = 0; inner < nInner; inner++)
+			{
+				sum += A(row,inner) * B(inner,col);
+			}
+			C(row,col) = sum;
+		}
+	}
+
+    return C;
+}
+
+template<typename T, size_t nRow, size_t nCol>
+array2D<T, nCol, nRow> 
+transpose(const array2D<T, nRow, nCol>& arr)
+{
+    array2D<T, nCol, nRow> tArr;
+
+    for(size_t i=0; i<nRow; i++)
+    {
+        for(size_t j=0; j<nCol; j++)
+        {
+            tArr(j,i) = arr(i,j);
+        }
+    }
+    return tArr;
+}
+
+template<typename T, size_t nRow, size_t nCol>
+iOstream& operator<<(iOstream& os, const array2D<T, nRow, nCol>& arr)
+{
+    os<<'[';
+    for(size_t i=0; i<nRow; i++)
+    {
+        os<<'[';
+        for(size_t j=0; j<nCol-1; j++)
+        {
+            os<<arr(i,j)<<' ';
+        }
+
+        if(i < nRow-1)    
+            os<<arr(i,nCol-1)<<"]\n";
+        else
+            os<<arr(i,nCol-1)<<"]";
+    }
+    os<<']';
+
+    return os;
+}
+
+}
+
+#endif //__array2D_hpp__
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/zAxis/zAxis.cpp b/src/phasicFlow/structuredData/zAxis/zAxis.cpp
index ee906c8f..81be8bf5 100644
--- a/src/phasicFlow/structuredData/zAxis/zAxis.cpp
+++ b/src/phasicFlow/structuredData/zAxis/zAxis.cpp
@@ -22,129 +22,60 @@ Licence:
 #include "zAxis.hpp"
 
 
-pFlow::zAxis::zAxis(const realx3 &p1, const realx3 &p2) :
-p1_(p1),
-p2_(p2)
+pFlow::zAxis::zAxis(const realx3 &p1, const realx3 &p2) 
+:
+	p1_(p1),
+	p2_(p2)
 {
-	n_ = p2-p1;
-	auto len = pFlow::length(n_);
-
-	if(len < smallValue )
+	if( equal(p1,p2))
 	{
+		fatalErrorInFunction<<
+		"points are equal "<< p1 <<" "<<p2<<endl;
 		fatalExit;
 	}
-	n_ /= len;
-
-	makeTransMatrix();
+	makeTransMatrix(p2-p1);
 }
 
-pFlow::realx3 pFlow::zAxis::transferToZ(const realx3 & p)
+pFlow::realx3 pFlow::zAxis::transferToZ(const realx3 & p)const
 {
-	real pp[4][1] = {p.x(), p.y(), p.z(), 1.0};
-	real pn[4][1];
-
-	MatMul(Trans_z_xz_P1_, pp, pn);
+	const auto p0 = p1_-p;
+	array2D<real,3,1> po {p0.x(), p0.y(), p0.z()}; 
 	
-	return realx3(
-		pn[0][0],
-		pn[1][0],
-		pn[2][0]
-		);
+	auto tp = MatMul(rotMat_, po);
+	
+	return realx3(tp(0,0), tp(1,0), tp(2,0));
 
 }
 
-pFlow::realx3 pFlow::zAxis::transferBackZ(const realx3 & p)
+pFlow::realx3 pFlow::zAxis::transferBackZ(const realx3 & p)const
 {
-	real pp[4][1] = { p.x(), p.y(), p.z(), 1.0 };
-	real pn[4][1];
 
-	MatMul(ITrans_P1_xz_z_, pp, pn);
+	auto rp = MatMul(invRotMat_,array2D<real,3,1>{p.x(), p.y(), p.z()});
 
-	return realx3(
-		pn[0][0],
-		pn[1][0],
-		pn[2][0]
-		);
+	return realx3(rp(0,0)+p1_.x(),rp(1,0)+p1_.y(), rp(2,0)+p1_.z());
 }
 
 
-void pFlow::zAxis::makeTransMatrix()
+
+
+void pFlow::zAxis::makeTransMatrix(const realx3& vec)
 {
-	
-	// transfering point p1 to the origin
-	real TransP1[4][4] =
-	{
-		1.0, 0.0, 0.0, -p1_.x(),
-		0.0, 1.0, 0.0, -p1_.y(),
-		0.0, 0.0, 1.0, -p1_.z(),
-		0.0, 0.0, 0.0, 1.0
-	};
+	const realx3 Unit = {0,0,1};
+	const auto v_unit = normalize(vec);
+	auto uvw = cross(v_unit,Unit);
 
-	// for transformation back
-	real ITransP1[4][4] =
-	{
-		1.0, 0.0, 0.0, p1_.x(),
-		0.0, 1.0, 0.0, p1_.y(),
-		0.0, 0.0, 1.0, p1_.z(),
-		0.0, 0.0, 0.0, 1.0
-	};
+	const auto rcos = dot(v_unit, Unit);
+    const auto rsin = uvw.length();
 
-	
-	real u = n_.x();
-	real v = n_.y();
-	real w = n_.z();
-	real u2v2 = sqrt(u*u + v*v);
+	if (abs(rsin) > smallValue)
+        uvw /= rsin;
+	const auto& [u, v, w] = uvw;
 
-	//correcting the transformation matrix in the case of coincidence with z - axis
-	if ( equal(w,1.0) )
-	{
-		assignMat(TransP1 , Trans_z_xz_P1_);
-		assignMat(ITransP1, ITrans_P1_xz_z_);
-		return;
-	}
+	rotMat_ = 
+	ArrayType({rcos, 0, 0, 0, rcos, 0, 0, 0, rcos})+
+	ArrayType({0, -w, v, w, 0, -u, -v, u, 0})+
+	(1-rcos)*ArrayType({u*u, u*v, u*w, u*v, v*v, w*v, u*w, v*w, w*w});
 
-	u2v2 = max(smallValue, u2v2);
+	invRotMat_ = transpose(rotMat_);
 
-	real TransXZ[4][4] =
-	{
-		u / u2v2, v / u2v2, 0.0, 0.0,
-		-v / u2v2, u / u2v2, 0.0, 0.0,
-		0.0, 0.0, 1.0, 0.0,
-		0.0, 0.0, 0.0, 1.0
-	};
-
-	real TransZ[4][4] =
-	{
-		w, 0.0, -u2v2, 0.0,
-		0.0, 1.0, 0.0, 0.0,
-		u2v2, 0.0, w, 0.0,
-		0.0, 0.0, 0.0, 1.0
-	};
-
-	real temp[4][4];
-	// creat transformation matrix to transfer point from line axis to z-axis
-	MatMul(TransXZ, TransP1, temp);
-	MatMul(TransZ, temp, Trans_z_xz_P1_);
-
-
-	real ITransXZ[4][4] =
-	{
-		u / u2v2, -v / u2v2, 0.0, 0.0,
-		+v / u2v2, u / u2v2, 0.0, 0.0,
-		0.0, 0.0, 1.0, 0.0,
-		0.0, 0.0, 0.0, 1.0
-	};
-
-	real ITransZ[4][4] =
-	{
-		w, 0.0, +u2v2, 0.0,
-		0.0, 1.0, 0.0, 0.0,
-		-u2v2, 0.0, w, 0.0,
-		0.0, 0.0, 0.0, 1.0
-	};
-
-	// creat transformation matrix to transfer point to from z-axis to line axis
-	MatMul(ITransXZ, ITransZ, temp);
-	MatMul(ITransP1, temp, ITrans_P1_xz_z_);
-
-}
+}  
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/zAxis/zAxis.hpp b/src/phasicFlow/structuredData/zAxis/zAxis.hpp
index 7dabf7cc..b791164b 100644
--- a/src/phasicFlow/structuredData/zAxis/zAxis.hpp
+++ b/src/phasicFlow/structuredData/zAxis/zAxis.hpp
@@ -29,78 +29,53 @@ p1 and p2
 #define __zAxis_hpp__
 
 #include "types.hpp"
+#include "array2D.hpp"
 
 namespace pFlow
 {
 
-template <typename T, int32 nRow, int32 nInner, int32 nCol >
-void MatMul(T(&A)[nRow][nInner], T(&B)[nInner][nCol], T(&C)[nRow][nCol]);
-
-template <typename T, int32 nRow, int32 nCol >
-void assignMat(T(&A)[nRow][nCol], T(&B)[nRow][nCol]);
 
 class zAxis
 {
+
+private:
+	
+	using ArrayType = array2D<real,3uL,3uL>;
+
+	/// the origin of the transformed coord
+	realx3 			p1_;
+
+	/// the direction vector of rotated coordinates 
+	realx3 			p2_;
+
+	/// transformation matrix to rotate original coordinates
+	/// to rotated coordinates 
+	ArrayType 		rotMat_;
+
+	/// rotation matrix to rotate back from rotated coordinates
+	/// to the original axis
+	ArrayType 		invRotMat_;
+
+	void makeTransMatrix(const realx3& v);
+
 public:
 	// constructors 
 	zAxis(const realx3 &lp1, const realx3 &lp2);
 	
+	inline
 	real length()const
 	{
 		return pFlow::length(p2_-p1_);
 	}
 
-	realx3 transferToZ(const realx3 & p);
+	realx3 transferToZ(const realx3 & p)const;
 
-	realx3 transferBackZ(const realx3 & p);
+	realx3 transferBackZ(const realx3 & p)const;
 
-private:
-	void makeTransMatrix();
-protected:
 
-	realx3 p1_;
-	realx3 p2_;
-	realx3 n_;
-
-	real Trans_z_xz_P1_[4][4];
-
-	real ITrans_P1_xz_z_[4][4];
 };
 
 
-
-template <typename T, int32 nRow, int32 nInner, int32 nCol >
-void MatMul(T(&A)[nRow][nInner], T(&B)[nInner][nCol], T(&C)[nRow][nCol])
-{
-
-	for (int32 row = 0; row < nRow; row++)
-	{
-		for (int32 col = 0; col < nCol; col++)
-		{
-			T sum = 0;
-			for (int inner = 0; inner < nInner; inner++)
-			{
-				sum += A[row][inner] * B[inner][col];
-			}
-			C[row][col] = sum;
-		}
-	}
-}
-
-template <typename T, int32 nRow, int32 nCol >
-void assignMat(T(&A)[nRow][nCol], T(&B)[nRow][nCol])
-{
-
-	for (int32 row = 0; row < nRow; row++)
-	{
-		for (int32 col = 0; col < nCol; col++)
-		{
-			B[row][col] = A[row][col];
-		}
-	}
-}
-
-
 }
 
 #endif