23 files changed, 4592 insertions, 0 deletions
diff --git a/openEMS/tools/AdrOp.cpp b/openEMS/tools/AdrOp.cpp
new file mode 100644
index 0000000..228f8de
--- /dev/null
+++ b/openEMS/tools/AdrOp.cpp
@@ -0,0 +1,554 @@
+/*
+*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "AdrOp.h"
+
+using namespace std;
+
+AdrOp::AdrOp(unsigned int muiImax, unsigned int muiJmax, unsigned int muiKmax, unsigned int muiLmax)
+{
+	//error-handling...
+	error = new ErrorMsg(9);
+	if (error==NULL)
+	{
+		fprintf(stderr,"Memory allocation failed!! exiting...");
+		exit(1);
+	}
+	error->SetMsg(1,"Adress Operator: Memory allocation failed!! exiting...");
+	error->SetMsg(2,"Adress Operator: Invalid Adress requested!! exiting...");
+	error->SetMsg(3,"Adress Operator: Invalid Position set!! exiting...");
+	error->SetMsg(4,"Adress Operator: Invalid jump or passing end of iteration!! exiting...");
+	error->SetMsg(5,"Adress Operator: 4D not yet implemented!! exiting...");
+	error->SetMsg(6,"Adress Operator: Position not set!! exiting...");
+	error->SetMsg(7,"Adress Operator: Cells not added to Adress Operator!! exiting...");
+	error->SetMsg(8,"Adress Operator: Invalid Node!! exiting...");
+	error->SetMsg(9,"Adress Operator: Grid invalid!! exiting...");
+
+	//if (muiImax<0) muiImax=0;
+	//if (muiJmax<0) muiJmax=0;
+	//if (muiKmax<0) muiKmax=0;
+	//if (muiLmax<0) muiLmax=0;
+
+	uiDimension=0;
+	if (muiImax>0) uiDimension++;
+	else exit(-1);
+	if (muiJmax>0) uiDimension++;
+	else exit(-2);
+	if (muiKmax>0) uiDimension++;
+	if ( (muiLmax>0) && (muiKmax>0) ) uiDimension++;
+//	cout << "\n-----Adress Operator created: Dimension: " << uiDimension << "----" <<endl;
+	uiImax=muiImax;
+	uiJmax=muiJmax;
+	uiKmax=muiKmax;
+	uiLmax=muiLmax=0;
+	if (uiDimension==2) uiSize=uiImax*uiJmax;
+	else if (uiDimension==3) uiSize=uiImax*uiJmax*uiKmax;
+	else if (uiDimension==4) uiSize=uiImax*uiJmax*uiKmax*uiLmax;
+	else uiSize=0;
+	bPosSet=false;
+	if (uiDimension==4) error->Error(5);
+	iIshift=iJshift=iKshift=0;
+	reflect=false;
+	uiTypeOffset=0;
+	clCellAdr=NULL;
+	dGrid[0]=NULL;
+	dGrid[1]=NULL;
+	dGrid[2]=NULL;
+	dGrid[3]=NULL;
+	dDeltaUnit=1;
+	bDebug=false;
+}
+
+AdrOp::AdrOp(AdrOp* origOP)
+{
+	clCellAdr=NULL;
+	error=NULL; // has to be done!!!
+
+	uiDimension=origOP->uiDimension;
+	uiSize=origOP->uiSize;
+	uiImax=origOP->uiImax;
+	uiJmax=origOP->uiJmax;
+	uiKmax=origOP->uiKmax;
+	uiLmax=origOP->uiLmax;
+	uiIpos=origOP->uiIpos;
+	uiJpos=origOP->uiJpos;
+	uiKpos=origOP->uiKpos;
+	uiLpos=origOP->uiLpos;
+	for (int ii=0; ii<4; ++ii) dGrid[ii]=origOP->dGrid[ii];
+	dDeltaUnit=origOP->dDeltaUnit;
+	iIshift=origOP->iIshift;
+	iJshift=origOP->iJshift;
+	iKshift=origOP->iKshift;
+	for (int ii=0; ii<3; ++ii) iCellShift[ii]=origOP->iCellShift[ii];
+	i=origOP->i;
+	j=origOP->j;
+	k=origOP->k;
+	l=origOP->l;
+	reflect=origOP->reflect;
+	uiTypeOffset=origOP->uiTypeOffset;
+
+	bPosSet=origOP->bPosSet;
+	bDebug=origOP->bDebug;
+//	return;
+	if (origOP->clCellAdr!=NULL) clCellAdr= new AdrOp(origOP->clCellAdr);
+}
+
+AdrOp::~AdrOp()
+{
+//	cerr << "\n------Adress Operator deconstructed-----\n" << endl;
+	delete error;
+	error=NULL;
+	delete clCellAdr;
+	clCellAdr=NULL;
+}
+
+unsigned int AdrOp::SetPos(unsigned int muiIpos, unsigned int muiJpos, unsigned int muiKpos, unsigned int muiLpos)
+{
+	if (bDebug) fprintf(stderr,"AdrOp Debug:: SetPos(%d,%d,%d,%d) Max(%d,%d,%d,%d) \n",muiIpos,muiJpos,muiKpos,muiLpos,uiImax,uiJmax,uiKmax,uiLmax);
+	bPosSet=false;
+	if (muiIpos<uiImax) uiIpos=muiIpos;
+	else error->Error(3);
+	if (muiJpos<uiJmax) uiJpos=muiJpos;
+	else error->Error(3);
+	if ((muiKpos>=uiKmax) && (uiDimension>2)) error->Error(3);
+	else if (uiDimension>2) uiKpos=muiKpos;
+	if ((muiLpos>=uiLmax) && (uiDimension>3)) error->Error(3);
+	else if (uiDimension>3) uiLpos=muiLpos;
+	bPosSet=true;
+//	cerr << "Position i:" << uiIpos << " j: " << uiJpos << " k: " << uiKpos << " l: " << 0 << " MAX: i:" << uiImax << " j: "  << uiJmax << " k: " << uiKmax << endl; //debug
+	ADRESSEXPENSE(0,0,0,0,uiDimension+1,18)
+	return GetPos();
+}
+
+bool  AdrOp::SetPosChecked(unsigned int muiIpos, unsigned int muiJpos, unsigned int muiKpos, unsigned int muiLpos)
+{
+	bPosSet=true;
+	if (muiIpos<uiImax) uiIpos=muiIpos;
+	else bPosSet=false;
+	if (muiJpos<uiJmax) uiJpos=muiJpos;
+	else bPosSet=false;
+	if ((muiKpos>=uiKmax) && (uiDimension>2)) bPosSet=false;
+	else if (uiDimension>2) uiKpos=muiKpos;
+	if ((muiLpos>=uiLmax) && (uiDimension>3)) bPosSet=false;
+	else if (uiDimension>3) uiLpos=muiLpos;
+	ADRESSEXPENSE(0,0,0,0,uiDimension+1,18)
+	return bPosSet;
+}
+
+void AdrOp::SetGrid(double *gridI,double *gridJ,double *gridK,double *gridL)
+{
+	dGrid[0]=gridI;
+	dGrid[1]=gridJ;
+	dGrid[2]=gridK;
+	dGrid[3]=gridL;
+	ADRESSEXPENSE(0,0,0,0,4,0)
+}
+
+bool AdrOp::CheckPos(unsigned int muiIpos, unsigned int muiJpos, unsigned int muiKpos, unsigned int muiLpos)
+{
+	bPosSet=true;
+	if ((muiIpos>=uiImax)) bPosSet=false;
+	if ((muiJpos>=uiJmax)) bPosSet=false;
+	if ((muiKpos>=uiKmax) && (uiDimension>2)) bPosSet=false;
+	if ((muiLpos>=uiLmax) && (uiDimension>3)) bPosSet=false;
+	ADRESSEXPENSE(0,0,0,0,uiDimension+1,18)
+	return bPosSet;
+}
+
+bool AdrOp::CheckRelativePos(int muiIrel,int muiJrel,int muiKrel, int muiLrel)
+{
+	bPosSet=true;
+	if ((muiIrel+(int)uiIpos<0) || (muiIrel+(int)uiIpos>=(int)uiImax)) bPosSet=false;
+	if ((muiJrel+(int)uiJpos<0) || (muiJrel+(int)uiJpos>=(int)uiJmax)) bPosSet=false;
+	if (((muiKrel+(int)uiKpos<0) || (muiKrel+(int)uiKpos>=(int)uiKmax)) && (uiDimension>2)) bPosSet=false;
+	if (((muiLrel+(int)uiLpos<0) || (muiLrel+(int)uiLpos>=(int)uiLmax)) && (uiDimension>3)) bPosSet=false;
+	ADRESSEXPENSE(2*uiDimension,0,0,0,uiDimension+1,18)
+	return bPosSet;
+}
+
+unsigned int AdrOp::GetPos(int muiIrel, int muiJrel, int muiKrel, int /*muiLrel*/)
+{
+	if (bPosSet==false) error->Error(6);
+	if (reflect)
+	{
+#if EXPENSE_LOG==1
+		if (muiIrel+(int)uiIpos<0) ADRESSEXPENSE(2,1,0,0,1,0)
+			if (muiIrel+(int)uiIpos>(int)uiImax-1) ADRESSEXPENSE(4,1,0,0,1,0)
+				if (muiJrel+(int)uiJpos<0) ADRESSEXPENSE(2,1,0,0,1,0)
+					if (muiJrel+(int)uiJpos>(int)uiJmax-1) ADRESSEXPENSE(4,1,0,0,1,0)
+						if (muiKrel+(int)uiKpos<0) ADRESSEXPENSE(2,1,0,0,1,0)
+							if (muiKrel+(int)uiKpos>(int)uiKmax-1) ADRESSEXPENSE(4,1,0,0,1,0)
+#endif
+
+								if (muiIrel+(int)uiIpos<0) muiIrel=-2*uiIpos-muiIrel-uiTypeOffset;
+		if (muiIrel+(int)uiIpos>(int)uiImax-1) muiIrel=2*(uiImax-1-uiIpos)-muiIrel+uiTypeOffset;
+		if (muiJrel+(int)uiJpos<0) muiJrel=-2*uiJpos-muiJrel-uiTypeOffset;
+		if (muiJrel+(int)uiJpos>(int)uiJmax-1) muiJrel=2*(uiJmax-1-uiJpos)-muiJrel+uiTypeOffset;
+		if (muiKrel+(int)uiKpos<0) muiKrel=-2*uiKpos-muiKrel-uiTypeOffset;
+		if (muiKrel+(int)uiKpos>(int)uiKmax-1) muiKrel=2*(uiKmax-1-uiKpos)-muiKrel+uiTypeOffset;
+	}
+	if (uiDimension==2)
+	{
+		ADRESSEXPENSE(7,1,0,0,0,7)
+		if ( (muiIrel+uiIpos<uiImax) && (muiJrel+uiJpos<uiJmax) )
+			return (muiIrel+uiIpos)+(muiJrel+uiJpos)*uiImax;
+		else error->Error(2);
+		return 0;
+	}
+	else if (uiDimension==3)
+	{
+		ADRESSEXPENSE(9,3,0,0,0,11)
+		if ( (muiIrel+uiIpos<uiImax) && (muiJrel+uiJpos<uiJmax) && (muiKrel+uiKpos<uiKmax) )
+			return (muiIrel+uiIpos) + (muiJrel+uiJpos)*uiImax + (muiKrel+uiKpos)*uiJmax*uiImax;
+		else error->Error(2);
+		return 0;
+	}
+	else return 0;
+}
+
+unsigned int AdrOp::GetPosFromNode(int ny, unsigned int uiNode)
+{
+	while (ny<0) ny+=uiDimension;
+	ny=ny%uiDimension;
+	unsigned int help=uiNode,i=0,j=0,k=0,l=0;
+	i=help%uiImax;
+	help=help/uiImax;
+	j=help%uiJmax;
+	help=help/uiJmax;
+	if (uiKmax>0)
+	{
+		k=help%uiKmax;
+		help=help/uiKmax;
+		l=help;
+	}
+	if (!CheckPos(i,j,k,l)) error->Error(8);
+	ADRESSEXPENSE(0,7,0,0,13,4)
+	switch (ny)
+	{
+	case 0:
+		{
+			return i;
+			break;
+		}
+	case 1:
+		{
+			return j;
+			break;
+		}
+	case 2:
+		{
+			return k;
+			break;
+		}
+	case 3:
+		{
+			return l;
+			break;
+		}
+	}
+	return 0;
+}
+
+double AdrOp::GetNodeVolume(unsigned int uiNode)
+{
+	for (unsigned int n=0; n<uiDimension; n++) if (dGrid[n]==NULL) error->Error(9);
+	double dVol=1;
+	unsigned int uiMax[4]={uiImax,uiJmax,uiKmax,uiLmax};
+	unsigned int uiPos[4]={GetPosFromNode(0,uiNode),GetPosFromNode(1,uiNode),GetPosFromNode(2,uiNode),GetPosFromNode(3,uiNode)};
+	for (unsigned int n=0; n<uiDimension; n++)
+	{
+		if ((uiPos[n]>0) && (uiPos[n]<uiMax[n]-1))
+		{
+			dVol*=0.5*dDeltaUnit*(dGrid[n][uiPos[n]+1]-dGrid[n][uiPos[n]-1]);
+			ADRESSEXPENSE(4,0,1,3,0,4)
+		}
+		else if ((uiPos[n]==0) && (uiPos[n]<uiMax[n]-1))
+		{
+			dVol*=dDeltaUnit*(dGrid[n][uiPos[n]+1]-dGrid[n][uiPos[n]]);
+			ADRESSEXPENSE(3,0,1,2,0,4)
+		}
+		else if ((uiPos[n]>0) && (uiPos[n]==uiMax[n]-1))
+		{
+			dVol*=dDeltaUnit*(dGrid[n][uiPos[n]]-dGrid[n][uiPos[n]-1]);
+			ADRESSEXPENSE(3,0,1,2,0,4)
+		}
+	}
+	return dVol;
+}
+
+double AdrOp::GetIndexWidth(int ny, int index)
+{
+	for (unsigned int n=0; n<uiDimension; n++) if (dGrid[n]==NULL) error->Error(9);
+	double width=0;
+	while (ny<0) ny+=uiDimension;
+	ny=ny%uiDimension;
+	unsigned int uiMax[4]={uiImax,uiJmax,uiKmax,uiLmax};
+	if ((index>0) && (index<(int)uiMax[ny]-1)) width= (dGrid[ny][index+1]-dGrid[ny][index-1])/2*dDeltaUnit;
+	else if ((index==0) && (index<(int)uiMax[ny]-1)) width=(dGrid[ny][index+1]-dGrid[ny][index])*dDeltaUnit;
+	else if ((index>0) && (index==(int)uiMax[ny]-1)) width=(dGrid[ny][index]-dGrid[ny][index-1])*dDeltaUnit;
+	else width= 0;
+	return width;
+}
+
+double AdrOp::GetIndexCoord(int ny, int index)
+{
+	for (unsigned int n=0; n<uiDimension; n++) if (dGrid[n]==NULL) error->Error(9);
+	while (ny<0) ny+=uiDimension;
+	ny=ny%uiDimension;
+	unsigned int uiMax[4]={uiImax,uiJmax,uiKmax,uiLmax};
+	if (index<0) index=0;
+	if (index>=(int)uiMax[ny]) index=uiMax[ny]-1;
+	return dGrid[ny][index]*dDeltaUnit;
+}
+
+double AdrOp::GetIndexDelta(int ny, int index)
+{
+	if (index<0) return GetIndexCoord(ny, 0) - GetIndexCoord(ny, 1);
+	unsigned int uiMax[4]={uiImax,uiJmax,uiKmax,uiLmax};
+	if (index>=(int)uiMax[ny]-1) return GetIndexCoord(ny, (int)uiMax[ny]-2) - GetIndexCoord(ny, (int)uiMax[ny]-1);
+	return GetIndexCoord(ny, index+1) - GetIndexCoord(ny, index);
+}
+
+
+unsigned int AdrOp::Shift(int ny, int step)
+{
+	if (bPosSet==false) error->Error(6);
+	while (ny<0) ny+=uiDimension;
+	ny=ny%uiDimension;
+	switch (ny)
+	{
+	case 0:
+		{
+			iIshift=step;
+//			if ((int)uiIpos+step<0) iIshift=-2*uiIpos-iIshift;
+//			else if ((int)uiIpos+step>=(int)uiImax) iIshift=-1*iIshift+2*(uiImax-1-uiIpos);
+			break;
+		}
+	case 1:
+		{
+			iJshift=step;
+//			if ((int)uiJpos+iJshift<0) iJshift=-2*uiJpos-iJshift;
+//			else if ((int)uiJpos+iJshift>=(int)uiJmax) iJshift=-1*iJshift+2*(uiJmax-1-uiJpos);
+			break;
+		}
+	case 2:
+		{
+			iKshift=step;
+//			if ((int)uiKpos+iKshift<0) iKshift=-2*uiKpos-iKshift;
+//			else if ((int)uiKpos+iKshift>=(int)uiKmax) iKshift=-1*iKshift+2*(uiKmax-1-uiKpos);
+			break;
+		}
+	}
+	ADRESSEXPENSE(1,1,0,0,2,3)
+	return GetPos(iIshift,iJshift,iKshift);
+}
+
+bool AdrOp::CheckShift(int ny, int step)
+{
+	while (ny<0) ny+=uiDimension;
+	ny=ny%uiDimension;
+	int shift[3]={0,0,0};
+	shift[ny]=step;
+	if (this->CheckPos(uiIpos+shift[0],uiJpos+shift[1],uiKpos+shift[2]))
+	{
+		this->Shift(ny,step);
+		return true;
+	}
+	else return false;
+}
+
+unsigned int AdrOp::GetShiftedPos(int ny, int step)
+{
+	if ((ny<0) || (ny>2))
+		return GetPos(iIshift,iJshift,iKshift);
+	int pos[3] = {iIshift, iJshift, iKshift};
+	pos[ny]+=step;
+	return GetPos(pos[0],pos[1],pos[2]);
+}
+
+void AdrOp::ResetShift()
+{
+	iIshift=iJshift=iKshift=0;
+}
+
+unsigned int AdrOp::Iterate(int jump)
+{
+	if (abs(jump)>=(int)uiImax) error->Error(4);
+	i=uiIpos+jump;
+	if (i>=uiImax)
+	{
+		i=i-uiImax;
+		j=uiJpos+1;
+		if (j>=uiJmax)
+		{
+			j=0;
+			if (uiDimension==3)
+			{
+				k=uiKpos+1;
+				if (k>=uiKmax) k=0;
+				uiKpos=k;
+			}
+		}
+		uiIpos=i;
+		uiJpos=j;
+		return GetPos();
+	}
+	else
+	{
+		uiIpos=i;
+		return GetPos();
+	}
+}
+
+unsigned int AdrOp::GetSize()
+{
+	return uiSize;
+}
+
+
+void AdrOp::SetReflection2Node()
+{
+	reflect=true;
+	uiTypeOffset=0;
+}
+
+void AdrOp::SetReflection2Cell()
+{
+	reflect=true;
+	uiTypeOffset=1;
+}
+
+void AdrOp::SetReflectionOff()
+{
+	reflect=false;
+}
+
+AdrOp* AdrOp::AddCellAdrOp()
+{
+	if (clCellAdr!=NULL) return clCellAdr;
+	if (uiDimension==3) clCellAdr = new AdrOp(uiImax-1,uiJmax-1,uiKmax-1);
+	else if (uiDimension==2) clCellAdr = new AdrOp(uiImax-1,uiJmax-1);
+	else clCellAdr=NULL;
+	if (clCellAdr!=NULL)
+	{
+		clCellAdr->SetPos(0,0,0);
+		clCellAdr->SetReflection2Cell();
+	}
+	iCellShift[0]=iCellShift[1]=iCellShift[2]=0;
+	return clCellAdr;
+}
+
+AdrOp* AdrOp::DeleteCellAdrOp()
+{
+	delete clCellAdr;
+	clCellAdr=NULL;
+	return NULL;
+}
+
+unsigned int AdrOp::ShiftCell(int ny, int step)
+{
+	if (clCellAdr==NULL) error->Error(7);
+	while (ny<0) ny+=uiDimension;
+	ny=ny%uiDimension;
+	iCellShift[ny]=step;
+	ADRESSEXPENSE(3,1,0,0,1,2)
+	return clCellAdr->GetPos(uiIpos+iCellShift[0],uiJpos+iCellShift[1],uiKpos+iCellShift[2]);
+}
+
+bool AdrOp::ShiftCellCheck(int ny, int step)
+{
+	return clCellAdr->CheckShift(ny,step);
+}
+
+void AdrOp::ResetCellShift()
+{
+	if (clCellAdr==NULL) error->Error(7);
+	iCellShift[0]=iCellShift[1]=iCellShift[2]=0;
+}
+
+unsigned int AdrOp::GetCellPos(bool incShift)
+{
+	if (bPosSet==false) error->Error(6);
+	if (clCellAdr==NULL) error->Error(7);
+#if EXPENSE_LOG==1
+	if (incShift) ADRESSEXPENSE(3,0,0,0,0,2)
+#endif
+		if (incShift) return clCellAdr->GetPos(uiIpos+iCellShift[0],uiJpos+iCellShift[1],uiKpos+iCellShift[2]);
+		else return clCellAdr->GetPos(uiIpos,uiJpos,uiKpos);
+}
+
+unsigned int AdrOp::GetCellPos(int i, int j, int k)
+{
+	if (bPosSet==false) error->Error(6);
+	return clCellAdr->GetPos(uiIpos+i,uiJpos+j,uiKpos+k);
+}
+
+
+double AdrOp::GetShiftCellVolume(int ny, int step)
+{
+	for (unsigned int n=0; n<uiDimension; n++) if (dGrid[n]==NULL) error->Error(9);
+	int uiMax[4]={(int)uiImax-1,(int)uiJmax-1,(int)uiKmax-1,(int)uiLmax-1};
+	while (ny<0) ny+=uiDimension;
+	ny=ny%uiDimension;
+	iCellShift[ny]=step;
+	int uiPos[4]={(int)uiIpos+iCellShift[0],(int)uiJpos+iCellShift[1],(int)uiKpos+iCellShift[2]};
+	double dVol=1;
+	for (unsigned int n=0; n<uiDimension; ++n)
+	{
+		if (uiMax[n]>0)
+		{
+			while ((uiPos[n]<0) || (uiPos[n]>=uiMax[n]))
+			{
+				if (uiPos[n]<0) uiPos[n]=-1*uiPos[n]-1;
+				else if (uiPos[n]>=uiMax[n]) uiPos[n]=uiMax[n]-(uiPos[n]-uiMax[n]+1);
+			}
+			dVol*=(dGrid[n][uiPos[n]+1]-dGrid[n][uiPos[n]])*dDeltaUnit;
+		}
+	}
+	return dVol;
+}
+
+
+deltaAdrOp::deltaAdrOp(unsigned int max)
+{
+	uiMax=max;
+}
+
+deltaAdrOp::~deltaAdrOp()
+{
+}
+
+void deltaAdrOp::SetMax(unsigned int max)
+{
+	uiMax=max;
+}
+
+unsigned int deltaAdrOp::GetAdr(int pos)
+{
+	if (uiMax==1) return 0;
+	if (pos<0)	pos=pos*-1;
+	else if (pos>(int)uiMax-1) pos=2*(uiMax-1)-pos+1;
+	if ((pos<0) || (pos>(int)uiMax-1))
+	{
+		fprintf(stderr," Error exiting... ");
+		getchar();
+		exit(-1);
+	}
+	return pos;
+}
+
+
diff --git a/openEMS/tools/AdrOp.h b/openEMS/tools/AdrOp.h
new file mode 100644
index 0000000..e4a690b
--- /dev/null
+++ b/openEMS/tools/AdrOp.h
@@ -0,0 +1,149 @@
+/*
+*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+/*!
+\class AdrOp
+\author Thorsten Liebig
+\version $Revision: 1.10 $
+\date $Date: 2006/10/29 18:50:44 $
+*/
+
+#ifndef ADROP_H
+#define ADROP_H
+
+#include <iostream>
+#include <stdio.h>
+#include <stdlib.h>
+#include <vector>
+#include "ExpenseLog.h"
+#include "ErrorMsg.h"
+
+class AdrOp
+{
+public:
+	///Constructor, define dimension/size here
+	AdrOp(unsigned int muiImax, unsigned int muiYmax, unsigned int muiKmax=0, unsigned int muiLmax=0);
+	///Copy-Constructor
+	AdrOp(AdrOp* origOP);
+	///Deconstructor
+	virtual ~AdrOp();
+	///Set the current n-dim position, get 1-dim array position as return value
+	/*!A position has to be set or all other methodes will case error! \n The methode will exit with error message if invalid position is set! \sa ErrorMsg */
+	unsigned int SetPos(unsigned int muiIpos, unsigned int muiJpos, unsigned int muiKpos=0, unsigned int muiLpos=0);
+
+	bool SetPosChecked(unsigned int muiIpos, unsigned int muiJpos, unsigned int muiKpos=0, unsigned int muiLpos=0);
+
+	void SetGrid(double *gridI,double *gridJ,double *gridK=NULL,double *gridL=NULL);
+	void SetGridDelta(double delta) {this->dDeltaUnit=delta;};
+
+	bool CheckPos(unsigned int muiIpos, unsigned int muiJpos, unsigned int muiKpos=0, unsigned int muiLpos=0);
+	bool CheckRelativePos(int muiIrel=0,int muiJrel=0,int muiKrel=0, int muiLrel=0);
+	///will return current 1-dim position, in addition to a relative n-dim shift
+	/*!In case of crossing the boundaries, activate reflection or an error will be invoked\sa SetReflection2Node \sa SetReflection2Cell \sa SetReflectionOff */
+	unsigned int GetPos(int muiIrel=0,int muiJrel=0,int muiKrel=0, int muiLrel=0);
+
+	double GetNodeVolume(unsigned int uiNode);
+
+	double GetIndexWidth(int ny, int index);
+	double GetIndexCoord(int ny, int index);
+	///Get the gird delta at the given  index of direction ny. (if index<0 return negative value as index=0 would give, if index>=max-1 returns negative value as index=max-2 would give)
+	double GetIndexDelta(int ny, int index);
+
+//	double GetCellVolume(unsigned int uiCell);
+
+	unsigned int GetPosFromNode(int ny, unsigned int uiNode);
+	///Set a shift in ny direction (e.g. 0 for i-direction)
+	/*!Shift set by this methode will be ignored by methode GetPos*/
+	unsigned int Shift(int ny, int step);
+	///Set a checked shift in ny direction (e.g. 0 for i-direction)
+	/*!Shift set by this methode will be ignored by methode GetPos*/
+	bool CheckShift(int ny, int step);
+	///Returns the current 1-dim position including shift by methode "Shift" + additional (transitory) shift
+	unsigned int GetShiftedPos(int ny=-1, int step=0);
+	///Reset shift set by "Shift"-methode
+	void ResetShift();
+	///Iterates through AdrOp; --- obsolete ---
+	unsigned int Iterate(int jump=1);
+	///Retruns size of array
+	unsigned int GetSize();
+	///Set mode to reflect by node
+	/*!1D-example (6 nodes): \image html node_reflect.PNG order: 0,1,2,3,4,5,4,3,...*/
+	void SetReflection2Node();
+	///Set mode to reflect by cell
+	/*!1D-example (5 cells): \image html cell_reflect.PNG order: 0,1,2,3,4,4,3,...*/
+	void SetReflection2Cell();
+	///Deactivate reflection (default)
+	void SetReflectionOff();
+	///Add a cell adress operator (dimensions: i-1 j-1 k-1 l-1)
+	/*!\image html cells_nodes.png */
+	AdrOp* AddCellAdrOp();
+
+	AdrOp* GetCellAdrOp() {return clCellAdr;};
+
+	///Deconstructed cell adress operator if no longer needed
+	AdrOp* DeleteCellAdrOp();
+	///Shift cell in ny dircetion; cell reflection is active
+	unsigned int ShiftCell(int ny, int step);
+	///Shift cell in ny dircetion; cell reflection is active; return check
+	bool ShiftCellCheck(int ny, int step);
+	///Reset cell shift
+	void ResetCellShift();
+	///Get current cell position from cell adress operator
+	unsigned int GetCellPos(bool incShift=true);
+	///Get cell position from cell adress operator
+	unsigned int GetCellPos(int i, int j, int k=0);
+	//get volume of cell; incl shift
+	double GetShiftCellVolume(int ny, int step);
+
+
+	void SetDebugOn() {this->bDebug=true;};
+	void SetDebugOff() {this->bDebug=false;};
+
+protected:
+	AdrOp *clCellAdr;
+	unsigned int uiDimension;
+	unsigned int uiSize;
+	unsigned int uiImax,uiJmax,uiKmax,uiLmax;
+	unsigned int uiIpos, uiJpos, uiKpos, uiLpos;
+	double *dGrid[4];
+	double dDeltaUnit;
+	int iIshift, iJshift, iKshift;
+	int iCellShift[3];
+	unsigned int i,j,k,l;
+	bool reflect;
+	unsigned int uiTypeOffset;
+
+	bool bPosSet;
+	bool bDebug;
+	ErrorMsg *error;
+};
+
+
+class deltaAdrOp
+{
+public:
+	deltaAdrOp(unsigned int max=0);
+	virtual ~deltaAdrOp();
+	void SetMax(unsigned int max);
+	unsigned int GetAdr(int pos);
+
+protected:
+	unsigned int uiMax;
+};
+
+
+#endif // ADROP_H
diff --git a/openEMS/tools/CMakeLists.txt b/openEMS/tools/CMakeLists.txt
new file mode 100644
index 0000000..8b62b8a
--- /dev/null
+++ b/openEMS/tools/CMakeLists.txt
@@ -0,0 +1,28 @@
+
+set(SOURCES
+  ${SOURCES}
+  ${CMAKE_CURRENT_SOURCE_DIR}/AdrOp.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/ErrorMsg.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/array_ops.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/global.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/hdf5_file_reader.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/hdf5_file_writer.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/sar_calculation.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/useful.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/vtk_file_writer.cpp
+  PARENT_SCOPE
+)
+
+#set(HEADERS
+#  constants.h
+#  array_ops.h
+#  global.h
+#  useful.h
+#  aligned_allocator.h
+#  hdf5_file_reader.h
+#  hdf5_file_writer.h
+#)
+
+# tools lib
+#add_library(tools STATIC ${SOURCES} )
+
diff --git a/openEMS/tools/ErrorMsg.cpp b/openEMS/tools/ErrorMsg.cpp
new file mode 100644
index 0000000..9a1974c
--- /dev/null
+++ b/openEMS/tools/ErrorMsg.cpp
@@ -0,0 +1,98 @@
+/*
+*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "ErrorMsg.h"
+
+ErrorMsg::ErrorMsg(unsigned int NoMessage)
+{
+	NoMsg=NoMessage;
+	if (NoMsg>0) Msg = new char*[NoMsg];
+	if (Msg==NULL)
+	{
+		fprintf(stderr,"Memory allocation failed!! exiting...");
+		exit(1);
+	}
+	for (unsigned int i=0; i<NoMsg; i++) Msg[i]=NULL;
+}
+
+ErrorMsg::~ErrorMsg()
+{
+	for (unsigned int i=0; i<NoMsg; i++)
+	{
+		delete[] Msg[i];
+		Msg[i]=NULL;
+	};
+	delete[] Msg;
+	Msg=NULL;
+}
+
+void ErrorMsg::SetMsg(unsigned int nr, const char *Message)
+{
+	if ((nr<1) || (nr>NoMsg) || (Message==NULL)) ownError();
+	Msg[nr-1] = new char[strlen(Message)+1];
+	if (Msg[nr-1]==NULL)
+	{
+		fprintf(stderr,"Memory allocation failed!! exiting...");
+		exit(1);
+	}
+	Msg[nr-1]=strcpy(Msg[nr-1],Message);
+}
+
+void ErrorMsg::Error(unsigned int nr,char *chAddMsg)
+{
+	if ((nr>0) && (nr<=NoMsg))
+	{
+		if (Msg[nr-1]!=NULL) fprintf(stderr,"%s",Msg[nr-1]);
+		else fprintf(stderr,"unkown error occured!! Error code: %d exiting...",nr);
+		if (chAddMsg!=NULL) fprintf(stderr,"%s",chAddMsg);
+		getchar();
+		exit(nr);
+	}
+	else
+	{
+		fprintf(stderr,"unkown error occured!! Error code: %d exiting...",nr);
+		getchar();
+		exit(nr);
+	}
+}
+
+void ErrorMsg::Error(unsigned int nr,int addNr)
+{
+	if ((nr>0) && (nr<=NoMsg))
+	{
+		if (Msg[nr-1]!=NULL) fprintf(stderr,"%s",Msg[nr-1]);
+		else fprintf(stderr,"unkown error occured!! Error code: %d exiting...",nr);
+		fprintf(stderr,"%d",addNr);
+		getchar();
+		exit(nr);
+	}
+	else
+	{
+		fprintf(stderr,"unkown error occured!! Error code: %d exiting...",nr);
+		getchar();
+		exit(nr);
+	}
+}
+
+void ErrorMsg::ownError(void)
+{
+	fprintf(stdout," Error occured by using Error Message class!! ... exiting...");
+	exit(-1);
+}
diff --git a/openEMS/tools/ErrorMsg.h b/openEMS/tools/ErrorMsg.h
new file mode 100644
index 0000000..fa63c93
--- /dev/null
+++ b/openEMS/tools/ErrorMsg.h
@@ -0,0 +1,50 @@
+/*
+*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+/*!
+\class ErrorMsg
+\author Thorsten Liebig
+\version $Revision: 1.2 $
+\date $Date: 2006/01/25 11:47:07 $
+*/
+
+#ifndef _ERRORMSG_H_
+#define _ERRORMSG_H_
+
+class ErrorMsg
+{
+public:
+	///Constructor defines number of error messages
+	ErrorMsg(unsigned int NoMessage=0);
+	///Deconstructor
+	virtual ~ErrorMsg();
+	///Methode for defining error messages
+	/*! \param nr Number of defining error message \param *Message Set error message string  \sa Error */
+	void SetMsg(unsigned int nr, const char *Message);
+	///Call an error message. Will exit the program!
+	/*! \param nr Number of called error message. default is 0 \sa SetMsg*/
+	void Error(unsigned int nr=0,char *chAddMsg=0);
+
+	void Error(unsigned int nr,int addNr);
+
+protected:
+	void ownError(void);
+	unsigned int NoMsg;
+	char **Msg;
+};
+
+#endif //_ERRORMSG_H_
diff --git a/openEMS/tools/ExpenseLog.cpp b/openEMS/tools/ExpenseLog.cpp
new file mode 100644
index 0000000..95a0d4c
--- /dev/null
+++ b/openEMS/tools/ExpenseLog.cpp
@@ -0,0 +1,145 @@
+#include "ExpenseLog.h"
+
+ExpenseModule::ExpenseModule(const char* moduleName)
+{
+	chModuleName=moduleName;
+	uiDoubleAdditions=uiDoubleMultiplications=uiIntAdditions=uiIntMultiplications=uiAssignments=uiBoolOp=0;
+	uiMrdDA=uiMrdDM=uiMrdIA=uiMrdIM=uiMrdAssign=uiMrdBO=0;
+}
+
+ExpenseModule::~ExpenseModule() {};
+
+void ExpenseModule::Clear()
+{
+	uiDoubleAdditions=uiDoubleMultiplications=uiIntAdditions=uiIntMultiplications=uiAssignments=uiBoolOp=0;
+	uiMrdDA=uiMrdDM=uiMrdIA=uiMrdIM=uiMrdAssign=uiMrdBO=0;
+}
+
+void ExpenseModule::AddDoubleAdditons(unsigned int number)
+{
+	uiDoubleAdditions+=number;
+	if (uiDoubleAdditions>=MRD)
+	{
+		uiDoubleAdditions-=MRD;
+		++uiMrdDA;
+	}
+}
+
+void ExpenseModule::AddDoubleMultiplications(unsigned int number)
+{
+	uiDoubleMultiplications+=number;
+	if (uiDoubleMultiplications>=MRD)
+	{
+		uiDoubleMultiplications-=MRD;
+		++uiMrdDM;
+	}
+}
+
+void ExpenseModule::AddIntAdditons(unsigned int number)
+{
+	uiIntAdditions+=number;
+	if (uiIntAdditions>=MRD)
+	{
+		uiIntAdditions-=MRD;
+		++uiMrdIA;
+	}
+}
+
+void ExpenseModule::AddIntMultiplications(unsigned int number)
+{
+	uiIntMultiplications+=number;
+	if (uiIntMultiplications>=MRD)
+	{
+		uiIntMultiplications-=MRD;
+		++uiMrdIM;
+	}
+}
+
+void ExpenseModule::AddAssignments(unsigned int number)
+{
+	uiAssignments+=number;
+	if (uiAssignments>=MRD)
+	{
+		uiAssignments-=MRD;
+		++uiMrdAssign;
+	}
+}
+
+void ExpenseModule::AddBoolOperations(unsigned int number)
+{
+	uiBoolOp+=number;
+	if (uiBoolOp>=MRD)
+	{
+		uiBoolOp-=MRD;
+		++uiMrdBO;
+	}
+}
+
+void ExpenseModule::AddOperations(unsigned int IntAdd, unsigned int IntMul, unsigned int DoubleAdd, unsigned int DoubleMul, unsigned int Assigns, unsigned int BoolOp)
+{
+	this->AddIntAdditons(IntAdd);
+	this->AddIntMultiplications(IntMul);
+	this->AddDoubleAdditons(DoubleAdd);
+	this->AddDoubleMultiplications(DoubleMul);
+	this->AddAssignments(Assigns);
+	this->AddBoolOperations(BoolOp);
+}
+
+void ExpenseModule::PrintfSelf(FILE* file)
+{
+	fprintf(file,"\n***********\n Module:  %s\n Additions:\n Double: %3.0d%9d\tInteger: %3.0d%9d",chModuleName,uiMrdDA,uiDoubleAdditions,uiMrdIA,uiIntAdditions);
+	fprintf(file,"\n\n Multiplications:\n Double: %3.0d%9d\tInteger: %3.0d%9d\n",uiMrdDM,uiDoubleMultiplications,uiMrdIM,uiIntMultiplications);
+	fprintf(file,"\n Assignments: %3.0d%9d\tBool Operations: %3.0d%9d\n",uiMrdAssign,uiAssignments,uiMrdBO,uiBoolOp);
+	fprintf(file,"\n***********\n");
+}
+
+
+
+
+/***********************************************************************************************************************/
+
+ExpenseLog::ExpenseLog(void)
+{
+}
+
+ExpenseLog::~ExpenseLog(void)
+{
+	for (size_t i=0; i<vModules.size(); ++i)
+	{
+		delete vModules.at(i);
+		vModules.at(i)=NULL;
+	}
+	vModules.clear();
+}
+
+ExpenseModule* ExpenseLog::AddModule(const char* name)
+{
+	ExpenseModule* newModule = new ExpenseModule(name);
+	vModules.push_back(newModule);
+	return newModule;
+}
+
+void ExpenseLog::PrintAll(FILE *file)
+{
+	double totalAdd=0,totalMul=0,totalBool=0,totalAssign=0;
+	fprintf(stderr,"\n ----------------\n Expense Log PrintOut\n Nr of Modules: %d\n",vModules.size());
+	for (size_t i=0; i<vModules.size(); ++i)
+	{
+		totalAdd+=((double)vModules.at(i)->uiIntAdditions+(double)vModules.at(i)->uiDoubleAdditions) + 1e9*((double)vModules.at(i)->uiMrdIA+(double)vModules.at(i)->uiMrdIA);
+		totalMul+=((double)vModules.at(i)->uiIntMultiplications+(double)vModules.at(i)->uiDoubleMultiplications) + 1e9*(double)(vModules.at(i)->uiMrdIM+vModules.at(i)->uiMrdIM);
+		totalBool+=(double)vModules.at(i)->uiBoolOp + 1e9*(double)vModules.at(i)->uiMrdBO;
+		totalAssign+=(double)vModules.at(i)->uiAssignments + 1e9*(double)vModules.at(i)->uiMrdAssign;
+		vModules.at(i)->PrintfSelf(file);
+	}
+	fprintf(stderr," Total:\n Additions: %e Multiplications: %e\n Bool Operations: %e Assignments: %e\n",totalAdd,totalMul,totalBool,totalAssign);
+	fprintf(stderr,"\n ----------------\n");
+}
+
+void ExpenseLog::ClearAll()
+{
+	for (size_t i=0; i<vModules.size(); ++i)
+	{
+		vModules.at(i)->Clear();
+	}
+}
+
diff --git a/openEMS/tools/ExpenseLog.h b/openEMS/tools/ExpenseLog.h
new file mode 100644
index 0000000..f8cc1ba
--- /dev/null
+++ b/openEMS/tools/ExpenseLog.h
@@ -0,0 +1,92 @@
+#pragma once
+#include <stdio.h>
+#include <stdlib.h>
+#include <vector>
+
+using namespace std;
+
+#define EXPENSE_LOG 0
+#define MRD 1000000000
+
+#if EXPENSE_LOG==1
+
+#define EXPENSE_DEFINE \
+ExpenseLog EL; \
+ExpenseModule* EngineExpense=EL.AddModule("Static Engine Expenses"); \
+ExpenseModule* PPExpense=EL.AddModule("Static Post Processing"); \
+ExpenseModule* AdrOpExpense=EL.AddModule("Adress Operator");
+#define EXTERN_EXPENSE_DEFINE extern ExpenseLog EL;
+#define ENGINEEXPENSE_DEFINE extern ExpenseModule* EngineExpense;
+#define POSTPROCEXPENSE_DEFINE extern ExpenseModule* PPExpense;
+#define ADREXPENSE_DEFINE extern ExpenseModule* AdrOpExpense;
+#define ENGINEEXPENSE(IA,IM,DA,DM,AS,BO) EngineExpense->AddOperations((IA),(IM),(DA),(DM),(AS),(BO));
+#define POSTPROCEXPENSE(IA,IM,DA,DM,AS,BO) PPExpense->AddOperations((IA),(IM),(DA),(DM),(AS),(BO));
+#define ADRESSEXPENSE(IA,IM,DA,DM,AS,BO) AdrOpExpense->AddOperations((IA),(IM),(DA),(DM),(AS),(BO));
+#define EXPENSEPRINT EL.PrintAll(stderr);
+#define EXPENSECLEAR EL.ClearAll();
+#else
+
+#define EXPENSE_DEFINE
+#define EXTERN_EXPENSE_DEFINE
+#define ENGINEEXPENSE_DEFINE
+#define POSTPROCEXPENSE_DEFINE
+#define ADREXPENSE_DEFINE
+#define ENGINEEXPENSE(IA,IM,DA,DM,AS,BO)
+#define POSTPROCEXPENSE(IA,IM,DA,DM,AS,BO)
+#define ADRESSEXPENSE(IA,IM,DA,DM,AS,BO)
+#define EXPENSEPRINT
+#define EXPENSECLEAR
+#endif
+
+class ExpenseModule
+{
+	friend class ExpenseLog;
+public:
+	ExpenseModule(const char* moduleName);
+	~ExpenseModule();
+
+	void Clear();
+
+	void AddDoubleAdditons(unsigned int number);
+	void AddDoubleMultiplications(unsigned int number);
+
+	void AddIntAdditons(unsigned int number);
+	void AddIntMultiplications(unsigned int number);
+
+	void AddAssignments(unsigned int number);
+	void AddBoolOperations(unsigned int number);
+
+	void AddOperations(unsigned int IntAdd, unsigned int IntMul, unsigned int DoubleAdd, unsigned int DoubleMul, unsigned int Assigns, unsigned int BoolOp);
+
+	void PrintfSelf(FILE* file=stdout);
+
+protected:
+	const char* chModuleName;
+	unsigned int uiDoubleAdditions;
+	unsigned int uiDoubleMultiplications;
+	unsigned int uiIntAdditions;
+	unsigned int uiIntMultiplications;
+	unsigned int uiAssignments;
+	unsigned int uiBoolOp;
+	unsigned int uiMrdDA;
+	unsigned int uiMrdDM;
+	unsigned int uiMrdIA;
+	unsigned int uiMrdIM;
+	unsigned int uiMrdAssign;
+	unsigned int uiMrdBO;
+};
+
+class ExpenseLog
+{
+public:
+	ExpenseLog(void);
+	~ExpenseLog(void);
+
+	ExpenseModule* AddModule(const char* name);
+	void PrintAll(FILE *file=stdout);
+	void ClearAll();
+protected:
+	vector<ExpenseModule*> vModules;
+};
+
+
diff --git a/openEMS/tools/aligned_allocator.h b/openEMS/tools/aligned_allocator.h
new file mode 100644
index 0000000..efabe7d
--- /dev/null
+++ b/openEMS/tools/aligned_allocator.h
@@ -0,0 +1,173 @@
+/*
+*	Copyright (C) 2010 Sebastian Held (sebastian.held@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+// based on http://blogs.msdn.com/b/vcblog/archive/2008/08/28/the-aligned_allocator.aspx
+// from Stephan T. Lavavej
+
+
+// The following headers are required for all allocators.
+#include <stddef.h>  // Required for size_t and ptrdiff_t and NULL
+#include <new>       // Required for placement new and std::bad_alloc
+#include <stdexcept> // Required for std::length_error
+
+
+#ifdef WIN32
+#define __MSVCRT_VERSION__ 0x0700
+#include <malloc.h>
+#define MEMALIGN( array, alignment, size ) !(*array = _aligned_malloc( size, alignment ))
+#define FREE( array ) _aligned_free( array )
+#else
+#define MEMALIGN( array, alignment, size ) posix_memalign( array, alignment, size )
+#define FREE( array ) free( array )
+#endif
+
+
+template <typename T> class aligned_allocator
+{
+public:
+	// The following will be the same for virtually all allocators.
+	typedef T * pointer;
+	typedef const T * const_pointer;
+	typedef T& reference;
+	typedef const T& const_reference;
+	typedef T value_type;
+	typedef size_t size_type;
+	typedef ptrdiff_t difference_type;
+
+	T * address(T& r) const
+	{
+		return &r;
+	}
+
+	const T * address(const T& s) const
+	{
+		return &s;
+	}
+
+	size_t max_size() const
+	{
+		// The following has been carefully written to be independent of
+		// the definition of size_t and to avoid signed/unsigned warnings.
+		return (static_cast<size_t>(0) - static_cast<size_t>(1)) / sizeof(T);
+	}
+
+	// The following must be the same for all allocators.
+	template <typename U> struct rebind
+	{
+		typedef aligned_allocator<U> other;
+	};
+
+	bool operator!=(const aligned_allocator& other) const
+	{
+		return !(*this == other);
+	}
+
+	void construct(T * const p, const T& t) const
+	{
+		void * const pv = static_cast<void *>(p);
+		new (pv) T(t);
+	}
+
+	void destroy(T * const p) const; // Defined below.
+
+	// Returns true if and only if storage allocated from *this
+	// can be deallocated from other, and vice versa.
+	// Always returns true for stateless allocators.
+	bool operator==(const aligned_allocator& other) const
+	{
+		return true;
+	}
+
+	// Default constructor, copy constructor, rebinding constructor, and destructor.
+	// Empty for stateless allocators.
+	aligned_allocator() { }
+	aligned_allocator(const aligned_allocator&) { }
+	template <typename U> aligned_allocator(const aligned_allocator<U>&) { }
+	~aligned_allocator() { }
+
+	// The following will be different for each allocator.
+	T * allocate(const size_t n) const
+	{
+//		std::cout << "Allocating " << n << (n == 1 ? " object" : "objects")	<< " of size " << sizeof(T) << "." << std::endl;
+		// The return value of allocate(0) is unspecified.
+		// aligned_allocator returns NULL in order to avoid depending
+		// on malloc(0)'s implementation-defined behavior
+		// (the implementation can define malloc(0) to return NULL,
+		// in which case the bad_alloc check below would fire).
+		// All allocators can return NULL in this case.
+		if (n == 0)
+		{
+			return NULL;
+		}
+
+		// All allocators should contain an integer overflow check.
+		// The Standardization Committee recommends that std::length_error
+		// be thrown in the case of integer overflow.
+		if (n > max_size())
+		{
+			throw std::length_error("aligned_allocator<T>::allocate() - Integer overflow.");
+		}
+
+		// Allocators should throw std::bad_alloc in the case of memory allocation failure.
+		void * pv;
+		if (MEMALIGN( &pv, 16, n * sizeof(T)))
+			throw std::bad_alloc();
+
+		return static_cast<T *>(pv);
+	}
+
+	void deallocate(T * const p, const size_t n) const
+	{
+//		std::cout << "Deallocating " << n << (n == 1 ? " object" : "objects") << " of size " << sizeof(T) << "." << std::endl;
+		// aligned_allocator wraps free().
+		UNUSED(n);
+		FREE(p);
+	}
+
+	// The following will be the same for all allocators that ignore hints.
+	template <typename U> T * allocate(const size_t n, const U * /* const hint */) const
+	{
+		return allocate(n);
+	}
+
+	// Allocators are not required to be assignable, so
+	// all allocators should have a private unimplemented
+	// assignment operator. Note that this will trigger the
+	// off-by-default (enabled under /Wall) warning C4626
+	// "assignment operator could not be generated because a
+	// base class assignment operator is inaccessible" within
+	// the STL headers, but that warning is useless.
+
+private:
+	aligned_allocator& operator=(const aligned_allocator&);
+};
+
+// A compiler bug causes it to believe that p->~T() doesn't reference p.
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable: 4100) // unreferenced formal parameter
+#endif
+
+// The definition of destroy() must be the same for all allocators.
+template <typename T> void aligned_allocator<T>::destroy(T * const p) const
+{
+	p->~T();
+}
+
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
diff --git a/openEMS/tools/array_ops.cpp b/openEMS/tools/array_ops.cpp
new file mode 100644
index 0000000..8710b3c
--- /dev/null
+++ b/openEMS/tools/array_ops.cpp
@@ -0,0 +1,144 @@
+/*
+*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "array_ops.h"
+#include <ostream>
+
+using namespace std;
+
+#ifdef WIN32
+#define __MSVCRT_VERSION__ 0x0700
+#include <malloc.h>
+//(void**)&array, 16, sizeof(typeof(f4vector**))*numLines[0]
+#define MEMALIGN( array, alignment, size ) !(*array = _aligned_malloc( size, alignment ))
+#define FREE( array ) _aligned_free( array )
+#else
+#define MEMALIGN( array, alignment, size ) posix_memalign( array, alignment, size )
+#define FREE( array ) free( array )
+#endif
+
+void Delete1DArray_v4sf(f4vector* array)
+{
+	if (array==NULL) return;
+	FREE( array );
+}
+
+
+void Delete3DArray_v4sf(f4vector*** array, const unsigned int* numLines)
+{
+	if (array==NULL) return;
+	unsigned int pos[3];
+	for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
+	{
+		for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
+		{
+			FREE( array[pos[0]][pos[1]] );
+			//delete[] array[pos[0]][pos[1]];
+		}
+		FREE( array[pos[0]] );
+		//delete[] array[pos[0]];
+	}
+	FREE( array );
+	//delete[] array;
+}
+
+void Delete_N_3DArray_v4sf(f4vector**** array, const unsigned int* numLines)
+{
+	if (array==NULL) return;
+	for (int n=0; n<3; ++n)
+	{
+		Delete3DArray_v4sf(array[n],numLines);
+	}
+	FREE( array );
+	//delete[] array;
+}
+
+f4vector* Create1DArray_v4sf(const unsigned int numLines)
+{
+	f4vector* array=NULL;
+	if (MEMALIGN( (void**)&array, 16, sizeof(typeof(f4vector))*numLines ))
+	{
+		cerr << "cannot allocate aligned memory" << endl;
+		exit(3);
+	}
+	for (unsigned int pos=0; pos<numLines; ++pos)
+	{
+		array[pos].f[0] = 0;
+		array[pos].f[1] = 0;
+		array[pos].f[2] = 0;
+		array[pos].f[3] = 0;
+	}
+	return array;
+}
+
+//! \brief this function allocates a 3D array, which is aligned to 16 byte
+f4vector*** Create3DArray_v4sf(const unsigned int* numLines)
+{
+	unsigned int numZ = ceil((double)numLines[2]/4.0);
+
+	f4vector*** array=NULL;
+	unsigned int pos[3];
+	if (MEMALIGN( (void**)&array, 16, sizeof(typeof(f4vector**))*numLines[0] ))
+	{
+		cerr << "cannot allocate aligned memory" << endl;
+		exit(3);
+	}
+	//array = new f4vector**[numLines[0]];
+	for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
+	{
+		if (MEMALIGN( (void**)&array[pos[0]], 16, sizeof(typeof(f4vector*))*numLines[1] ))
+		{
+			cerr << "cannot allocate aligned memory" << endl;
+			exit(3);
+		}
+		//array[pos[0]] = new f4vector*[numLines[1]];
+		for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
+		{
+			if (MEMALIGN( (void**)&array[pos[0]][pos[1]], 16, sizeof(typeof(f4vector))*numZ ))
+			{
+				cerr << "cannot allocate aligned memory" << endl;
+				exit(3);
+			}
+			//array[pos[0]][pos[1]] = new f4vector[numZ];
+			for (pos[2]=0; pos[2]<numZ; ++pos[2])
+			{
+				array[pos[0]][pos[1]][pos[2]].f[0] = 0;
+				array[pos[0]][pos[1]][pos[2]].f[1] = 0;
+				array[pos[0]][pos[1]][pos[2]].f[2] = 0;
+				array[pos[0]][pos[1]][pos[2]].f[3] = 0;
+			}
+		}
+	}
+	return array;
+}
+
+f4vector**** Create_N_3DArray_v4sf(const unsigned int* numLines)
+{
+	f4vector**** array=NULL;
+	if (MEMALIGN( (void**)&array, 16, sizeof(typeof(f4vector***))*3 ))
+	{
+		cerr << "cannot allocate aligned memory" << endl;
+		exit(3);
+	}
+	//array = new f4vector***[3];
+	for (int n=0; n<3; ++n)
+	{
+		array[n]=Create3DArray_v4sf(numLines);
+	}
+	return array;
+}
+
diff --git a/openEMS/tools/array_ops.h b/openEMS/tools/array_ops.h
new file mode 100644
index 0000000..5c2fd7f
--- /dev/null
+++ b/openEMS/tools/array_ops.h
@@ -0,0 +1,196 @@
+/*
+*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef ARRAY_OPS_H
+#define ARRAY_OPS_H
+
+#ifdef __SIZEOF_FLOAT__
+#if __SIZEOF_FLOAT__ != 4
+#error wrong size of float
+#endif
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <string>
+#include <math.h>
+#include "constants.h"
+
+typedef float v4sf __attribute__ ((vector_size (16))); // vector of four single floats
+typedef   int v4si __attribute__ ((vector_size (4*sizeof(int)))); // vector of four single ints
+
+union f4vector
+{
+	v4sf v;
+	float f[4];
+};
+
+void Delete1DArray_v4sf(f4vector* array);
+void Delete3DArray_v4sf(f4vector*** array, const unsigned int* numLines);
+void Delete_N_3DArray_v4sf(f4vector**** array, const unsigned int* numLines);
+f4vector* Create1DArray_v4sf(const unsigned int numLines);
+f4vector*** Create3DArray_v4sf(const unsigned int* numLines);
+f4vector**** Create_N_3DArray_v4sf(const unsigned int* numLines);
+
+// *************************************************************************************
+// templates
+// *************************************************************************************
+template <typename T>
+T** Create2DArray(const unsigned int* numLines)
+{
+	T** array=NULL;
+	unsigned int pos[3];
+	array = new T*[numLines[0]];
+	for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
+	{
+		array[pos[0]] = new T[numLines[1]];
+		for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
+		{
+			array[pos[0]][pos[1]] = 0;
+		}
+	}
+	return array;
+}
+
+template <typename T>
+void Delete2DArray(T** array, const unsigned int* numLines)
+{
+	if (array==NULL) return;
+	unsigned int pos[3];
+	for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
+	{
+		delete[] array[pos[0]];
+	}
+	delete[] array;
+}
+
+template <typename T>
+inline T& Access_N_3DArray(T**** array, unsigned int n, unsigned int* pos)
+{
+	return array[n][pos[0]][pos[1]][pos[2]];
+}
+
+template <typename T>
+inline T& Access_N_3DArray(T**** array, unsigned int n, unsigned int x, unsigned int y, unsigned int z )
+{
+	return array[n][x][y][z];
+}
+
+template <typename T>
+T*** Create3DArray(const unsigned int* numLines)
+{
+	T*** array=NULL;
+	unsigned int pos[3];
+	array = new T**[numLines[0]];
+	for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
+	{
+		array[pos[0]] = new T*[numLines[1]];
+		for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
+		{
+			array[pos[0]][pos[1]] = new T[numLines[2]];
+			for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
+			{
+				array[pos[0]][pos[1]][pos[2]] = 0;
+			}
+		}
+	}
+	return array;
+}
+
+template <typename T>
+T*** Copy3DArray(T*** array_in, T*** array_out, const unsigned int* numLines)
+{
+	if (array_out==NULL)
+		array_out = Create3DArray<T>(numLines);
+	unsigned int pos[3];
+	for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
+		for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
+			for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
+				array_out[pos[0]][pos[1]][pos[2]] = array_in[pos[0]][pos[1]][pos[2]];
+	return array_out;
+}
+
+template <typename T>
+T**** Create_N_3DArray(const unsigned int* numLines)
+{
+	T**** array=NULL;
+	array = new T***[3];
+	for (int n=0; n<3; ++n)
+	{
+		array[n]=Create3DArray<T>( numLines );
+	}
+	return array;
+}
+
+template <typename T>
+T**** Copy_N_3DArray(T**** array_in, T**** array_out, const unsigned int* numLines)
+{
+	if (array_out==NULL)
+		array_out = Create_N_3DArray<T>(numLines);
+	for (int n=0; n<3; ++n)
+		array_out[n]=Copy3DArray<T>( array_in[n], array_out[n], numLines);
+	return array_out;
+}
+
+template <typename T>
+void Delete3DArray(T*** array, const unsigned int* numLines)
+{
+	if (!array) return;
+	unsigned int pos[3];
+	for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
+	{
+		for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
+		{
+			delete[] array[pos[0]][pos[1]];
+		}
+		delete[] array[pos[0]];
+	}
+	delete[] array;
+}
+
+template <typename T>
+void Delete_N_3DArray(T**** array, const unsigned int* numLines)
+{
+	if (!array) return;
+	for (int n=0; n<3; ++n)
+	{
+		Delete3DArray<T>(array[n],numLines);
+	}
+	delete[] array;
+}
+
+template <typename T>
+void Dump_N_3DArray2File(std::ostream &file, T**** array, const unsigned int* numLines)
+{
+	unsigned int pos[3];
+	for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
+	{
+		for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
+		{
+			for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
+			{
+				file << pos[0] << "\t" << pos[1] << "\t" << pos[2];
+				for (int n=0; n<3; ++n)
+					file << "\t" << (float)array[n][pos[0]][pos[1]][pos[2]];
+				file << std::endl;
+			}
+		}
+	}
+}
+
+#endif // ARRAY_OPS_H
diff --git a/openEMS/tools/constants.h b/openEMS/tools/constants.h
new file mode 100644
index 0000000..9707fdf
--- /dev/null
+++ b/openEMS/tools/constants.h
@@ -0,0 +1,29 @@
+/*
+*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef CONSTANTS_H
+#define CONSTANTS_H
+
+#define FDTD_FLOAT float
+
+#define __EPS0__ 8.85418781762e-12
+#define __MUE0__ 1.256637062e-6
+#define __C0__ 299792458
+#define __Z0__ 376.730313461
+#define PI 3.141592653589793238462643383279
+
+#endif // CONSTANTS_H
diff --git a/openEMS/tools/global.cpp b/openEMS/tools/global.cpp
new file mode 100644
index 0000000..3e97fac
--- /dev/null
+++ b/openEMS/tools/global.cpp
@@ -0,0 +1,78 @@
+/*
+*	Copyright (C) 2010 Sebastian Held <sebastian.held@gmx.de>
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <cstring>
+#include <iostream>
+#include "global.h"
+
+using namespace std;
+
+// create global object
+Global g_settings;
+
+Global::Global()
+{
+	m_showProbeDiscretization = false;
+	m_nativeFieldDumps = false;
+	m_VerboseLevel = 0;
+}
+
+void Global::ShowArguments(ostream& ostr, string front)
+{
+	ostr << front << "--showProbeDiscretization\tShow probe discretization information" << endl;
+	ostr << front << "--nativeFieldDumps\t\tDump all fields using the native field components" << endl;
+	ostr << front << "-v,-vv,-vvv\t\t\tSet debug level: 1 to 3" << endl;
+}
+
+//! \brief This function initializes the object
+bool Global::parseCommandLineArgument( const char *argv )
+{
+	if (!argv)
+		return false;
+
+	if (strcmp(argv,"--showProbeDiscretization")==0)
+	{
+		cout << "openEMS - showing probe discretization information" << endl;
+		m_showProbeDiscretization = true;
+		return true;
+	}
+	else if (strcmp(argv,"--nativeFieldDumps")==0)
+	{
+		cout << "openEMS - dumping all fields using the native field components" << endl;
+		m_nativeFieldDumps = true;
+		return true;
+	}
+	else if (strcmp(argv,"-v")==0)
+	{
+		cout << "openEMS - verbose level 1" << endl;
+		m_VerboseLevel = 1;
+		return true;
+	}
+	else if (strcmp(argv,"-vv")==0)
+	{
+		cout << "openEMS - verbose level 2" << endl;
+		m_VerboseLevel = 2;
+		return true;
+	}
+	else if (strcmp(argv,"-vvv")==0)
+	{
+		cout << "openEMS - verbose level 3" << endl;
+		m_VerboseLevel = 3;
+		return true;
+	}
+	return false;
+}
diff --git a/openEMS/tools/global.h b/openEMS/tools/global.h
new file mode 100644
index 0000000..c3bb524
--- /dev/null
+++ b/openEMS/tools/global.h
@@ -0,0 +1,65 @@
+/*
+*	Copyright (C) 2010 Sebastian Held <sebastian.held@gmx.de>
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef GLOBAL_H
+#define GLOBAL_H
+
+#include <sstream>
+
+#define UNUSED(x) (void)(x);
+
+class Global
+{
+public:
+	Global();
+
+	//! Show all possible (global) command line arguments
+	void ShowArguments(std::ostream& ostr, std::string front=std::string());
+
+	//! Parse the given command line arguments
+	bool parseCommandLineArgument( const char *argv );
+
+	bool showProbeDiscretization() const {return m_showProbeDiscretization;}
+
+	//! Returns true if native field dumps are requested...
+	bool NativeFieldDumps() const {return m_nativeFieldDumps;}
+	//! Set dumps to use native fields.
+	void SetNativeFieldDumps(bool val) {m_nativeFieldDumps=val;}
+
+	//! Set the verbose level
+	void SetVerboseLevel(int level) {m_VerboseLevel=level;m_SavedVerboseLevel=level;}
+	//! Get the verbose level
+	int GetVerboseLevel() const {return m_VerboseLevel;}
+
+	//! Set a new verbose level temporarily, restore it with RestoreVerboseLevel()
+	void SetTempVerboseLevel(int level) {m_SavedVerboseLevel=m_VerboseLevel;m_VerboseLevel=level;}
+	//! Restore the temporarily overwritten verbose level
+	void RestoreVerboseLevel() {m_VerboseLevel=m_SavedVerboseLevel;}
+
+protected:
+	bool m_showProbeDiscretization;
+	bool m_nativeFieldDumps;
+	int m_VerboseLevel;
+	int m_SavedVerboseLevel;
+};
+
+extern Global g_settings;
+
+// declare a parameter as unused
+#define UNUSED(x) (void)(x);
+
+#endif // GLOBAL_H
diff --git a/openEMS/tools/hdf5_file_reader.cpp b/openEMS/tools/hdf5_file_reader.cpp
new file mode 100644
index 0000000..a6e1359
--- /dev/null
+++ b/openEMS/tools/hdf5_file_reader.cpp
@@ -0,0 +1,689 @@
+/*
+*	Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "hdf5_file_reader.h"
+#include "../tools/array_ops.h"
+#include <hdf5.h>
+
+#include <sstream>
+#include <iostream>
+#include <iomanip>
+
+using namespace std;
+
+HDF5_File_Reader::HDF5_File_Reader(string filename)
+{
+	m_filename = filename;
+	//suppress hdf5 error output
+	//H5Eset_auto(NULL, NULL);
+}
+
+HDF5_File_Reader::~HDF5_File_Reader()
+{
+}
+
+bool HDF5_File_Reader::IsValid()
+{
+	htri_t val = H5Fis_hdf5(m_filename.c_str());
+	if (val<0)
+	{
+		cerr << "HDF5_File_Reader::IsValid: the given file """ << m_filename << """ is not accessible..." << endl;
+		return false;
+	}
+	if (val==0)
+	{
+		cerr << "HDF5_File_Reader::IsValid: the given file """ << m_filename << """ is invalid..." << endl;
+		return false;
+	}
+	if (val==0)
+	cerr << "HDF5_File_Reader::IsValid: the given file """ << m_filename << """ is valid..." << endl;
+	return true;
+}
+
+bool HDF5_File_Reader::OpenGroup(hid_t &file, hid_t &group, string groupName)
+{
+	file = H5Fopen( m_filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT );
+	if (file==-1)
+	{
+		cerr << "HDF5_File_Reader::OpenGroup: opening the given file """ << m_filename << """ failed" << endl;
+		return 0;
+	}
+	if (H5Lexists(file, groupName.c_str(), H5P_DEFAULT)<=0)
+	{
+		H5Fclose(file);
+		return 0;
+	}
+
+	group = H5Gopen(file, groupName.c_str(), H5P_DEFAULT );
+	if (group<0)
+	{
+		cerr << "HDF5_File_Reader::OpenGroup: can't open group """ << groupName << """" << endl;
+		H5Fclose(file);
+		return 0;
+	}
+	return true;
+}
+
+bool HDF5_File_Reader::ReadAttribute(string grp_name, string attr_name, vector<float> &attr_values)
+{
+	vector<double> d_attr_values;
+	if (ReadAttribute(grp_name, attr_name, d_attr_values)==false)
+		return false;
+	attr_values.resize(d_attr_values.size(),0);
+	for (size_t n=0;n<d_attr_values.size();++n)
+		attr_values.at(n)=d_attr_values.at(n);
+	return true;
+}
+
+bool HDF5_File_Reader::ReadAttribute(string grp_name, string attr_name, vector<double> &attr_values)
+{
+	attr_values.clear();
+
+	hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT );
+	if (hdf5_file==-1)
+	{
+		cerr << "HDF5_File_Reader::OpenGroup: opening the given file """ << m_filename << """ failed" << endl;
+		return 0;
+	}
+
+	if (H5Lexists(hdf5_file, grp_name.c_str(), H5P_DEFAULT)<=0)
+	{
+		H5Fclose(hdf5_file);
+		return false;
+	}
+
+	hid_t attr = H5Aopen_by_name(hdf5_file, grp_name.c_str(), attr_name.c_str(), H5P_DEFAULT, H5P_DEFAULT);
+	if (attr==-1)
+	{
+		cerr << "HDF5_File_Reader::ReadAttribute: Opening the given Attribute: """ << attr_name << """ failed" << endl;
+		H5Fclose(hdf5_file);
+		return false;
+	}
+
+	hid_t type = H5Aget_type(attr);
+	if (type<0)
+	{
+		cerr << "HDF5_File_Reader::ReadAttribute: Dataset type error" << endl;
+		H5Aclose(attr);
+		H5Fclose(hdf5_file);
+		return false;
+	}
+
+	attr_values.clear();
+	if (H5Tget_class(type)==H5T_FLOAT)
+	{
+		size_t numVal = H5Aget_storage_size(attr)/H5Tget_size(type);
+		hid_t datatype=-1;
+		void *value=NULL;
+		float *f_value=NULL;
+		double *d_value=NULL;
+		if (H5Tget_size(type)==sizeof(float))
+		{
+			f_value = new float[numVal];
+			value = f_value;
+			datatype = H5T_NATIVE_FLOAT;
+		}
+		if (H5Tget_size(type)==sizeof(double))
+		{
+			d_value = new double[numVal];
+			value = d_value;
+			datatype = H5T_NATIVE_DOUBLE;
+		}
+		if (H5Aread(attr, datatype, value)<0)
+		{
+		  cerr << "HDF5_File_Reader::ReadAttribute: Reading the given Attribute failed" << endl;
+		  H5Aclose(attr);
+		  H5Fclose(hdf5_file);
+		  return false;
+		}
+		if (f_value)
+			for (size_t n=0;n<numVal;++n)
+				attr_values.push_back(f_value[n]);
+		if (d_value)
+			for (size_t n=0;n<numVal;++n)
+				attr_values.push_back(d_value[n]);
+		delete[] f_value;
+		delete[] d_value;
+	}
+	else
+	{
+		cerr << "HDF5_File_Reader::ReadAttribute: Attribute type not supported" << endl;
+		H5Aclose(attr);
+		H5Fclose(hdf5_file);
+		return false;
+	}
+	H5Aclose(attr);
+	H5Fclose(hdf5_file);
+	return true;
+}
+
+bool HDF5_File_Reader::ReadDataSet(string ds_name, hsize_t &nDim, hsize_t* &dims, float* &data)
+{
+	double* d_data;
+	if (ReadDataSet(ds_name, nDim, dims, d_data)==false)
+		return false;
+	hsize_t data_size = 1;
+	for (unsigned int d=0;d<nDim;++d)
+		data_size*=dims[d];
+	data = new float[data_size];
+	for (size_t n=0;n<data_size;++n)
+		data[n]=d_data[n];
+	delete[] d_data;
+	return true;
+}
+
+bool HDF5_File_Reader::ReadDataSet(string ds_name, hsize_t &nDim, hsize_t* &dims, double* &data)
+{
+	if (IsValid()==false)
+		return false;
+
+	hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT );
+	if (hdf5_file==-1)
+	{
+		cerr << "HDF5_File_Reader::ReadDataSet: opening the given file """ << m_filename << """ failed" << endl;
+		return false;
+	}
+
+	hid_t dataset = H5Dopen(hdf5_file, ds_name.c_str(), H5P_DEFAULT );
+	if (dataset<0)
+	{
+		cerr << "HDF5_File_Reader::ReadDataSet: dataset not found" << endl;
+		H5Fclose(hdf5_file);
+		return false;
+	}
+	hid_t type = H5Dget_type(dataset);
+	if (type<0)
+	{
+		cerr << "HDF5_File_Reader::ReadDataSet: dataset type error" << endl;
+		H5Dclose(dataset);
+		H5Fclose(hdf5_file);
+		return false;
+	}
+	if (H5Tget_class(type)!=H5T_FLOAT)
+	{
+		cerr << "HDF5_File_Reader::ReadDataSet: dataset type not a float" << endl;
+		H5Dclose(dataset);
+		H5Fclose(hdf5_file);
+		return false;
+	}
+
+	hid_t space = H5Dget_space(dataset);
+	nDim = H5Sget_simple_extent_ndims(space);
+	dims = new hsize_t[nDim];
+	H5Sget_simple_extent_dims(space, dims, NULL );
+	hsize_t data_size = 1;
+	for (unsigned int d=0;d<nDim;++d)
+		data_size*=dims[d];
+
+	void *value=NULL;
+	float *f_value=NULL;
+	data = new double[data_size];
+	if (H5Tget_size(type)==sizeof(float))
+	{
+		f_value = new float[data_size];
+		value = f_value;
+	}
+	else
+		value = data;
+
+	if (H5Dread(dataset,type,H5S_ALL,H5S_ALL,H5P_DEFAULT,value)<0)
+	{
+		cerr << "HDF5_File_Reader::ReadDataSet: error reading data" << endl;
+		H5Dclose(dataset);
+		H5Fclose(hdf5_file);
+		delete[] data;
+		delete[] f_value;
+		data=NULL;
+		return false;
+	}
+	if (f_value)
+		for (size_t n=0;n<data_size;++n)
+			data[n]=f_value[n];
+	delete[] f_value;
+	H5Dclose(dataset);
+	H5Fclose(hdf5_file);
+	return true;
+}
+
+bool HDF5_File_Reader::ReadMesh(float** lines, unsigned int* numLines, int &meshType)
+{
+	if (IsValid()==false)
+		return false;
+
+	hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT );
+	if (hdf5_file==-1)
+	{
+		cerr << "HDF5_File_Reader::ReadMesh: opening the given file """ << m_filename << """ failed" << endl;
+		return false;
+	}
+
+	vector<string> Mesh_Names;
+	if (H5Lexists(hdf5_file, "/Mesh/x", H5P_DEFAULT) && H5Lexists(hdf5_file, "/Mesh/y", H5P_DEFAULT) && H5Lexists(hdf5_file, "/Mesh/z", H5P_DEFAULT))
+	{
+		meshType = 0;
+		Mesh_Names.push_back("/Mesh/x");
+		Mesh_Names.push_back("/Mesh/y");
+		Mesh_Names.push_back("/Mesh/z");
+	}
+	else if (H5Lexists(hdf5_file, "/Mesh/rho", H5P_DEFAULT) && H5Lexists(hdf5_file, "/Mesh/alpha", H5P_DEFAULT) && H5Lexists(hdf5_file, "/Mesh/z", H5P_DEFAULT))
+	{
+		meshType = 1;
+		Mesh_Names.push_back("/Mesh/rho");
+		Mesh_Names.push_back("/Mesh/alpha");
+		Mesh_Names.push_back("/Mesh/z");
+	}
+	else
+	{
+		cerr << "HDF5_File_Reader::ReadMesh: no falid mesh information found" << endl;
+		H5Fclose(hdf5_file);
+		return false;
+	}
+
+	for (int n=0;n<3;++n)
+	{
+		hsize_t nDim;
+		hsize_t* dims=NULL;
+		float* data=NULL;
+		ReadDataSet(Mesh_Names.at(n), nDim, dims, data);
+		if (nDim!=1)
+		{
+			cerr << "HDF5_File_Reader::ReadMesh: mesh dimension error" << endl;
+			delete[] dims;
+			delete[] data;
+			H5Fclose(hdf5_file);
+			return false;
+		}
+		numLines[n]=dims[0];
+		delete[] dims;
+		lines[n]=data;
+	}
+	H5Fclose(hdf5_file);
+	return true;
+}
+
+unsigned int HDF5_File_Reader::GetNumTimeSteps()
+{
+	if (IsValid()==false)
+		return false;
+
+	hid_t hdf5_file;
+	hid_t TD_grp;
+	if (OpenGroup(hdf5_file, TD_grp, "/FieldData/TD")==false)
+		return false;
+
+	hsize_t numObj;
+	if (H5Gget_num_objs(TD_grp,&numObj)<0)
+	{
+		cerr << "HDF5_File_Reader::GetNumTimeSteps: can't read number of datasets" << endl;
+		H5Gclose(TD_grp);
+		H5Fclose(hdf5_file);
+		return 0;
+	}
+	H5Gclose(TD_grp);
+	H5Fclose(hdf5_file);
+	return numObj;
+}
+
+bool HDF5_File_Reader::ReadTimeSteps(vector<unsigned int> &timestep, vector<string> &names)
+{
+	if (IsValid()==false)
+		return false;
+
+	hid_t hdf5_file;
+	hid_t TD_grp;
+	if (OpenGroup(hdf5_file, TD_grp, "/FieldData/TD")==false)
+		return false;
+
+	hsize_t numObj;
+	if (H5Gget_num_objs(TD_grp,&numObj)<0)
+	{
+		cerr << "HDF5_File_Reader::ReadTimeSteps: can't read number of datasets" << endl;
+		H5Gclose(TD_grp);
+		H5Fclose(hdf5_file);
+		return false;
+	}
+	char name[100];
+	timestep.clear();
+	timestep.resize(numObj,0);
+	names.clear();
+	names.resize(numObj);
+	for (hsize_t n=0;n<numObj;++n)
+	{
+		if (H5Gget_objtype_by_idx(TD_grp, n)  != H5G_DATASET)
+		{
+			cerr << "HDF5_File_Reader::ReadTimeSteps: invalid timestep found!" << endl;
+			H5Gclose(TD_grp);
+			H5Fclose(hdf5_file);
+			return false;
+		}
+		H5Gget_objname_by_idx(TD_grp, n, name, 100 );
+		istringstream is(name);
+		unsigned int num;
+		if (is >> num)
+		{
+			timestep.at(n)=num;
+			names.at(n)=name;
+		}
+		else
+		{
+			cerr << "HDF5_File_Reader::ReadTimeSteps: invalid timestep format found!" << endl;
+			H5Gclose(TD_grp);
+			H5Fclose(hdf5_file);
+			return false;
+		}
+	}
+	H5Gclose(TD_grp);
+	H5Fclose(hdf5_file);
+	return true;
+}
+
+float**** HDF5_File_Reader::GetTDVectorData(size_t idx, float &time, unsigned int data_size[])
+{
+	if (IsValid()==false)
+		return NULL;
+
+	hid_t hdf5_file;
+	hid_t TD_grp;
+	if (OpenGroup(hdf5_file, TD_grp, "/FieldData/TD")==false)
+		return NULL;
+
+	hsize_t numObj;
+	if (H5Gget_num_objs(TD_grp,&numObj)<0)
+	{
+		cerr << "HDF5_File_Reader::GetTDVectorData: can't read number of datasets" << endl;
+		H5Gclose(TD_grp);
+		H5Fclose(hdf5_file);
+		return NULL;
+	}
+	if (idx>=numObj)
+	{
+		H5Gclose(TD_grp);
+		H5Fclose(hdf5_file);
+		return NULL;
+	}
+
+	if (H5Gget_objtype_by_idx(TD_grp, idx)  != H5G_DATASET)
+	{
+		cerr << "HDF5_File_Reader::GetTDVectorData: invalid timestep found!" << endl;
+		H5Gclose(TD_grp);
+		H5Fclose(hdf5_file);
+		return NULL;
+	}
+
+	char name[100];
+	H5Gget_objname_by_idx(TD_grp, idx, name, 100 );
+	string ds_name = "/FieldData/TD/" + string(name);
+
+	hid_t attr = H5Aopen_by_name(hdf5_file, ds_name.c_str(), "time", H5P_DEFAULT, H5P_DEFAULT);
+	if (attr<0)
+	{
+		cerr << "HDF5_File_Reader::GetTDVectorData: time attribute not found!" << endl;
+		H5Gclose(TD_grp);
+		H5Fclose(hdf5_file);
+		return NULL;
+	}
+	if (H5Aread(attr, H5T_NATIVE_FLOAT, &time)<0)
+	{
+		cerr << "HDF5_File_Reader::GetTDVectorData: can't read time attribute!" << endl;
+		H5Aclose(attr);
+		H5Gclose(TD_grp);
+		H5Fclose(hdf5_file);
+		return NULL;
+	}
+
+	hsize_t nDim;
+	hsize_t* dims=NULL;
+	double* data=NULL;
+	ReadDataSet(ds_name, nDim, dims, data);
+	if (nDim!=4)
+	{
+		cerr << "HDF5_File_Reader::GetTDVectorData: data dimension invalid" << endl;
+		delete[] dims;
+		H5Aclose(attr);
+		H5Gclose(TD_grp);
+		H5Fclose(hdf5_file);
+		return NULL;
+	}
+	if (dims[0]!=3)
+	{
+		cerr << "HDF5_File_Reader::GetTDVectorData: vector data dimension invalid" << endl;
+		delete[] dims;
+		H5Aclose(attr);
+		H5Gclose(TD_grp);
+		H5Fclose(hdf5_file);
+		return NULL;
+	}
+	data_size[0]=dims[3];
+	data_size[1]=dims[2];
+	data_size[2]=dims[1];
+	delete[] dims;
+	data_size[3]=3;
+	size_t pos = 0;
+	float**** field = Create_N_3DArray<float>(data_size);
+	for (unsigned int d=0;d<3;++d)
+		for (unsigned int k=0;k<data_size[2];++k)
+			for (unsigned int j=0;j<data_size[1];++j)
+				for (unsigned int i=0;i<data_size[0];++i)
+				{
+					field[d][i][j][k]=data[pos++];
+				}
+	delete[] data;
+	H5Aclose(attr);
+	H5Gclose(TD_grp);
+	H5Fclose(hdf5_file);
+	return field;
+}
+
+unsigned int HDF5_File_Reader::GetNumFrequencies()
+{
+	vector<float> frequencies;
+	if (ReadFrequencies(frequencies)==false)
+		return 0;
+	return frequencies.size();
+}
+
+bool HDF5_File_Reader::ReadFrequencies(vector<float> &frequencies)
+{
+	if (IsValid()==false)
+		return false;
+
+	return ReadAttribute("/FieldData/FD","frequency",frequencies);
+}
+
+bool HDF5_File_Reader::ReadFrequencies(vector<double> &frequencies)
+{
+	if (IsValid()==false)
+		return false;
+
+	return ReadAttribute("/FieldData/FD","frequency",frequencies);
+}
+
+
+complex<float>**** HDF5_File_Reader::GetFDVectorData(size_t idx, unsigned int data_size[])
+{
+	hsize_t nDim;
+	hsize_t* dims=NULL;
+	double* data=NULL;
+	stringstream ds_name;
+
+	// read real values
+	ds_name << "/FieldData/FD/f" << idx << "_real";
+	if (ReadDataSet(ds_name.str(), nDim, dims, data) == false)
+		return NULL;
+	if (nDim!=4)
+	{
+		cerr << "HDF5_File_Reader::GetFDVectorData: data dimension invalid" << endl;
+		delete[] dims;
+		delete[] data;
+		return NULL;
+	}
+	if (dims[0]!=3)
+	{
+		cerr << "HDF5_File_Reader::GetFDVectorData: vector data dimension invalid" << endl;
+		delete[] dims;
+		delete[] data;
+		return NULL;
+	}
+	data_size[0]=dims[3];
+	data_size[1]=dims[2];
+	data_size[2]=dims[1];
+	delete[] dims;
+	data_size[3]=3;
+	size_t pos = 0;
+	complex<float>**** field = Create_N_3DArray<complex<float> >(data_size);
+	for (unsigned int d=0;d<3;++d)
+		for (unsigned int k=0;k<data_size[2];++k)
+			for (unsigned int j=0;j<data_size[1];++j)
+				for (unsigned int i=0;i<data_size[0];++i)
+				{
+					field[d][i][j][k]=data[pos++];
+				}
+	delete[] data;
+
+	// read imaginary values
+	ds_name.str("");
+	ds_name << "/FieldData/FD/f" << idx << "_imag";
+	if (ReadDataSet(ds_name.str(), nDim, dims, data) == false)
+	{
+		Delete_N_3DArray<complex<float> >(field, data_size);
+		return NULL;
+	}
+	if (nDim!=4)
+	{
+		cerr << "HDF5_File_Reader::GetFDVectorData: data dimension invalid" << endl;
+		delete[] dims;
+		delete[] data;
+		Delete_N_3DArray<complex<float> >(field, data_size);
+		return NULL;
+	}
+	if (dims[0]!=3)
+	{
+		cerr << "HDF5_File_Reader::GetFDVectorData: vector data dimension invalid" << endl;
+		delete[] dims;
+		delete[] data;
+		Delete_N_3DArray<complex<float> >(field, data_size);
+		return NULL;
+	}
+	if ((data_size[0]!=dims[3]) || (data_size[1]!=dims[2]) || (data_size[2]!=dims[1]))
+	{
+		cerr << "HDF5_File_Reader::GetFDVectorData: data dimension mismatch" << endl;
+		delete[] dims;
+		delete[] data;
+		Delete_N_3DArray<complex<float> >(field, data_size);
+		return NULL;
+	}
+	delete[] dims;
+
+	pos = 0;
+	complex<double> I(0,1);
+	for (unsigned int d=0;d<3;++d)
+		for (unsigned int k=0;k<data_size[2];++k)
+			for (unsigned int j=0;j<data_size[1];++j)
+				for (unsigned int i=0;i<data_size[0];++i)
+				{
+					field[d][i][j][k]+= I*data[pos++];
+				}
+	delete[] data;
+
+	return field;
+}
+
+bool HDF5_File_Reader::CalcFDVectorData(vector<float> &frequencies, vector<complex<float>****> &FD_data, unsigned int data_size[4])
+{
+	FD_data.clear();
+
+	if (GetNumTimeSteps()<=0)
+	{
+		cerr << "HDF5_File_Reader::CalcFDVectorData: error, no TD data found..." << endl;
+		return false;
+	}
+
+	float time;
+	//read first TD data
+	float**** field = this->GetTDVectorData(0,time,data_size);
+	if (field==NULL)
+	{
+		cerr << "HDF5_File_Reader::CalcFDVectorData: error, no TD data found..." << endl;
+		return false;
+	}
+
+	//init
+	FD_data.resize(frequencies.size(), NULL);
+	for (size_t fn=0;fn<frequencies.size();++fn)
+		FD_data.at(fn) = Create_N_3DArray<complex<float> >(data_size);
+
+	size_t ts=0;
+	unsigned int pos[3];
+	complex<float> PI_2_I(0.0,-2.0*M_PI);
+	complex<float> exp_jwt_2_dt;
+	float time_diff=0;
+	float time_old =0;
+	complex<float>**** field_fd = NULL;
+	while (field)
+	{
+		if ((ts>1) && abs(time_diff - (time - time_old))>1e15)
+		{
+			cerr << "HDF5_File_Reader::CalcFDVectorData: time interval error..." << endl;
+			for (size_t fn=0;fn<frequencies.size();++fn)
+				Delete_N_3DArray(FD_data.at(fn),data_size);
+			FD_data.clear();
+			return false;
+		}
+		time_diff = time - time_old;
+		for (size_t fn=0;fn<frequencies.size();++fn)
+		{
+			exp_jwt_2_dt = exp( (complex<float>)(PI_2_I * frequencies.at(fn) * time) );
+			field_fd = FD_data.at(fn);
+			for (pos[0]=0; pos[0]<data_size[0]; ++pos[0])
+			{
+				for (pos[1]=0; pos[1]<data_size[1]; ++pos[1])
+				{
+					for (pos[2]=0; pos[2]<data_size[2]; ++pos[2])
+					{
+						field_fd[0][pos[0]][pos[1]][pos[2]] += field[0][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
+						field_fd[1][pos[0]][pos[1]][pos[2]] += field[1][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
+						field_fd[2][pos[0]][pos[1]][pos[2]] += field[2][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
+					}
+				}
+			}
+		}
+		++ts;
+		Delete_N_3DArray(field,data_size);
+		time_old = time;
+		field = this->GetTDVectorData(ts,time,data_size);
+	}
+
+	// finalize data
+	time_diff*=2;
+	for (size_t fn=0;fn<frequencies.size();++fn)
+	{
+		field_fd = FD_data.at(fn);
+		for (pos[0]=0; pos[0]<data_size[0]; ++pos[0])
+		{
+			for (pos[1]=0; pos[1]<data_size[1]; ++pos[1])
+			{
+				for (pos[2]=0; pos[2]<data_size[2]; ++pos[2])
+				{
+					field_fd[0][pos[0]][pos[1]][pos[2]] *= time_diff;
+					field_fd[1][pos[0]][pos[1]][pos[2]] *= time_diff;
+					field_fd[2][pos[0]][pos[1]][pos[2]] *= time_diff;
+				}
+			}
+		}
+	}
+	return true;
+}
diff --git a/openEMS/tools/hdf5_file_reader.h b/openEMS/tools/hdf5_file_reader.h
new file mode 100644
index 0000000..f190240
--- /dev/null
+++ b/openEMS/tools/hdf5_file_reader.h
@@ -0,0 +1,78 @@
+/*
+*	Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef HDF5_FILE_READER_H
+#define HDF5_FILE_READER_H
+
+#include <string>
+#include <vector>
+#include <complex>
+#include <hdf5.h>
+
+class HDF5_File_Reader
+{
+public:
+	HDF5_File_Reader(std::string filename);
+	virtual ~HDF5_File_Reader();
+
+	bool ReadMesh(float** lines, unsigned int* numLines, int &meshType);
+
+	//! Get the number of timesteps stored at /FieldData/TD/<NUMBER_OF_TS>
+	unsigned int GetNumTimeSteps();
+	bool ReadTimeSteps(std::vector<unsigned int> &timestep, std::vector<std::string> &names);
+
+	/*!
+	  Get time-domain data stored at /FieldData/TD/<NUMBER_OF_TS>
+	  \param[in]  idx	time step index to extract
+	  \param[out] time	time attribute for the given timestep
+	  \param[out] data_size data size found
+	  \return field data found in given timestep, caller must delete array, returns NULL if timestep was not found
+	  */
+	float**** GetTDVectorData(size_t idx, float &time, unsigned int data_size[4]);
+
+	unsigned int GetNumFrequencies();
+	bool ReadFrequencies(std::vector<float> &frequencies);
+	bool ReadFrequencies(std::vector<double> &frequencies);
+
+	/*!
+	  Get frequency-domain data stored at "/FieldData/FD/f<idx>_real" and "/FieldData/FD/f<idx>_imag"
+	  \param[in]  idx	frequency index to extract
+	  \param[out] data_size data size found
+	  \return complex field data found for the given frequency, caller must delete array, returns NULL if frequency was not found
+	  */
+	std::complex<float>**** GetFDVectorData(size_t idx, unsigned int data_size[4]);
+
+	/*!
+	  Calculate
+	  */
+	bool CalcFDVectorData(std::vector<float> &frequencies, std::vector<std::complex<float>****> &FD_data, unsigned int data_size[4]);
+
+	bool ReadAttribute(std::string grp_name, std::string attr_name, std::vector<double> &attr_values);
+	bool ReadAttribute(std::string grp_name, std::string attr_name, std::vector<float> &attr_values);
+
+	bool IsValid();
+
+protected:
+	std::string m_filename;
+
+	bool ReadDataSet(std::string ds_name, hsize_t &nDim, hsize_t* &dims, double* &data);
+	bool ReadDataSet(std::string ds_name, hsize_t &nDim, hsize_t* &dims, float* &data);
+
+	bool OpenGroup(hid_t &file, hid_t &group, std::string groupName);
+};
+
+#endif // HDF5_FILE_READER_H
diff --git a/openEMS/tools/hdf5_file_writer.cpp b/openEMS/tools/hdf5_file_writer.cpp
new file mode 100644
index 0000000..e3a2113
--- /dev/null
+++ b/openEMS/tools/hdf5_file_writer.cpp
@@ -0,0 +1,515 @@
+/*
+*	Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+using namespace std;
+
+#include "hdf5_file_writer.h"
+#include <boost/algorithm/string.hpp>
+#include <hdf5.h>
+
+#include <sstream>
+#include <iostream>
+#include <iomanip>
+
+HDF5_File_Writer::HDF5_File_Writer(string filename)
+{
+	m_filename = filename;
+	m_Group = "/";
+	hid_t hdf5_file = H5Fcreate(m_filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+	if (hdf5_file<0)
+	{
+		cerr << "HDF5_File_Writer::HDF5_File_Writer: Error, creating the given file """ << m_filename << """ failed" << endl;
+	}
+	H5Fclose(hdf5_file);
+}
+
+HDF5_File_Writer::~HDF5_File_Writer()
+{
+}
+
+hid_t HDF5_File_Writer::OpenGroup(hid_t hdf5_file, string group)
+{
+	if (hdf5_file<0)
+	{
+		cerr << "HDF5_File_Writer::CreateGroup: Error, invalid file id" << endl;
+		return -1;
+	}
+
+	vector<string> results;
+	boost::split(results, group, boost::is_any_of("/"));
+
+	hid_t grp=H5Gopen(hdf5_file,"/", H5P_DEFAULT);
+	if (grp<0)
+	{
+		cerr << "HDF5_File_Writer::OpenGroup: Error, opening root group " << endl;
+		return -1;
+	}
+
+	for (size_t n=0;n<results.size();++n)
+	{
+		hid_t old_grp = grp;
+		if (!results.at(n).empty())
+		{
+			if (H5Lexists(grp, results.at(n).c_str(), H5P_DEFAULT))
+			{
+			  grp = H5Gopen(grp, results.at(n).c_str(), H5P_DEFAULT);
+				H5Gclose(old_grp);
+				if (grp<0)
+				{
+					cerr << "HDF5_File_Writer::OpenGroup: Error, failed to open existing group" << endl;
+					return -1;
+				}
+			}
+			else
+			{
+			  grp = H5Gcreate(grp,results.at(n).c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+				H5Gclose(old_grp);
+				if (grp<0)
+				{
+					cerr << "HDF5_File_Writer::OpenGroup: Error, creating group """ << group << """ failed" << endl;
+					return -1;
+				}
+			}
+		}
+	}
+	return grp;
+}
+
+void HDF5_File_Writer::SetCurrentGroup(std::string group, bool createGrp)
+{
+	m_Group = group;
+	if (createGrp==false)
+		return;
+
+	hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT );
+	if (hdf5_file<0)
+	{
+		cerr << "HDF5_File_Writer::SetCurrentGroup: Error, opening the given file """ << m_filename << """ failed" << endl;
+		return;
+	}
+	hid_t hdf5_group = OpenGroup(hdf5_file, m_Group);
+	if (hdf5_group<0)
+		cerr << "HDF5_File_Writer::WriteData: Error opening group" << endl;
+	H5Gclose(hdf5_group);
+	H5Fclose(hdf5_file);
+}
+
+bool HDF5_File_Writer::WriteRectMesh(unsigned int const* numLines, double const* const* discLines, int MeshType, double scaling)
+{
+	float* array[3];
+	for (int n=0; n<3; ++n)
+	{
+		array[n] = new float[numLines[n]];
+		for (unsigned int i=0; i<numLines[n]; ++i)
+			array[n][i]=discLines[n][i];
+
+	}
+	bool success = WriteRectMesh(numLines,array,MeshType,scaling);
+	for (int n=0; n<3; ++n)
+		delete[] array[n];
+	return success;
+}
+
+bool HDF5_File_Writer::WriteRectMesh(unsigned int const* numLines, float const* const* discLines, int MeshType, float scaling)
+{
+	hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT );
+	if (hdf5_file<0)
+	{
+		cerr << "HDF5_File_Writer::WriteRectMesh: Error, opening the given file """ << m_filename << """ failed" << endl;
+		return false;
+	}
+
+	if (H5Lexists(hdf5_file, "/Mesh", H5P_DEFAULT))
+	{
+		cerr << "HDF5_File_Writer::WriteRectMesh: Error, group ""/Mesh"" already exists" << endl;
+		H5Fclose(hdf5_file);
+		return false;
+	}
+
+	hid_t mesh_grp = H5Gcreate(hdf5_file,"/Mesh", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+	if (mesh_grp<0)
+	{
+		cerr << "HDF5_File_Writer::WriteRectMesh: Error, creating group ""/Mesh"" failed" << endl;
+		H5Fclose(hdf5_file);
+		return false;
+	}
+
+	string names[] = {"x","y","z"};
+	if (MeshType==1)
+	{
+		names[0]="rho";
+		names[1]="alpha";
+	}
+	if (MeshType==2)
+	{
+		names[0]="r";
+		names[1]="theta";
+		names[2]="phi";
+	}
+
+	for (int n=0; n<3; ++n)
+	{
+		hsize_t dims[1]={numLines[n]};
+		hid_t space = H5Screate_simple(1, dims, NULL);
+		hid_t dataset = H5Dcreate(mesh_grp, names[n].c_str(), H5T_NATIVE_FLOAT, space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+		float* array = new float[numLines[n]];
+		for (unsigned int i=0; i<numLines[n]; ++i)
+		{
+			if ((MeshType==1) && (n==1)) //check for alpha-direction
+				array[i] = discLines[n][i];
+			else if ((MeshType==2) && (n>0)) //check for theta/phi-direction
+				array[i] = discLines[n][i];
+			else
+				array[i] = discLines[n][i] * scaling;
+		}
+		if (H5Dwrite(dataset, H5T_NATIVE_FLOAT, space, H5P_DEFAULT, H5P_DEFAULT, array))
+		{
+			cerr << "HDF5_File_Writer::WriteRectMesh: Error, writing to dataset failed" << endl;
+			delete[] array;
+			H5Dclose(dataset);
+			H5Sclose(space);
+			H5Gclose(mesh_grp);
+			H5Fclose(hdf5_file);
+			return false;
+		}
+		delete[] array;
+		H5Dclose(dataset);
+		H5Sclose(space);
+	}
+	H5Gclose(mesh_grp);
+	H5Fclose(hdf5_file);
+	return true;
+}
+
+bool HDF5_File_Writer::WriteScalarField(std::string dataSetName, float const* const* const* field, size_t datasize[3])
+{
+	size_t pos = 0;
+	size_t size = datasize[0]*datasize[1]*datasize[2];
+	size_t n_size[3]={datasize[2],datasize[1],datasize[0]};
+	float* buffer = new float[size];
+	for (size_t k=0;k<datasize[2];++k)
+		for (size_t j=0;j<datasize[1];++j)
+			for (size_t i=0;i<datasize[0];++i)
+			{
+				buffer[pos++]=field[i][j][k];
+			}
+	bool success = WriteData(dataSetName,buffer,3,n_size);
+	delete[] buffer;
+	return success;
+}
+
+bool HDF5_File_Writer::WriteScalarField(std::string dataSetName, double const* const* const* field, size_t datasize[3])
+{
+	size_t pos = 0;
+	size_t size = datasize[0]*datasize[1]*datasize[2];
+	size_t n_size[3]={datasize[2],datasize[1],datasize[0]};
+	double* buffer = new double[size];
+	for (size_t k=0;k<datasize[2];++k)
+		for (size_t j=0;j<datasize[1];++j)
+			for (size_t i=0;i<datasize[0];++i)
+			{
+				buffer[pos++]=field[i][j][k];
+			}
+	bool success = WriteData(dataSetName,buffer,3,n_size);
+	delete[] buffer;
+	return success;
+}
+
+bool HDF5_File_Writer::WriteScalarField(std::string dataSetName, std::complex<float> const* const* const* field, size_t datasize[3])
+{
+	size_t pos = 0;
+	size_t size = datasize[0]*datasize[1]*datasize[2];
+	size_t n_size[3]={datasize[2],datasize[1],datasize[0]};
+	float* buffer = new float[size];
+	for (size_t k=0;k<datasize[2];++k)
+		for (size_t j=0;j<datasize[1];++j)
+			for (size_t i=0;i<datasize[0];++i)
+			{
+				buffer[pos++]=real(field[i][j][k]);
+			}
+	bool success = WriteData(dataSetName + "_real",buffer,3,n_size);
+
+	pos = 0;
+	for (size_t k=0;k<datasize[2];++k)
+		for (size_t j=0;j<datasize[1];++j)
+			for (size_t i=0;i<datasize[0];++i)
+			{
+				buffer[pos++]=imag(field[i][j][k]);
+			}
+	success &= WriteData(dataSetName + "_imag",buffer,3,n_size);
+
+	delete[] buffer;
+	return success;
+}
+
+bool HDF5_File_Writer::WriteScalarField(std::string dataSetName, std::complex<double> const* const* const* field, size_t datasize[3])
+{
+	size_t pos = 0;
+	size_t size = datasize[0]*datasize[1]*datasize[2];
+	size_t n_size[3]={datasize[2],datasize[1],datasize[0]};
+	double* buffer = new double[size];
+	for (size_t k=0;k<datasize[2];++k)
+		for (size_t j=0;j<datasize[1];++j)
+			for (size_t i=0;i<datasize[0];++i)
+			{
+				buffer[pos++]=real(field[i][j][k]);
+			}
+	bool success = WriteData(dataSetName + "_real",buffer,3,n_size);
+
+	pos = 0;
+	for (size_t k=0;k<datasize[2];++k)
+		for (size_t j=0;j<datasize[1];++j)
+			for (size_t i=0;i<datasize[0];++i)
+			{
+				buffer[pos++]=imag(field[i][j][k]);
+			}
+	success &= WriteData(dataSetName + "_imag",buffer,3,n_size);
+
+	delete[] buffer;
+	return success;
+}
+
+bool HDF5_File_Writer::WriteVectorField(std::string dataSetName, float const* const* const* const* field, size_t datasize[3])
+{
+	size_t pos = 0;
+	size_t size = datasize[0]*datasize[1]*datasize[2]*3;
+	float* buffer = new float[size];
+	for (int n=0;n<3;++n)
+		for (size_t k=0;k<datasize[2];++k)
+			for (size_t j=0;j<datasize[1];++j)
+				for (size_t i=0;i<datasize[0];++i)
+				{
+					buffer[pos++]=field[n][i][j][k];
+				}
+	size_t n_size[4]={3,datasize[2],datasize[1],datasize[0]};
+	bool success = WriteData(dataSetName,buffer,4,n_size);
+	delete[] buffer;
+	return success;
+}
+
+bool HDF5_File_Writer::WriteVectorField(std::string dataSetName, double const* const* const* const* field, size_t datasize[3])
+{
+	size_t pos = 0;
+	size_t size = datasize[0]*datasize[1]*datasize[2]*3;
+	double* buffer = new double[size];
+	for (int n=0;n<3;++n)
+		for (size_t k=0;k<datasize[2];++k)
+			for (size_t j=0;j<datasize[1];++j)
+				for (size_t i=0;i<datasize[0];++i)
+				{
+					buffer[pos++]=field[n][i][j][k];
+				}
+	size_t n_size[4]={3,datasize[2],datasize[1],datasize[0]};
+	bool success = WriteData(dataSetName,buffer,4,n_size);
+	delete[] buffer;
+	return success;
+}
+
+bool HDF5_File_Writer::WriteVectorField(std::string dataSetName, std::complex<float> const* const* const* const* field, size_t datasize[3])
+{
+	size_t pos = 0;
+	size_t size = datasize[0]*datasize[1]*datasize[2]*3;
+	size_t n_size[4]={3,datasize[2],datasize[1],datasize[0]};
+	float* buffer = new float[size];
+	for (int n=0;n<3;++n)
+		for (size_t k=0;k<datasize[2];++k)
+			for (size_t j=0;j<datasize[1];++j)
+				for (size_t i=0;i<datasize[0];++i)
+				{
+					buffer[pos++]=real(field[n][i][j][k]);
+				}
+	bool success = WriteData(dataSetName + "_real",buffer,4,n_size);
+
+	pos = 0;
+	for (int n=0;n<3;++n)
+		for (size_t k=0;k<datasize[2];++k)
+			for (size_t j=0;j<datasize[1];++j)
+				for (size_t i=0;i<datasize[0];++i)
+				{
+					buffer[pos++]=imag(field[n][i][j][k]);
+				}
+	success &= WriteData(dataSetName + "_imag",buffer,4,n_size);
+
+	delete[] buffer;
+	return success;
+}
+
+bool HDF5_File_Writer::WriteVectorField(std::string dataSetName, std::complex<double> const* const* const* const* field, size_t datasize[3])
+{
+	size_t pos = 0;
+	size_t size = datasize[0]*datasize[1]*datasize[2]*3;
+	size_t n_size[4]={3,datasize[2],datasize[1],datasize[0]};
+	double* buffer = new double[size];
+	for (int n=0;n<3;++n)
+		for (size_t k=0;k<datasize[2];++k)
+			for (size_t j=0;j<datasize[1];++j)
+				for (size_t i=0;i<datasize[0];++i)
+				{
+					buffer[pos++]=real(field[n][i][j][k]);
+				}
+	bool success = WriteData(dataSetName + "_real",buffer,4,n_size);
+
+	pos = 0;
+	for (int n=0;n<3;++n)
+		for (size_t k=0;k<datasize[2];++k)
+			for (size_t j=0;j<datasize[1];++j)
+				for (size_t i=0;i<datasize[0];++i)
+				{
+					buffer[pos++]=imag(field[n][i][j][k]);
+				}
+	success &= WriteData(dataSetName + "_imag",buffer,4,n_size);
+
+	delete[] buffer;
+	return success;
+}
+
+bool HDF5_File_Writer::WriteData(std::string dataSetName, float const* field_buf, size_t dim, size_t* datasize)
+{
+	return WriteData(dataSetName, H5T_NATIVE_FLOAT, field_buf,dim, datasize);
+}
+
+bool HDF5_File_Writer::WriteData(std::string dataSetName, double const* field_buf, size_t dim, size_t* datasize)
+{
+	return WriteData(dataSetName, H5T_NATIVE_DOUBLE, field_buf,dim, datasize);
+}
+
+bool HDF5_File_Writer::WriteData(std::string dataSetName,  hid_t mem_type, void const* field_buf, size_t dim, size_t* datasize)
+{
+	hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT );
+	if (hdf5_file<0)
+	{
+		cerr << "HDF5_File_Writer::WriteData: Error, opening the given file """ << m_filename << """ failed" << endl;
+		return false;
+	}
+
+	hid_t group = OpenGroup(hdf5_file,m_Group);
+	if (group<0)
+	{
+		cerr << "HDF5_File_Writer::WriteData: Error opening group" << endl;
+		H5Fclose(hdf5_file);
+		return false;
+	}
+
+	hsize_t dims[dim];
+	for (size_t n=0;n<dim;++n)
+		dims[n]=datasize[n];
+	hid_t space = H5Screate_simple(dim, dims, NULL);
+	hid_t dataset = H5Dcreate(group, dataSetName.c_str(), mem_type, space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+	if (H5Dwrite(dataset, mem_type, space, H5P_DEFAULT, H5P_DEFAULT, field_buf))
+	{
+		cerr << "HDF5_File_Writer::WriteData: Error, writing to dataset failed" << endl;
+		H5Dclose(dataset);
+		H5Sclose(space);
+		H5Gclose(group);
+		H5Fclose(hdf5_file);
+		return false;
+	}
+	H5Dclose(dataset);
+	H5Sclose(space);
+	H5Gclose(group);
+	H5Fclose(hdf5_file);
+	return true;
+}
+
+bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name, void const* value, hsize_t size, hid_t mem_type)
+{
+	hid_t hdf5_file = H5Fopen( m_filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT );
+	if (hdf5_file<0)
+	{
+		cerr << "HDF5_File_Writer::WriteAtrribute: Error, opening the given file """ << m_filename << """ failed" << endl;
+		return false;
+	}
+
+	if (H5Lexists(hdf5_file, locName.c_str(), H5P_DEFAULT)<0)
+	{
+		cerr << "HDF5_File_Writer::WriteAtrribute: Error, failed to find location: """ << locName << """" << endl;
+		H5Fclose(hdf5_file);
+		return false;
+	}
+	hid_t loc = H5Oopen(hdf5_file, locName.c_str(), H5P_DEFAULT);
+	if (loc<0)
+	{
+		cerr << "HDF5_File_Writer::WriteAtrribute: Error, failed to open location: """ << locName << """" << endl;
+		H5Fclose(hdf5_file);
+		return false;
+	}
+
+	hid_t dataspace_id = H5Screate_simple(1, &size, NULL);
+
+	/* Create a dataset attribute. */
+	hid_t attribute_id = H5Acreate(loc, attr_name.c_str(), mem_type, dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
+	if (attribute_id<0)
+	{
+		cerr << "HDF5_File_Writer::WriteAtrribute: Error, failed to create the attrbute" << endl;
+		H5Sclose(dataspace_id);
+		H5Oclose(loc);
+		H5Fclose(hdf5_file);
+		return false;
+	}
+
+	/* Write the attribute data. */
+	if (H5Awrite(attribute_id, mem_type, value)<0)
+	{
+		cerr << "HDF5_File_Writer::WriteAtrribute: Error, failed to write the attrbute" << endl;
+		H5Aclose(attribute_id);
+		H5Sclose(dataspace_id);
+		H5Oclose(loc);
+		H5Fclose(hdf5_file);
+		return false;
+	}
+	H5Aclose(attribute_id);
+	H5Sclose(dataspace_id);
+	H5Oclose(loc);
+	H5Fclose(hdf5_file);
+	return true;
+}
+
+bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name, float const* value, hsize_t size)
+{
+	return WriteAtrribute(locName,attr_name, value, size, H5T_NATIVE_FLOAT);
+}
+
+bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name, double const* value, hsize_t size)
+{
+	return WriteAtrribute(locName,attr_name, value, size, H5T_NATIVE_DOUBLE);
+}
+
+bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name, vector<float> values)
+{
+	float val[values.size()];
+	for (size_t n=0;n<values.size();++n)
+		val[n]=values.at(n);
+	return HDF5_File_Writer::WriteAtrribute(locName, attr_name,val,values.size(),H5T_NATIVE_FLOAT);
+}
+
+bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name, vector<double> values)
+{
+	double val[values.size()];
+	for (size_t n=0;n<values.size();++n)
+		val[n]=values.at(n);
+	return HDF5_File_Writer::WriteAtrribute(locName, attr_name, val, values.size(), H5T_NATIVE_DOUBLE);
+}
+
+bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name, float value)
+{
+	return HDF5_File_Writer::WriteAtrribute(locName, attr_name,&value,1, H5T_NATIVE_FLOAT);
+}
+
+bool HDF5_File_Writer::WriteAtrribute(std::string locName, std::string attr_name, double value)
+{
+	return HDF5_File_Writer::WriteAtrribute(locName, attr_name,&value,1, H5T_NATIVE_DOUBLE);
+}
diff --git a/openEMS/tools/hdf5_file_writer.h b/openEMS/tools/hdf5_file_writer.h
new file mode 100644
index 0000000..7ad975c
--- /dev/null
+++ b/openEMS/tools/hdf5_file_writer.h
@@ -0,0 +1,68 @@
+/*
+*	Copyright (C) 2011 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef HDF5_FILE_WRITER_H
+#define HDF5_FILE_WRITER_H
+
+#include <string>
+#include <vector>
+#include <complex>
+#include <hdf5.h>
+
+class HDF5_File_Writer
+{
+public:
+	HDF5_File_Writer(std::string filename);
+	~HDF5_File_Writer();
+
+	bool WriteRectMesh(unsigned int const* numLines, double const* const* discLines, int MeshType=0, double scaling=1);
+	bool WriteRectMesh(unsigned int const* numLines, float const* const* discLines, int MeshType=0, float scaling=1);
+
+	bool WriteScalarField(std::string dataSetName, float const* const* const* field, size_t datasize[3]);
+	bool WriteScalarField(std::string dataSetName, double const* const* const* field, size_t datasize[3]);
+
+	bool WriteScalarField(std::string dataSetName, std::complex<float> const* const* const* field, size_t datasize[3]);
+	bool WriteScalarField(std::string dataSetName, std::complex<double> const* const* const* field, size_t datasize[3]);
+
+	bool WriteVectorField(std::string dataSetName, float const* const* const* const* field, size_t datasize[3]);
+	bool WriteVectorField(std::string dataSetName, double const* const* const* const* field, size_t datasize[3]);
+
+	bool WriteVectorField(std::string dataSetName, std::complex<float> const* const* const* const* field, size_t datasize[3]);
+	bool WriteVectorField(std::string dataSetName, std::complex<double> const* const* const* const* field, size_t datasize[3]);
+
+	bool WriteData(std::string dataSetName, float const* field_buf, size_t dim, size_t* datasize);
+	bool WriteData(std::string dataSetName, double const* field_buf, size_t dim, size_t* datasize);
+
+	bool WriteAtrribute(std::string locName, std::string attr_name, void const* value, hsize_t size, hid_t mem_type);
+	bool WriteAtrribute(std::string locName, std::string attr_name, float const* value, hsize_t size);
+	bool WriteAtrribute(std::string locName, std::string attr_name, double const* value, hsize_t size);
+	bool WriteAtrribute(std::string locName, std::string attr_name, std::vector<float> values);
+	bool WriteAtrribute(std::string locName, std::string attr_name, std::vector<double> values);
+	bool WriteAtrribute(std::string locName, std::string attr_name, float value);
+	bool WriteAtrribute(std::string locName, std::string attr_name, double value);
+
+	void SetCurrentGroup(std::string group, bool createGrp=true);
+
+protected:
+	std::string m_filename;
+	std::string m_Group;
+
+	hid_t OpenGroup(hid_t hdf5_file, std::string group);
+	bool WriteData(std::string dataSetName, hid_t mem_type, void const* field_buf, size_t dim, size_t* datasize);
+};
+
+#endif // HDF5_FILE_WRITER_H
diff --git a/openEMS/tools/sar_calculation.cpp b/openEMS/tools/sar_calculation.cpp
new file mode 100644
index 0000000..343aa07
--- /dev/null
+++ b/openEMS/tools/sar_calculation.cpp
@@ -0,0 +1,656 @@
+/*
+*	Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#include "sar_calculation.h"
+#include "cfloat"
+#include "array_ops.h"
+#include "global.h"
+
+using namespace std;
+
+SAR_Calculation::SAR_Calculation()
+{
+	m_Vx_Used = NULL;
+	m_Vx_Valid = NULL;
+	m_DebugLevel = 0;
+	SetAveragingMethod(SIMPLE, true);
+	Reset();
+}
+
+void SAR_Calculation::Reset()
+{
+	Delete3DArray(m_Vx_Used,m_numLines);
+	m_Vx_Used = NULL;
+	Delete3DArray(m_Vx_Valid,m_numLines);
+	m_Vx_Valid = NULL;
+
+	m_avg_mass = 0;
+	m_numLines[0]=m_numLines[1]=m_numLines[2]=0;
+	m_cellWidth[0]=m_cellWidth[1]=m_cellWidth[2]=NULL;
+
+	m_cell_volume = NULL;
+	m_cell_density = NULL;
+	m_cell_conductivity = NULL;
+	m_E_field = NULL;
+	m_J_field = NULL;
+
+	Delete3DArray(m_Vx_Used,m_numLines);
+	m_Vx_Used = NULL;
+	Delete3DArray(m_Vx_Valid,m_numLines);
+	m_Vx_Valid = NULL;
+}
+
+void SAR_Calculation::SetAveragingMethod(string method, bool silent)
+{
+	if (method.compare("IEEE_C95_3")==0)
+		return SetAveragingMethod(IEEE_C95_3, silent);
+	if (method.compare("IEEE_62704")==0)
+		return SetAveragingMethod(IEEE_62704, silent);
+	if (method.compare("Simple")==0)
+		return SetAveragingMethod(SIMPLE, silent);
+
+	cerr << __func__ << ": Error, """ << method << """ is an unknown averaging method..." << endl;
+	// unknown method, fallback to simple
+	SetAveragingMethod(SIMPLE, false);
+}
+
+void SAR_Calculation::SetAveragingMethod(SARAveragingMethod method, bool silent)
+{
+	if (method==IEEE_62704)
+	{
+		m_massTolerance = 0.0001;
+		m_maxMassIterations = 100;
+		m_maxBGRatio = 0.1;
+		m_markPartialAsUsed = false;
+		m_UnusedRelativeVolLimit = 1.05;
+		m_IgnoreFaceValid = false;
+		if (!silent)
+			cerr << __func__ << ": Setting averaging method to IEEE_62704" << endl;
+		return;
+	}
+	else if (method==IEEE_C95_3)
+	{
+		m_massTolerance = 0.05;
+		m_maxMassIterations = 100;
+		m_maxBGRatio = 1;
+		m_markPartialAsUsed = true;
+		m_UnusedRelativeVolLimit = 1;
+		m_IgnoreFaceValid = false;
+		if (!silent)
+			cerr << __func__ << ": Setting averaging method to IEEE_C95_3" << endl;
+		return;
+	}
+	else if (method==SIMPLE)
+	{
+		m_massTolerance = 0.05;
+		m_maxMassIterations = 100;
+		m_maxBGRatio = 1;
+		m_markPartialAsUsed = true;
+		m_UnusedRelativeVolLimit = 1;
+		m_IgnoreFaceValid = true;
+		if (!silent)
+			cerr << __func__ << ": Setting averaging method to SIMPLE" << endl;
+		return;
+	}
+
+	cerr << __func__ << ": Error, unknown averaging method..." << endl;
+	// unknown method, fallback to simple
+	SetAveragingMethod(SIMPLE, false);
+}
+
+void SAR_Calculation::SetNumLines(unsigned int numLines[3])
+{
+	Delete3DArray(m_Vx_Used,m_numLines);
+	m_Vx_Used = NULL;
+	Delete3DArray(m_Vx_Valid,m_numLines);
+	m_Vx_Valid = NULL;
+
+	for (int n=0;n<3;++n)
+		m_numLines[n]=numLines[n];
+}
+
+void SAR_Calculation::SetCellWidth(float* cellWidth[3])
+{
+	for (int n=0;n<3;++n)
+		m_cellWidth[n]=cellWidth[n];
+}
+
+float*** SAR_Calculation::CalcSAR(float*** SAR)
+{
+	if (CheckValid()==false)
+	{
+		cerr << "SAR_Calculation::CalcSAR: SAR calculation is invalid due to missing values... Abort..." << endl;
+		return NULL;
+	}
+	if (m_avg_mass<=0)
+		return CalcLocalSAR(SAR);
+	return CalcAveragedSAR(SAR);
+}
+
+bool SAR_Calculation::CheckValid()
+{
+	for (int n=0;n<3;++n)
+		if (m_cellWidth[n]==NULL)
+			return false;
+	if (m_E_field==NULL)
+		return false;
+	if ((m_J_field==NULL) && (m_cell_conductivity==NULL))
+		return false;
+	if (m_cell_density==NULL)
+		return false;
+	if (m_avg_mass<0)
+		return false;
+	return true;
+}
+
+double SAR_Calculation::CalcSARPower()
+{
+	if (CheckValid()==false)
+	{
+		cerr << "SAR_Calculation::CalcSARPower: SAR calculation is invalid due to missing values... Abort..." << endl;
+		return 0;
+	}
+	double power=0;
+	unsigned int pos[3];
+	for (pos[0]=0; pos[0]<m_numLines[0]; ++pos[0])
+	{
+		for (pos[1]=0; pos[1]<m_numLines[1]; ++pos[1])
+		{
+			for (pos[2]=0; pos[2]<m_numLines[2]; ++pos[2])
+			{
+				power += CalcLocalPowerDensity(pos)*CellVolume(pos);
+			}
+		}
+	}
+	return power;
+}
+
+double SAR_Calculation::CalcLocalPowerDensity(unsigned int pos[3])
+{
+	double l_pow=0;
+	if (m_cell_conductivity==NULL)
+	{
+		l_pow  = abs(m_E_field[0][pos[0]][pos[1]][pos[2]]) * abs(m_J_field[0][pos[0]][pos[1]][pos[2]]);
+		l_pow += abs(m_E_field[1][pos[0]][pos[1]][pos[2]]) * abs(m_J_field[1][pos[0]][pos[1]][pos[2]]);
+		l_pow += abs(m_E_field[2][pos[0]][pos[1]][pos[2]]) * abs(m_J_field[2][pos[0]][pos[1]][pos[2]]);
+	}
+	else
+	{
+		l_pow  = m_cell_conductivity[pos[0]][pos[1]][pos[2]]*abs(m_E_field[0][pos[0]][pos[1]][pos[2]]) * abs(m_E_field[0][pos[0]][pos[1]][pos[2]]);
+		l_pow += m_cell_conductivity[pos[0]][pos[1]][pos[2]]*abs(m_E_field[1][pos[0]][pos[1]][pos[2]]) * abs(m_E_field[1][pos[0]][pos[1]][pos[2]]);
+		l_pow += m_cell_conductivity[pos[0]][pos[1]][pos[2]]*abs(m_E_field[2][pos[0]][pos[1]][pos[2]]) * abs(m_E_field[2][pos[0]][pos[1]][pos[2]]);
+	}
+	l_pow*=0.5;
+	return l_pow;
+}
+
+float*** SAR_Calculation::CalcLocalSAR(float*** SAR)
+{
+	unsigned int pos[3];
+	m_Valid=0;
+	m_Used=0;
+	m_Unused=0;
+	m_AirVoxel=0;
+	for (pos[0]=0; pos[0]<m_numLines[0]; ++pos[0])
+	{
+		for (pos[1]=0; pos[1]<m_numLines[1]; ++pos[1])
+		{
+			for (pos[2]=0; pos[2]<m_numLines[2]; ++pos[2])
+			{
+				if (m_cell_density[pos[0]][pos[1]][pos[2]]>0)
+				{
+					++m_Valid;
+					SAR[pos[0]][pos[1]][pos[2]]	= CalcLocalPowerDensity(pos)/m_cell_density[pos[0]][pos[1]][pos[2]];
+				}
+				else
+				{
+					++m_AirVoxel;
+					SAR[pos[0]][pos[1]][pos[2]] = 0;
+				}
+			}
+		}
+	}
+	return SAR;
+}
+
+int SAR_Calculation::FindFittingCubicalMass(unsigned int pos[3], float box_size, unsigned int start[3], unsigned int stop[3],
+					float partial_start[3], float partial_stop[3], double &mass, double &volume, double &bg_ratio, int disabledFace, bool ignoreFaceValid)
+{
+	unsigned int mass_iterations = 0;
+	double old_mass=0;
+	double old_box_size=0;
+	bool face_valid;
+	bool mass_valid;
+	bool voxel_valid;
+
+	//iterate over cubical sizes to find fitting volume to mass
+	while (mass_iterations<m_maxMassIterations)
+	{
+		// calculate cubical mass
+		face_valid = GetCubicalMass(pos, box_size/2, start, stop, partial_start, partial_stop, mass, volume, bg_ratio, disabledFace);
+
+		// check if found mass is valid
+		mass_valid = abs(mass-m_avg_mass)<=m_massTolerance*m_avg_mass;
+		voxel_valid = mass_valid && (face_valid==true) && (bg_ratio<m_maxBGRatio);
+
+		if ((face_valid==false) && (mass<m_avg_mass*(1.0-m_massTolerance)) && (ignoreFaceValid==false))
+		{
+			// this is an invalid cube with a too small total mass --> increasing the box will not yield a valid cube
+			return 1;
+		}
+		else if ((face_valid==false  || bg_ratio>=m_maxBGRatio) && (mass_valid))
+		{
+			// this is an invalid cube with a valid total mass
+			if (ignoreFaceValid)
+				return 0;
+			return 2;
+		}
+		if (voxel_valid)
+		{
+			// valid cube found
+			return 0;
+		}
+
+		// if no valid or finally invalid cube is found, calculate an alternaive cube size
+		if (mass_iterations==0)
+		{
+			// on first interation, try a relative resize
+			old_box_size=box_size;
+			box_size*=pow(m_avg_mass/mass,1.0/3.0);
+		}
+		else
+		{
+			// on later interations, try a newton approach
+			float new_box_size = box_size - (mass-m_avg_mass)/(mass-old_mass)*(box_size-old_box_size);
+			old_box_size = box_size;
+			box_size = new_box_size;
+		}
+		old_mass=mass;
+
+		++mass_iterations;
+	}
+
+	// m_maxMassIterations iterations are exhausted...
+	return -1;
+}
+
+bool SAR_Calculation::GetCubicalMass(unsigned int pos[3], double box_size, unsigned int start[3], unsigned int stop[3],
+									 float partial_start[3], float partial_stop[3], double &mass, double &volume, double &bg_ratio, int disabledFace)
+{
+	if ((box_size<=0) || isnan(box_size) || isinf(box_size))
+	{
+		cerr << "SAR_Calculation::GetCubicalMass: critical error: invalid averaging box size!! EXIT" << endl;
+		exit(-1);
+	}
+	bool face_valid=true;
+	for (int n=0;n<3;++n)
+	{
+		// check start position
+		start[n]=pos[n];
+		partial_start[n]=1;
+		float dist=m_cellWidth[n][pos[n]]/2;
+		if (disabledFace==2*n)
+			dist=box_size;
+		else
+		{
+			while (dist<box_size)
+			{
+				// box is outside of domain
+				if (start[n]==0)
+				{
+					partial_start[n]=-1;
+					break;
+				}
+				--start[n];
+				dist+=m_cellWidth[n][start[n]];
+			}
+
+			if (partial_start[n]!=-1)
+			{   // box ends inside stop[n] voxel
+				partial_start[n]=1-(dist-box_size)/m_cellWidth[n][start[n]];
+			}
+			if ((partial_start[n]!=-1) && (pos[n]==start[n]))
+				partial_start[n]=2*box_size/m_cellWidth[n][start[n]];
+		}
+
+		// check stop position
+		stop[n]=pos[n];
+		partial_stop[n]=1;
+		dist=m_cellWidth[n][pos[n]]/2;
+		if (disabledFace==2*n+1)
+			dist=box_size;
+		else
+		{
+			while (dist<box_size)
+			{
+				// box is outside of domain
+				if (stop[n]==m_numLines[n]-1)
+				{
+					partial_stop[n]=-1;
+					break;
+				}
+				++stop[n];
+				dist+=m_cellWidth[n][stop[n]];
+			}
+
+			if (partial_stop[n]!=-1)
+			{   // box ends inside stop[n] voxel
+				partial_stop[n]=1-(dist-box_size)/m_cellWidth[n][stop[n]];
+			}
+			if ((partial_stop[n]!=-1) && (pos[n]==stop[n]))
+				partial_stop[n]=2*box_size/m_cellWidth[n][stop[n]];
+		}
+	}
+
+	for (int n=0;n<3;++n)
+	{
+		if (partial_start[n]==-1)
+			face_valid=false;
+		if (partial_stop[n]==-1)
+			face_valid=false;
+	}
+
+	mass = 0;
+	volume = 0;
+	double bg_volume=0;
+	double weight[3];
+	unsigned int f_pos[3];
+	bool face_intersect[6] = {false,false,false,false,false,false};
+	for (f_pos[0]=start[0];f_pos[0]<=stop[0];++f_pos[0])
+	{
+		weight[0]=1;
+		if (f_pos[0]==start[0])
+			weight[0]*=abs(partial_start[0]);
+		if (f_pos[0]==stop[0])
+			weight[0]*=abs(partial_stop[0]);
+
+		for (f_pos[1]=start[1];f_pos[1]<=stop[1];++f_pos[1])
+		{
+			weight[1]=1;
+			if (f_pos[1]==start[1])
+				weight[1]*=abs(partial_start[1]);
+			if (f_pos[1]==stop[1])
+				weight[1]*=abs(partial_stop[1]);
+
+			for (f_pos[2]=start[2];f_pos[2]<=stop[2];++f_pos[2])
+			{
+				weight[2]=1;
+				if (f_pos[2]==start[2])
+					weight[2]*=abs(partial_start[2]);
+				if (f_pos[2]==stop[2])
+					weight[2]*=abs(partial_stop[2]);
+
+				mass += CellMass(f_pos)*weight[0]*weight[1]*weight[2];
+				volume += CellVolume(f_pos)*weight[0]*weight[1]*weight[2];
+
+				if (m_cell_density[f_pos[0]][f_pos[1]][f_pos[2]]==0)
+					bg_volume += CellVolume(f_pos)*weight[0]*weight[1]*weight[2];
+				else
+				{
+					for (int n=0;n<3;++n)
+					{
+						if (start[n]==f_pos[n])
+							face_intersect[2*n]=true;
+						if (stop[n]==f_pos[n])
+							face_intersect[2*n+1]=true;
+					}
+				}
+			}
+		}
+	}
+	//check if all bounds have intersected a material boundary
+	for (int n=0;n<6;++n)
+		face_valid *= face_intersect[n];
+
+	bg_ratio = bg_volume/volume;
+
+	return face_valid;
+}
+
+float SAR_Calculation::CalcCubicalSAR(float*** SAR, unsigned int pos[3], unsigned int start[3], unsigned int stop[3], float partial_start[3], float partial_stop[3], bool assignUsed)
+{
+	double power_mass=0;
+	double mass=0;
+	double weight[3];
+	unsigned int f_pos[3];
+	for (f_pos[0]=start[0];f_pos[0]<=stop[0];++f_pos[0])
+	{
+		weight[0]=1;
+		if (f_pos[0]==start[0])
+			weight[0]*=abs(partial_start[0]);
+		if (f_pos[0]==stop[0])
+			weight[0]*=abs(partial_stop[0]);
+
+		for (f_pos[1]=start[1];f_pos[1]<=stop[1];++f_pos[1])
+		{
+			weight[1]=1;
+			if (f_pos[1]==start[1])
+				weight[1]*=abs(partial_start[1]);
+			if (f_pos[1]==stop[1])
+				weight[1]*=abs(partial_stop[1]);
+
+			for (f_pos[2]=start[2];f_pos[2]<=stop[2];++f_pos[2])
+			{
+				weight[2]=1;
+				if (f_pos[2]==start[2])
+					weight[2]*=abs(partial_start[2]);
+				if (f_pos[2]==stop[2])
+					weight[2]*=abs(partial_stop[2]);
+
+				if (m_cell_density[f_pos[0]][f_pos[1]][f_pos[2]]>=0)
+				{
+					mass += CellMass(f_pos)*weight[0]*weight[1]*weight[2];
+					power_mass += CalcLocalPowerDensity(f_pos) * CellVolume(f_pos)*weight[0]*weight[1]*weight[2];
+				}
+			}
+		}
+	}
+	float vx_SAR = power_mass/mass;
+	if (SAR==NULL)
+		return vx_SAR;
+
+	SAR[pos[0]][pos[1]][pos[2]]=vx_SAR;
+
+	if (assignUsed==false)
+		return vx_SAR;
+
+	// assign SAR to all used voxel
+	bool is_partial[3];
+	for (f_pos[0]=start[0];f_pos[0]<=stop[0];++f_pos[0])
+	{
+		if ( ((f_pos[0]==start[0]) && (partial_start[0]!=1)) || ((f_pos[0]==stop[0]) && (partial_stop[0]!=1)) )
+			is_partial[0]=true;
+		else
+			is_partial[0]=false;
+
+		for (f_pos[1]=start[1];f_pos[1]<=stop[1];++f_pos[1])
+		{
+			if ( ((f_pos[1]==start[1]) && (partial_start[1]!=1)) || ((f_pos[1]==stop[1]) && (partial_stop[1]!=1)) )
+				is_partial[1]=true;
+			else
+				is_partial[1]=false;
+
+			for (f_pos[2]=start[2];f_pos[2]<=stop[2];++f_pos[2])
+			{
+				if ( ((f_pos[2]==start[2]) && (partial_start[2]!=1)) || ((f_pos[2]==stop[2]) && (partial_stop[2]!=1)) )
+					is_partial[2]=true;
+				else
+					is_partial[2]=false;
+
+				if ( (!is_partial[0] && !is_partial[1] && !is_partial[2]) || m_markPartialAsUsed)
+				{
+					if (!m_Vx_Valid[f_pos[0]][f_pos[1]][f_pos[2]] && (m_cell_density[f_pos[0]][f_pos[1]][f_pos[2]]>0))
+					{
+						m_Vx_Used[f_pos[0]][f_pos[1]][f_pos[2]]=true;
+						SAR[f_pos[0]][f_pos[1]][f_pos[2]]=max(SAR[f_pos[0]][f_pos[1]][f_pos[2]], vx_SAR);
+					}
+				}
+			}
+		}
+	}
+	return vx_SAR;
+}
+
+
+float*** SAR_Calculation::CalcAveragedSAR(float*** SAR)
+{
+	unsigned int pos[3];
+	m_Vx_Used = Create3DArray<bool>(m_numLines);
+	m_Vx_Valid = Create3DArray<bool>(m_numLines);
+
+	double voxel_volume;
+	double total_mass;
+	unsigned int start[3];
+	unsigned int stop[3];
+	float partial_start[3];
+	float partial_stop[3];
+	double bg_ratio;
+	int EC=0;
+
+	// debug counter
+	unsigned int cnt_case1=0;
+	unsigned int cnt_case2=0;
+	unsigned int cnt_NoConvergence=0;
+
+	m_Valid=0;
+	m_Used=0;
+	m_Unused=0;
+	m_AirVoxel=0;
+
+	for (pos[0]=0; pos[0]<m_numLines[0]; ++pos[0])
+	{
+		for (pos[1]=0; pos[1]<m_numLines[1]; ++pos[1])
+		{
+			for (pos[2]=0; pos[2]<m_numLines[2]; ++pos[2])
+			{
+				if (m_cell_density[pos[0]][pos[1]][pos[2]]==0)
+				{
+					SAR[pos[0]][pos[1]][pos[2]] = 0;
+					++m_AirVoxel;
+					continue;
+				}
+
+				// guess an initial box size and find a fitting cube
+				EC = FindFittingCubicalMass(pos, pow(m_avg_mass/m_cell_density[pos[0]][pos[1]][pos[2]],1.0/3.0)/2, start, stop,
+											partial_start, partial_stop, total_mass, voxel_volume, bg_ratio, -1, m_IgnoreFaceValid);
+
+				if (EC==0)
+				{
+					m_Vx_Valid[pos[0]][pos[1]][pos[2]] = true;
+					m_Vx_Used[pos[0]][pos[1]][pos[2]] = true;
+					++m_Valid;
+					CalcCubicalSAR(SAR, pos, start, stop, partial_start, partial_stop, true);
+				}
+				else if (EC==1)
+					++cnt_case1;
+				else if (EC==2)
+					++cnt_case2;
+				else if (EC==-1)
+					++cnt_NoConvergence;
+				else
+					cerr << "other EC" << EC << endl;
+			}
+		}
+	}
+	if (cnt_NoConvergence>0)
+	{
+		cerr << "SAR_Calculation::CalcAveragedSAR: Warning, for some voxel a valid averaging cube could not be found (no convergence)... " << endl;
+	}
+	if (m_DebugLevel>0)
+	{
+		cerr << "Number of invalid cubes (case 1): " << cnt_case1 << endl;
+		cerr << "Number of invalid cubes (case 2): " << cnt_case2 << endl;
+		cerr << "Number of invalid cubes (failed to converge): " << cnt_NoConvergence << endl;
+	}
+
+	// count all used and unused etc. + special handling of unused voxels!!
+	m_Used=0;
+	m_Unused=0;
+	for (pos[0]=0;pos[0]<m_numLines[0];++pos[0])
+	{
+		for (pos[1]=0;pos[1]<m_numLines[1];++pos[1])
+		{
+			for (pos[2]=0;pos[2]<m_numLines[2];++pos[2])
+			{
+				if (!m_Vx_Valid[pos[0]][pos[1]][pos[2]] && m_Vx_Used[pos[0]][pos[1]][pos[2]])
+					++m_Used;
+				if ((m_cell_density[pos[0]][pos[1]][pos[2]]>0) && !m_Vx_Valid[pos[0]][pos[1]][pos[2]] && !m_Vx_Used[pos[0]][pos[1]][pos[2]])
+				{
+					++m_Unused;
+
+					SAR[pos[0]][pos[1]][pos[2]] = 0;
+					double unused_volumes[6];
+					float unused_SAR[6];
+
+					double min_Vol=FLT_MAX;
+
+					// special handling of unused voxels:
+					for (int n=0;n<6;++n)
+					{
+						EC = FindFittingCubicalMass(pos, pow(m_avg_mass/m_cell_density[pos[0]][pos[1]][pos[2]],1.0/3.0)/2, start, stop,
+													partial_start, partial_stop, total_mass, unused_volumes[n], bg_ratio, n, true);
+						if ((EC!=0) && (EC!=2))
+						{
+							// this should not happen
+							cerr << "SAR_Calculation::CalcAveragedSAR: Error handling unused voxels, can't find fitting cubical averaging volume' " << endl;
+							unused_SAR[n]=0;
+						}
+						else
+						{
+							unused_SAR[n]=CalcCubicalSAR(NULL, pos, start, stop, partial_start, partial_stop, false);
+							min_Vol = min(min_Vol,unused_volumes[n]);
+						}
+					}
+					for (int n=0;n<6;++n)
+					{
+						if (unused_volumes[n]<=m_UnusedRelativeVolLimit*min_Vol)
+							SAR[pos[0]][pos[1]][pos[2]] = max(SAR[pos[0]][pos[1]][pos[2]],unused_SAR[n]);
+					}
+				}
+			}
+		}
+	}
+
+	if (m_Valid+m_Used+m_Unused+m_AirVoxel!=m_numLines[0]*m_numLines[1]*m_numLines[2])
+	{
+		cerr << "SAR_Calculation::CalcAveragedSAR: critical error, mismatch in voxel status count... EXIT" << endl;
+		exit(1);
+	}
+
+	if (m_DebugLevel>0)
+		cerr << "SAR_Calculation::CalcAveragedSAR: Stats: Valid=" << m_Valid << " Used=" << m_Used << " Unused=" << m_Unused << " Air-Voxel=" << m_AirVoxel << endl;
+
+	return SAR;
+}
+
+double SAR_Calculation::CellVolume(unsigned int pos[3])
+{
+	if (m_cell_volume)
+		return m_cell_volume[pos[0]][pos[1]][pos[2]];
+
+	double volume=1;
+	for (int n=0;n<3;++n)
+		volume*=m_cellWidth[n][pos[n]];
+	return volume;
+}
+
+double SAR_Calculation::CellMass(unsigned int pos[3])
+{
+	return m_cell_density[pos[0]][pos[1]][pos[2]]*CellVolume(pos);
+}
+
diff --git a/openEMS/tools/sar_calculation.h b/openEMS/tools/sar_calculation.h
new file mode 100644
index 0000000..58d9ab8
--- /dev/null
+++ b/openEMS/tools/sar_calculation.h
@@ -0,0 +1,122 @@
+/*
+*	Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef SAR_CALCULATION_H
+#define SAR_CALCULATION_H
+
+#include <complex>
+
+class SAR_Calculation
+{
+public:
+	SAR_Calculation();
+
+	enum SARAveragingMethod { IEEE_C95_3, IEEE_62704, SIMPLE};
+
+	//! Reset and initialize all values (will keep all SAR settings)
+	void Reset();
+
+	//! Set the debug level
+	void SetDebugLevel(int level) {m_DebugLevel=level;}
+
+	//! Set the used averaging method
+	void SetAveragingMethod(SARAveragingMethod method, bool silent=false);
+
+	//! Set the used averaging method
+	void SetAveragingMethod(std::string method, bool silent=false);
+
+	//! Set number of lines in all direcitions. (mandatory information)
+	void SetNumLines(unsigned int numLines[3]);
+	//! Set cell width in all direcitions. (mandatory information for averaging)
+	void SetCellWidth(float* cellWidth[3]);
+
+	//! Set the averaging mash. (mandatory information for averaging)
+	void SetAveragingMass(float mass) {m_avg_mass=mass;}
+
+	//! Set the cell volumes (optional for speedup)
+	void SetCellVolumes(float*** cell_volume) {m_cell_volume=cell_volume;}
+
+	//! Set the cell densities (mandatory information)
+	void SetCellDensities(float*** cell_density) {m_cell_density=cell_density;}
+
+	//! Set the cell conductivities (mandatory if no current density field is given)
+	void SetCellCondictivity(float*** cell_conductivity) {m_cell_conductivity=cell_conductivity;}
+
+	//! Set the electric field (mandatory information)
+	void SetEField(std::complex<float>**** field) {m_E_field=field;}
+	//! Set the current density field (mandatory if no conductivity distribution is given)
+	void SetJField(std::complex<float>**** field) {m_J_field=field;}
+
+	//! Calculate the SAR, requires a preallocated 3D array
+	float*** CalcSAR(float*** SAR);
+
+	//! Calculate the total power dumped
+	double CalcSARPower();
+
+protected:
+	unsigned int m_numLines[3];
+	float* m_cellWidth[3];
+
+	float m_avg_mass;
+	float*** m_cell_volume;
+	float*** m_cell_density;
+	float*** m_cell_conductivity;
+	std::complex<float>**** m_E_field;
+	std::complex<float>**** m_J_field;
+
+	bool*** m_Vx_Used;
+	bool*** m_Vx_Valid;
+
+	unsigned int m_Valid;
+	unsigned int m_Used;
+	unsigned int m_Unused;
+	unsigned int m_AirVoxel;
+
+	int m_DebugLevel;
+
+	/*********** SAR calculation parameter and settings ***********/
+	float m_massTolerance;
+	unsigned int m_maxMassIterations;
+	float m_maxBGRatio;
+	bool m_markPartialAsUsed;
+	float m_UnusedRelativeVolLimit;
+	bool m_IgnoreFaceValid;
+
+	/*********** SAR calculations methods ********/
+	double CalcLocalPowerDensity(unsigned int pos[3]);
+
+	//! Calculate the local SAR
+	float*** CalcLocalSAR(float*** SAR);
+
+	/****** start SAR averaging and all necessary methods ********/
+	//! Calculate the averaged SAR
+	float*** CalcAveragedSAR(float*** SAR);
+
+	int FindFittingCubicalMass(unsigned int pos[3], float box_size, unsigned int start[3], unsigned int stop[3],
+						float partial_start[3], float partial_stop[3], double &mass, double &volume, double &bg_ratio, int disabledFace=-1, bool ignoreFaceValid=false);
+	bool GetCubicalMass(unsigned int pos[3], double box_size, unsigned int start[3], unsigned int stop[3],
+						float partial_start[3], float partial_stop[3], double &mass, double &volume, double &bg_ratio, int disabledFace=-1);
+
+	float CalcCubicalSAR(float*** SAR, unsigned int pos[3], unsigned int start[3], unsigned int stop[3], float partial_start[3], float partial_stop[3], bool assignUsed=false);
+	/****** end SAR averaging and all necessary methods ********/
+
+	bool CheckValid();
+	double CellVolume(unsigned int pos[3]);
+	double CellMass(unsigned int pos[3]);
+};
+
+#endif // SAR_CALCULATION_H
diff --git a/openEMS/tools/useful.cpp b/openEMS/tools/useful.cpp
new file mode 100644
index 0000000..50aacfb
--- /dev/null
+++ b/openEMS/tools/useful.cpp
@@ -0,0 +1,185 @@
+/*
+*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "useful.h"
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+#include <climits>
+#include <stdio.h>
+#include <stdlib.h>
+#include <iomanip>
+#include <sstream>
+#include <boost/algorithm/string.hpp>
+#include <iostream>
+
+unsigned int CalcNyquistNum(double fmax, double dT)
+{
+	if (fmax==0) return UINT_MAX;
+	if (dT==0) return 1;
+	double T0 = 1/fmax;
+	return floor(T0/2/dT);
+}
+
+double CalcNyquistFrequency(unsigned int nyquist, double dT)
+{
+	if (nyquist==0) return 0;
+	if (dT==0) return 0;
+	return floor(1/(double)nyquist/2/dT);
+}
+
+std::vector<unsigned int> AssignJobs2Threads(unsigned int jobs, unsigned int nrThreads, bool RemoveEmpty)
+{
+	std::vector<unsigned int> jpt; //jobs per thread
+
+	unsigned int ui_jpt = jobs/nrThreads;
+	for (unsigned int n=0; n<nrThreads; ++n)
+	{
+		jpt.push_back(ui_jpt);
+		jobs-=ui_jpt;
+	}
+
+	for (unsigned int n=0; n<nrThreads; ++n)
+	{
+		if (jobs>0)
+		{
+			++jpt.at(n);
+			--jobs;
+		}
+	}
+
+	if (jobs>0)
+		std::cerr << "AssignJobs2Threads: Error, " << jobs << " remain to be assigned, this should not have happend..." << std::endl;
+
+	if (RemoveEmpty)
+	{
+		while (jpt.back()==0)
+			jpt.pop_back();
+	}
+
+	return jpt;
+}
+
+std::vector<float> SplitString2Float(std::string str, std::string delimiter)
+{
+	std::vector<float> v_f;
+	std::vector<std::string> results;
+	boost::split(results, str, boost::is_any_of(delimiter));
+
+	for (size_t n=0;n<results.size();++n)
+	{
+		std::istringstream is(results.at(n));
+		float num;
+		if (is >> num)
+			v_f.push_back(num);
+	}
+	return v_f;
+}
+
+std::vector<double> SplitString2Double(std::string str, std::string delimiter)
+{
+	std::vector<double> v_f;
+	std::vector<std::string> results;
+	boost::split(results, str, boost::is_any_of(delimiter));
+
+	for (size_t n=0;n<results.size();++n)
+	{
+		std::istringstream is(results.at(n));
+		double num;
+		if (is >> num)
+			v_f.push_back(num);
+	}
+	return v_f;
+}
+
+bool CrossProd(const double *v1, const double *v2, double* out)
+{
+	int nP,nPP;
+	for (int n=0;n<3;++n)
+	{
+		nP = (n+1)%3;
+		nPP = (n+2)%3;
+		out[n] = v1[nP]*v2[nPP] - v1[nPP]*v2[nP];
+	}
+	return ((out[0]+out[1]+out[2])>0);
+}
+
+double ScalarProd(const double *v1, const double *v2)
+{
+	double out=0;
+	for (int n=0;n<3;++n)
+		out+=v1[n]*v2[n];
+	return out;
+}
+
+double Determinant(const double *mat)
+{
+	return mat[0]*mat[4]*mat[8]+mat[1]*mat[5]*mat[6]+mat[2]*mat[3]*mat[7]-mat[2]*mat[4]*mat[6]-mat[1]*mat[3]*mat[8]-mat[0]*mat[5]*mat[7];
+}
+
+double* Invert(const double* in, double* out)
+{
+	double det = Determinant(in);
+	out[0] = (in[4]*in[8]-in[5]*in[7])/det;
+	out[1] = (in[2]*in[7]-in[1]*in[8])/det;
+	out[2] = (in[1]*in[5]-in[2]*in[4])/det;
+	out[3] = (in[5]*in[6]-in[3]*in[8])/det;
+	out[4] = (in[0]*in[8]-in[2]*in[6])/det;
+	out[5] = (in[2]*in[3]-in[0]*in[5])/det;
+	out[6] = (in[3]*in[7]-in[4]*in[6])/det;
+	out[7] = (in[1]*in[6]-in[0]*in[7])/det;
+	out[8] = (in[0]*in[4]-in[1]*in[3])/det;
+	return out;
+}
+
+int LinePlaneIntersection(const double *p0, const double *p1, const double *p2, const double *l_start, const double *l_stop, double* is_point, double &dist)
+{
+	dist = 0;
+	double mat[9];
+	for (int n=0;n<3;++n)
+	{
+		is_point[n] = 0;
+		mat[3*n] = l_start[n]-l_stop[n];
+		mat[3*n+1] = p1[n]-p0[n];
+		mat[3*n+2] = p2[n]-p0[n];
+	}
+	double det = Determinant(mat);
+	if (fabs(det)<1e-50)
+		return -1;
+
+	double inv_mat[9];
+	Invert(mat, inv_mat);
+
+	double t=0,u=0,v=0;
+	for (int n=0;n<3;++n)
+	{
+		t+=inv_mat[n]*(l_start[n]-p0[n]);
+		u+=inv_mat[3+n]*(l_start[n]-p0[n]);
+		v+=inv_mat[6+n]*(l_start[n]-p0[n]);
+	}
+	dist = t;
+
+	for (int n=0;n<3;++n)
+		is_point[n] = l_start[n]*(1-dist) + l_stop[n]*dist;
+
+	if ((u<0) || (u>1) || (v<0) || (v>1))
+		return 1;
+	if ((t<0) || (t>1))
+		return 2;
+
+	return 0;
+}
diff --git a/openEMS/tools/useful.h b/openEMS/tools/useful.h
new file mode 100644
index 0000000..dec13ab
--- /dev/null
+++ b/openEMS/tools/useful.h
@@ -0,0 +1,44 @@
+/*
+*	Copyright (C) 2010 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef USEFUL_H
+#define USEFUL_H
+
+#include <vector>
+#include <string>
+
+//! Calc the nyquist number of timesteps for a given frequency and timestep
+unsigned int CalcNyquistNum(double fmax, double dT);
+
+//! Calc the highest frequency allowed for a given nyquist number of timesteps and timestep
+double CalcNyquistFrequency(unsigned int nyquist, double dT);
+
+//! Calculate an optimal job distribution to a given number of threads. Will return a vector with the jobs for each thread.
+std::vector<unsigned int> AssignJobs2Threads(unsigned int jobs, unsigned int nrThreads, bool RemoveEmpty=false);
+
+std::vector<float> SplitString2Float(std::string str, std::string delimiter=",");
+std::vector<double> SplitString2Double(std::string str, std::string delimiter=",");
+
+bool CrossProd(const double* v1, const double* v2, double* out);
+double ScalarProd(const double* v1, const double* v2);
+
+double Determinant(const double* mat);
+double* Invert(const double* in, double* out);
+
+int LinePlaneIntersection(const double *p0, const double* p1, const double* p2, const double* l_start, const double* l_stop, double* is_point, double &dist);
+
+#endif // USEFUL_H
diff --git a/openEMS/tools/vtk_file_writer.cpp b/openEMS/tools/vtk_file_writer.cpp
new file mode 100644
index 0000000..79c40f3
--- /dev/null
+++ b/openEMS/tools/vtk_file_writer.cpp
@@ -0,0 +1,340 @@
+/*
+*	Copyright (C) 2011,2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+using namespace std;
+
+#include "vtk_file_writer.h"
+
+#include <vtkRectilinearGrid.h>
+#include <vtkRectilinearGridWriter.h>
+#include <vtkXMLRectilinearGridWriter.h>
+#include <vtkStructuredGrid.h>
+#include <vtkStructuredGridWriter.h>
+#include <vtkXMLStructuredGridWriter.h>
+#include <vtkZLibDataCompressor.h>
+#include <vtkFloatArray.h>
+#include <vtkDoubleArray.h>
+#include <vtkFieldData.h>
+#include <vtkPointData.h>
+
+#include <sstream>
+#include <iomanip>
+
+
+VTK_File_Writer::VTK_File_Writer(string filename, int meshType)
+{
+	SetFilename(filename);
+	m_MeshType = meshType;
+	m_NativeDump = false;
+	m_Binary = true;
+	m_Compress = true;
+	m_AppendMode = false;
+	m_ActiveTS = false;
+
+	if (m_MeshType==0) //cartesian mesh
+		m_GridData = vtkRectilinearGrid::New();
+	else if (m_MeshType==1) //cylindrical mesh
+		m_GridData = vtkStructuredGrid::New();
+	else
+	{
+		cerr << "VTK_File_Writer::VTK_File_Writer: Error, unknown mesh type: " << m_MeshType << endl;
+		m_GridData=NULL;
+	}
+}
+
+VTK_File_Writer::~VTK_File_Writer()
+{
+	if (m_GridData)
+		m_GridData->Delete();
+	m_GridData = NULL;
+}
+
+void VTK_File_Writer::SetMeshLines(double const* const* lines, unsigned int const* count, double scaling)
+{
+	if (m_MeshType==0) //cartesian mesh
+	{
+		vtkRectilinearGrid* RectGrid = dynamic_cast<vtkRectilinearGrid*>(m_GridData);
+		if (RectGrid==NULL)
+		{
+			cerr << "VTK_File_Writer::SetMeshLines: Error, grid invalid, this should not have happend! " << endl;
+			exit(1);
+		}
+		RectGrid->SetDimensions(count[0],count[1],count[2]);
+		vtkDoubleArray *Coords[3];
+		for (int n=0;n<3;++n)
+		{
+			m_MeshLines[n].clear();
+			m_MeshLines[n].reserve(count[n]);
+			Coords[n] = vtkDoubleArray::New();
+			for (unsigned int i=0; i<count[n]; i++)
+			{
+				Coords[n]->InsertNextValue(lines[n][i]*scaling);
+				m_MeshLines[n].push_back(lines[n][i]*scaling);
+			}
+		}
+		RectGrid->SetXCoordinates(Coords[0]);
+		RectGrid->SetYCoordinates(Coords[1]);
+		RectGrid->SetZCoordinates(Coords[2]);
+		for (int n=0;n<3;++n)
+			Coords[n]->Delete();
+	}
+	else if (m_MeshType==1) //cylindrical mesh
+	{
+		vtkStructuredGrid* StructGrid = dynamic_cast<vtkStructuredGrid*>(m_GridData);
+		if (StructGrid==NULL)
+		{
+			cerr << "VTK_File_Writer::SetMeshLines: Error, grid invalid, this should not have happend! " << endl;
+			exit(1);
+		}
+
+		for (int n=0;n<3;++n)
+		{
+			m_MeshLines[n].clear();
+			m_MeshLines[n].reserve(count[n]);
+			double scale=1;
+			if (n!=1)
+				scale*=scaling;
+			for (unsigned int i=0; i<count[n]; i++)
+				m_MeshLines[n].push_back(lines[n][i]*scale);
+		}
+
+		StructGrid->SetDimensions(count[0],count[1],count[2]);
+		vtkPoints *points = vtkPoints::New();
+		points->SetNumberOfPoints(count[0]*count[1]*count[2]);
+		double r[3];
+		int id=0;
+		for (unsigned int k=0; k<count[2]; ++k)
+			for (unsigned int j=0; j<count[1]; ++j)
+				for (unsigned int i=0; i<count[0]; ++i)
+				{
+					r[0] = lines[0][i] * cos(lines[1][j]) * scaling;
+					r[1] = lines[0][i] * sin(lines[1][j]) * scaling;
+					r[2] = lines[2][k] * scaling;
+					points->SetPoint(id++,r);
+				}
+		StructGrid->SetPoints(points);
+		points->Delete();
+	}
+	else
+	{
+		cerr << "VTK_File_Writer::SetMeshLines: Error, unknown mesh type: " << m_MeshType << endl;
+	}
+}
+
+void VTK_File_Writer::AddScalarField(string fieldname, double const* const* const* field)
+{
+	vtkDoubleArray* array = vtkDoubleArray::New();
+	array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
+	array->SetName(fieldname.c_str());
+	int id=0;
+	for (unsigned int k=0;k<m_MeshLines[2].size();++k)
+	{
+		for (unsigned int j=0;j<m_MeshLines[1].size();++j)
+		{
+			for (unsigned int i=0;i<m_MeshLines[0].size();++i)
+			{
+				array->SetTuple1(id++,field[i][j][k]);
+			}
+		}
+	}
+	m_GridData->GetPointData()->AddArray(array);
+	array->Delete();
+}
+
+void VTK_File_Writer::AddScalarField(string fieldname, float const* const* const* field)
+{
+	vtkFloatArray* array = vtkFloatArray::New();
+	array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
+	array->SetName(fieldname.c_str());
+	int id=0;
+	for (unsigned int k=0;k<m_MeshLines[2].size();++k)
+	{
+		for (unsigned int j=0;j<m_MeshLines[1].size();++j)
+		{
+			for (unsigned int i=0;i<m_MeshLines[0].size();++i)
+			{
+				array->SetTuple1(id++,field[i][j][k]);
+			}
+		}
+	}
+	m_GridData->GetPointData()->AddArray(array);
+	array->Delete();
+}
+
+void VTK_File_Writer::AddVectorField(string fieldname, double const* const* const* const* field)
+{
+	vtkDoubleArray* array = vtkDoubleArray::New();
+	array->SetNumberOfComponents(3);
+	array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
+	array->SetName(fieldname.c_str());
+	int id=0;
+	double out[3];
+	for (unsigned int k=0;k<m_MeshLines[2].size();++k)
+	{
+		for (unsigned int j=0;j<m_MeshLines[1].size();++j)
+		{
+			double cos_a = cos(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
+			double sin_a = sin(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
+			for (unsigned int i=0;i<m_MeshLines[0].size();++i)
+			{
+				if ((m_MeshType==0) || (m_NativeDump))
+					array->SetTuple3(id++,field[0][i][j][k],field[1][i][j][k],field[2][i][j][k]);
+				else
+				{
+					out[0] = field[0][i][j][k] * cos_a - field[1][i][j][k] * sin_a;
+					out[1] = field[0][i][j][k] * sin_a + field[1][i][j][k] * cos_a;
+					out[2] = field[2][i][j][k];
+					array->SetTuple3(id++,out[0],out[1],out[2]);
+				}
+			}
+		}
+	}
+	m_GridData->GetPointData()->AddArray(array);
+	array->Delete();
+}
+
+void VTK_File_Writer::AddVectorField(string fieldname, float const* const* const* const* field)
+{
+	vtkFloatArray* array = vtkFloatArray::New();
+	array->SetNumberOfComponents(3);
+	array->SetNumberOfTuples(m_MeshLines[0].size()*m_MeshLines[1].size()*m_MeshLines[2].size());
+	array->SetName(fieldname.c_str());
+	int id=0;
+	float out[3];
+	for (unsigned int k=0;k<m_MeshLines[2].size();++k)
+	{
+		for (unsigned int j=0;j<m_MeshLines[1].size();++j)
+		{
+			float cos_a = cos(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
+			float sin_a = sin(m_MeshLines[1].at(j)); //needed only for m_MeshType==1 (cylindrical mesh)
+			for (unsigned int i=0;i<m_MeshLines[0].size();++i)
+			{
+				if ((m_MeshType==0) || (m_NativeDump))
+					array->SetTuple3(id++,field[0][i][j][k],field[1][i][j][k],field[2][i][j][k]);
+				else
+				{
+					out[0] = field[0][i][j][k] * cos_a - field[1][i][j][k] * sin_a;
+					out[1] = field[0][i][j][k] * sin_a + field[1][i][j][k] * cos_a;
+					out[2] = field[2][i][j][k];
+					array->SetTuple3(id++,out[0],out[1],out[2]);
+				}
+			}
+		}
+	}
+	m_GridData->GetPointData()->AddArray(array);
+	array->Delete();
+}
+
+
+int VTK_File_Writer::GetNumberOfFields() const
+{
+	return m_GridData->GetPointData()->GetNumberOfArrays();
+}
+
+void VTK_File_Writer::ClearAllFields()
+{
+	while (m_GridData->GetPointData()->GetNumberOfArrays()>0)
+	{
+		const char* name = m_GridData->GetPointData()->GetArrayName(0);
+		m_GridData->GetPointData()->RemoveArray(name);
+	}
+}
+
+bool VTK_File_Writer::Write()
+{
+	return WriteXML();
+}
+
+string VTK_File_Writer::GetTimestepFilename(int pad_length) const
+{
+	if (m_ActiveTS==false)
+		return m_filename;
+
+	stringstream ss;
+	ss << m_filename << "_" << std::setw( pad_length ) << std::setfill( '0' ) << m_timestep;
+
+	return ss.str();
+}
+
+
+bool VTK_File_Writer::WriteASCII()
+{
+	vtkDataWriter* writer = NULL;
+	if (m_MeshType==0) //cartesian mesh
+		writer = vtkRectilinearGridWriter::New();
+	else if (m_MeshType==1) //cylindrical mesh
+		writer = vtkStructuredGridWriter::New();
+	else
+	{
+		cerr << "VTK_File_Writer::WriteASCII: Error, unknown mesh type: " << m_MeshType << endl;
+		return false;
+	}
+
+	writer->SetHeader(m_header.c_str());
+#if VTK_MAJOR_VERSION>=6
+	writer->SetInputData(m_GridData);
+#else
+	writer->SetInput(m_GridData);
+#endif
+
+	string filename = GetTimestepFilename() + ".vtk";
+	writer->SetFileName(filename.c_str());
+	if (m_Binary)
+		writer->SetFileTypeToBinary();
+	else
+		writer->SetFileTypeToASCII();
+
+	writer->Write();
+	writer->Delete();
+	return true;
+}
+
+bool VTK_File_Writer::WriteXML()
+{
+	vtkXMLStructuredDataWriter* writer = NULL;
+	if (m_MeshType==0) //cartesian mesh
+		writer = vtkXMLRectilinearGridWriter::New();
+	else if (m_MeshType==1) //cylindrical mesh
+		writer = vtkXMLStructuredGridWriter::New();
+	else
+	{
+		cerr << "VTK_File_Writer::WriteXML: Error, unknown mesh type: " << m_MeshType << endl;
+		return false;
+	}
+
+#if VTK_MAJOR_VERSION>=6
+	writer->SetInputData(m_GridData);
+#else
+	writer->SetInput(m_GridData);
+#endif
+
+	string filename = GetTimestepFilename() + "." + writer->GetDefaultFileExtension();
+	writer->SetFileName(filename.c_str());
+	if (m_Compress)
+		writer->SetCompressor(vtkZLibDataCompressor::New());
+	else
+		writer->SetCompressor(NULL);
+
+	if (m_Binary)
+		writer->SetDataModeToBinary();
+	else
+		writer->SetDataModeToAscii();
+
+	writer->Write();
+	writer->Delete();
+	return true;
+}
diff --git a/openEMS/tools/vtk_file_writer.h b/openEMS/tools/vtk_file_writer.h
new file mode 100644
index 0000000..428d859
--- /dev/null
+++ b/openEMS/tools/vtk_file_writer.h
@@ -0,0 +1,94 @@
+/*
+*	Copyright (C) 2011,2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef VTK_FILE_WRITER_H
+#define VTK_FILE_WRITER_H
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include <string>
+#include <vector>
+#include <complex>
+
+class vtkDataSet;
+
+class VTK_File_Writer
+{
+public:
+	VTK_File_Writer(std::string filename, int meshType=0);
+	virtual ~VTK_File_Writer();
+
+	//! Set the filename
+	virtual void SetFilename(std::string filename) {m_filename=filename;}
+	//! Set the header information. May not be supported by all file types or setting.
+	virtual void SetHeader(std::string header) {m_header=header;}
+
+	//! Tell write to append data. May fail if filename has changed or filetype doesn't support this.
+	virtual void SetAppendMode(bool val) {m_AppendMode=val;}
+	//! Set binary flag (if the file type supports it)
+	virtual void SetBinary(bool val) {m_Binary=val;}
+	//! Set compression flag (if the file type supports it)
+	virtual void SetCompress(bool val) {m_Compress=val;}
+
+	void SetNativeDump(bool val) {m_NativeDump=val;}
+
+	virtual void SetMeshLines(double const* const* lines, unsigned int const* count, double scaling=1);
+
+	virtual void AddScalarField(std::string fieldname, double const* const* const* field);
+	virtual void AddScalarField(std::string fieldname, float const* const* const* field);
+	virtual void AddVectorField(std::string fieldname, double const* const* const* const* field);
+	virtual void AddVectorField(std::string fieldname, float const* const* const* const* field);
+
+	virtual int GetNumberOfFields() const;
+	virtual void ClearAllFields();
+
+	//! Get if timestep file series is active. \sa SetTimestepActive
+	virtual bool GetTimestepActive() {return m_ActiveTS;}
+	//! Set the timestep file series flag. \sa GetTimestepActive \sa SetTimestep
+	virtual void SetTimestepActive(bool val) {m_ActiveTS = val;}
+	//! Set the current timestep, this will set the timestep flag to true. \sa SetTimestepActive
+	virtual void SetTimestep(unsigned int ts) {m_timestep=ts;SetTimestepActive(true);}
+
+	virtual bool Write();
+
+	virtual bool WriteASCII();
+	virtual bool WriteXML();
+
+protected:
+	std::string m_filename;
+	std::string m_header;
+
+	//timestep properties
+	bool m_ActiveTS;
+	unsigned int m_timestep;
+
+	vtkDataSet* m_GridData;
+
+	//mesh information
+	int m_MeshType;
+	std::vector<double> m_MeshLines[3];
+	bool m_NativeDump;
+
+	bool m_AppendMode;
+	bool m_Binary;
+	bool m_Compress;
+
+	virtual std::string GetTimestepFilename(int pad_length=10) const;
+};
+
+#endif // VTK_FILE_Writer_H