diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp
index 8314cb6..3c45a16 100644
--- a/FDTD/operator.cpp
+++ b/FDTD/operator.cpp
@@ -102,6 +102,14 @@ unsigned int Operator::CalcNyquistNum(double fmax)
 	return floor(T0/2/dT);
 }
 
+string Operator::GetDirName(int ny) const
+{
+	if (ny==0) return "x";
+	if (ny==1) return "y";
+	if (ny==2) return "z";
+	return "";
+}
+
 double Operator::GetMeshDelta(int n, const unsigned int* pos, bool dualMesh) const
 {
 	if ((n<0) || (n>2)) return 0.0;
diff --git a/FDTD/operator.h b/FDTD/operator.h
index ee1a9b4..5a7aa8d 100644
--- a/FDTD/operator.h
+++ b/FDTD/operator.h
@@ -67,6 +67,9 @@ public:
 	void DumpOperator2File(string filename);
 	void DumpMaterial2File(string filename);
 
+	//! Get the name for the given direction: 0 -> x, 1 -> y, 2 -> z
+	virtual string GetDirName(int ny) const;
+
 	virtual double GetGridDelta() const {return gridDelta;}
 	//! Get the mesh delta times the grid delta for a 3D position
 	virtual double GetMeshDelta(int n, const int* pos, bool dualMesh=false) const;
diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp
index e67e5d2..0f0f9bb 100644
--- a/FDTD/operator_cylinder.cpp
+++ b/FDTD/operator_cylinder.cpp
@@ -68,6 +68,14 @@ inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny) const
 	return numLines[ny];
 }
 
+string Operator_Cylinder::GetDirName(int ny) const
+{
+	if (ny==0) return "rho";
+	if (ny==1) return "alpha";
+	if (ny==2) return "z";
+	return "";
+}
+
 double Operator_Cylinder::GetMeshDelta(int n, int* pos, bool dualMesh) const
 {
 	double delta = Operator::GetMeshDelta(n,pos,dualMesh);
diff --git a/FDTD/operator_cylinder.h b/FDTD/operator_cylinder.h
index 62fa279..df90261 100644
--- a/FDTD/operator_cylinder.h
+++ b/FDTD/operator_cylinder.h
@@ -35,6 +35,9 @@ public:
 
 	virtual unsigned int GetNumberOfLines(int ny) const;
 
+	//! Get the name for the given direction: 0 -> rho, 1 -> alpha, 2 -> z
+	virtual string GetDirName(int ny) const;
+
 	//! Get the mesh delta times the grid delta for a 3D position, including radius corrected alpha-mesh width
 	virtual double GetMeshDelta(int n, int* pos, bool dualMesh=false) const;
 
diff --git a/FDTD/processfields.cpp b/FDTD/processfields.cpp
index 75a8257..0ef6cb6 100644
--- a/FDTD/processfields.cpp
+++ b/FDTD/processfields.cpp
@@ -51,7 +51,8 @@ ProcessFields::~ProcessFields()
 void ProcessFields::InitProcess()
 {
 	if (Enabled==false) return;
-	string names[] = {"x","y","z"};
+	//get the correct direction names for all coordinate systems
+	string names[] = {Op->GetDirName(0),Op->GetDirName(1),Op->GetDirName(2)};
 	if (m_fileType==HDF5_FILETYPE)
 	{
 		unsigned int* NrLines;
diff --git a/matlab/PlotHDF5FieldData.m b/matlab/PlotHDF5FieldData.m
index 0c90613..b1b5134 100644
--- a/matlab/PlotHDF5FieldData.m
+++ b/matlab/PlotHDF5FieldData.m
@@ -11,7 +11,11 @@ end
 mesh = ReadHDF5Mesh(file);
 fields = ReadHDF5FieldData(file);
 
-[X Y Z] = meshgrid(double(mesh.lines{1}),double(mesh.lines{2}),double(mesh.lines{3}));
+if (mesh.type==0)
+    [X Y Z] = meshgrid(double(mesh.lines{1}),double(mesh.lines{2}),double(mesh.lines{3}));
+else
+    disp(['PlotHDF5FieldData:: Error: unknown mesh type ' num2str(mesh.type)]);
+end
 
 max_amp = 0;
 
diff --git a/matlab/ReadHDF5Mesh.m b/matlab/ReadHDF5Mesh.m
index 472b94b..5988e90 100644
--- a/matlab/ReadHDF5Mesh.m
+++ b/matlab/ReadHDF5Mesh.m
@@ -13,4 +13,10 @@ end
 hdf_mesh.names = names;
 for n=1:numel(names)
     hdf_mesh.lines{n} = hdf5read(file,names{n});
+end
+
+if (strcmp(names{1},'/mesh/rho'))
+    hdf_mesh.type=1;
+else
+    hdf_mesh.type=0;
 end
\ No newline at end of file