commit 6fa5b4140a6d8d06853f9ac70366ff1ccef11f69
Author: Thorsten Liebig <thorsten.liebig@gmx.de>
Date:   Sun Feb 28 22:48:03 2010 +0100

    Initial commit

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..c82e591
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,9 @@
+release/
+debug/
+doc/
+obj/
+Makefile*
+*~
+*.so*
+*.pro.user
+openEMS
diff --git a/FDTD/cartoperator.cpp b/FDTD/cartoperator.cpp
new file mode 100644
index 0000000..e361323
--- /dev/null
+++ b/FDTD/cartoperator.cpp
@@ -0,0 +1,448 @@
+#include "cartoperator.h"
+
+CartOperator::CartOperator()
+{
+    Init();
+}
+
+CartOperator::~CartOperator()
+{
+	Reset();
+}
+
+void CartOperator::Init()
+{
+    CSX = NULL;
+    MainOp=NULL;
+    DualOp=NULL;
+    for (int n=0;n<3;++n)
+    {
+		discLines[n]=NULL;
+        EC_C[n]=NULL;
+        EC_G[n]=NULL;
+        EC_L[n]=NULL;
+        EC_R[n]=NULL;
+		vv[n]=NULL;
+		vi[n]=NULL;
+		iv[n]=NULL;
+		ii[n]=NULL;
+	}
+}
+
+void CartOperator::Reset()
+{
+	for (int n=0;n<3;++n)
+	{
+		delete[] EC_C[n];
+		delete[] EC_G[n];
+		delete[] EC_L[n];
+		delete[] EC_R[n];
+		delete[] vv[n];
+		delete[] vi[n];
+		delete[] iv[n];
+		delete[] ii[n];
+	}
+	delete MainOp;
+	delete DualOp;
+	Init();
+}
+
+
+void CartOperator::SetGeometryCSX(ContinuousStructure* geo)
+{
+    if (geo==NULL) return;
+
+    Reset();
+    CSX = geo;
+
+    CSRectGrid* grid=CSX->GetGrid();
+    for (int n=0;n<3;++n)
+    {
+        discLines[n] = grid->GetLines(n,discLines[n],numLines[n],true);
+		if (numLines[n]<3) {cerr << "CartOperator::SetGeometryCSX: you need at least 3 disc-lines in every direction (3D!)!!!" << endl; Reset(); return;}
+    }
+    MainOp = new AdrOp(numLines[0],numLines[1],numLines[2]);
+    MainOp->SetGrid(discLines[0],discLines[1],discLines[2]);
+	if (grid->GetDeltaUnit()<=0)  {cerr << "CartOperator::SetGeometryCSX: grid delta unit must not be <=0 !!!" << endl; Reset(); return;}
+	else gridDelta=grid->GetDeltaUnit();
+	MainOp->SetGridDelta(1);
+    MainOp->AddCellAdrOp();
+}
+
+int CartOperator::CalcECOperator()
+{
+    if (Calc_EC()==0)
+        return -1;
+
+	CalcTimestep();
+
+	for (int n=0;n<3;++n)
+	{
+		delete[] vv[n];
+		vv[n] = new FDTD_FLOAT[MainOp->GetSize()];
+		delete[] vi[n];
+		vi[n] = new FDTD_FLOAT[MainOp->GetSize()];
+		delete[] iv[n];
+		iv[n] = new FDTD_FLOAT[MainOp->GetSize()];
+		delete[] ii[n];
+		ii[n] = new FDTD_FLOAT[MainOp->GetSize()];
+
+		for (unsigned int i=0;i<MainOp->GetSize();++i)
+		{
+			vv[n][i] = (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
+			vi[n][i] = (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
+
+			ii[n][i] = (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
+			iv[n][i] = (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
+//			cerr << iv[n][i] << endl;
+		}
+	}
+
+	//cleanup
+	for (int n=0;n<3;++n)
+	{
+		delete[] EC_C[n];EC_C[n]=NULL;
+		delete[] EC_G[n];EC_G[n]=NULL;
+		delete[] EC_L[n];EC_L[n]=NULL;
+		delete[] EC_R[n];EC_R[n]=NULL;
+	}
+
+	//Always apply PEC to all boundary's
+	bool PEC[6]={0,0,0,0,0,0};
+	ApplyElectricBC(PEC);
+
+    return 0;
+}
+
+void CartOperator::ApplyElectricBC(bool* dirs)
+{
+	if (dirs==NULL) return;
+	unsigned int pos[3];
+	unsigned int ipos;
+	for (int n=0;n<3;++n)
+	{
+		int nP = (n+1)%3;
+		int nPP = (n+2)%3;
+		for (pos[nP]=0;pos[nP]<numLines[nP];++pos[nP])
+		{
+			for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP])
+			{
+				pos[n]=0;
+				ipos=MainOp->SetPos(pos[0],pos[1],pos[2]);
+				vv[n][ipos] *= (FDTD_FLOAT)!dirs[2*n];
+				vi[n][ipos] *= (FDTD_FLOAT)!dirs[2*n];
+
+				pos[n]=numLines[n]-1;
+				ipos=MainOp->SetPos(pos[0],pos[1],pos[2]);
+				vv[n][ipos] *= (FDTD_FLOAT)!dirs[2*n+1];
+				vi[n][ipos] *= (FDTD_FLOAT)!dirs[2*n+1];
+			}
+		}
+	}
+}
+
+void CartOperator::ApplyMagneticBC(bool* dirs)
+{
+	if (dirs==NULL) return;
+	unsigned int pos[3];
+	unsigned int ipos;
+	for (int n=0;n<3;++n)
+	{
+		int nP = (n+1)%3;
+		int nPP = (n+2)%3;
+		for (pos[nP]=0;pos[nP]<numLines[nP];++pos[nP])
+		{
+			for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP])
+			{
+				pos[n]=0;
+				ipos=MainOp->SetPos(pos[0],pos[1],pos[2]);
+				ii[n][ipos] *= (FDTD_FLOAT)!dirs[2*n];
+				iv[n][ipos] *= (FDTD_FLOAT)!dirs[2*n];
+
+				pos[n]=numLines[n]-2;
+				ipos=MainOp->SetPos(pos[0],pos[1],pos[2]);
+				ii[n][ipos] *= (FDTD_FLOAT)!dirs[2*n+1];
+				iv[n][ipos] *= (FDTD_FLOAT)!dirs[2*n+1];
+			}
+		}
+	}
+}
+
+
+bool CartOperator::Calc_ECPos(int n, unsigned int* pos, double* inEC)
+{
+	double coord[3];
+	double shiftCoord[3];
+	int nP = (n+1)%3;
+	int nPP = (n+2)%3;
+	coord[0] = discLines[0][pos[0]];
+	coord[1] = discLines[1][pos[1]];
+	coord[2] = discLines[2][pos[2]];
+	unsigned int ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
+	double delta=MainOp->GetIndexDelta(n,pos[n]);
+	double deltaP=MainOp->GetIndexDelta(nP,pos[nP]);
+	double deltaPP=MainOp->GetIndexDelta(nPP,pos[nPP]);
+	double delta_M=MainOp->GetIndexDelta(n,pos[n]-1);
+	double deltaP_M=MainOp->GetIndexDelta(nP,pos[nP]-1);
+	double deltaPP_M=MainOp->GetIndexDelta(nPP,pos[nPP]-1);
+
+	//******************************* epsilon,kappa averaging *****************************//
+	//shift up-right
+	shiftCoord[n] = coord[n]+delta*0.5;
+	shiftCoord[nP] = coord[nP]+deltaP*0.25;
+	shiftCoord[nPP] = coord[nPP]+deltaPP*0.25;
+	CSProperties* prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
+	if (prop)
+	{
+		CSPropMaterial* mat = prop->ToMaterial();
+		inEC[0] = mat->GetEpsilon(n)*fabs(deltaP*deltaPP);
+		inEC[1] = mat->GetKappa(n)*fabs(deltaP*deltaPP);
+	}
+	else
+	{
+		inEC[0] = 1*fabs(deltaP*deltaPP);
+		inEC[1] = 0;
+	}
+	//shift up-left
+	shiftCoord[n] = coord[n]+delta*0.5;
+	shiftCoord[nP] = coord[nP]-deltaP_M*0.25;
+	shiftCoord[nPP] = coord[nPP]+deltaPP*0.25;
+	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
+	if (prop)
+	{
+		CSPropMaterial* mat = prop->ToMaterial();
+		inEC[0] += mat->GetEpsilon(n)*fabs(deltaP*deltaPP);
+		inEC[1] += mat->GetKappa(n)*fabs(deltaP*deltaPP);
+	}
+	else
+	{
+		inEC[0] += 1*fabs(deltaP*deltaPP);
+		inEC[1] += 0;
+	}
+
+	//shift down-right
+	shiftCoord[n] = coord[n]+delta*0.5;
+	shiftCoord[nP] = coord[nP]+deltaP*0.25;
+	shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25;
+	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
+	if (prop)
+	{
+		CSPropMaterial* mat = prop->ToMaterial();
+		inEC[0] += mat->GetEpsilon(n)*fabs(deltaP*deltaPP);
+		inEC[1] += mat->GetKappa(n)*fabs(deltaP*deltaPP);
+	}
+	else
+	{
+		inEC[0] += 1*fabs(deltaP*deltaPP);
+		inEC[1] += 0;
+	}
+
+	//shift down-left
+	shiftCoord[n] = coord[n]+delta*0.5;
+	shiftCoord[nP] = coord[nP]-deltaP_M*0.25;
+	shiftCoord[nPP] = coord[nPP]-deltaPP_M*0.25;
+	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
+	if (prop)
+	{
+		CSPropMaterial* mat = prop->ToMaterial();
+		inEC[0] += mat->GetEpsilon(n)*fabs(deltaP*deltaPP);
+		inEC[1] += mat->GetKappa(n)*fabs(deltaP*deltaPP);
+	}
+	else
+	{
+		inEC[0] += 1*fabs(deltaP*deltaPP);
+		inEC[1] += 0;
+	}
+
+	inEC[0]*=gridDelta/fabs(delta)/4*__EPS0__;
+	inEC[1]*=gridDelta/fabs(delta)/4;
+
+	//******************************* mu,sigma averaging *****************************//
+	//shift down
+	shiftCoord[n] = coord[n]-delta_M*0.25;
+	shiftCoord[nP] = coord[nP]+deltaP*0.5;
+	shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
+	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
+	if (prop)
+	{
+		CSPropMaterial* mat = prop->ToMaterial();
+		inEC[2] = fabs(delta_M) / mat->GetMue(n);
+		if (mat->GetSigma(n))
+			inEC[3] = fabs(delta_M) / mat->GetSigma(n);
+		else
+			inEC[3] = 0;
+	}
+	else
+	{
+		inEC[2] = fabs(delta_M);
+		inEC[3] = 0;
+	}
+	//shift up
+	shiftCoord[n] = coord[n]+delta*0.25;
+	shiftCoord[nP] = coord[nP]+deltaP*0.5;
+	shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
+	prop = CSX->GetPropertyByCoordPriority(shiftCoord,CSProperties::MATERIAL);
+	if (prop)
+	{
+		CSPropMaterial* mat = prop->ToMaterial();
+		inEC[2] += mat->GetMue(n)*fabs(delta);
+		if (mat->GetSigma(n))
+			inEC[3] += fabs(delta)/mat->GetSigma(n);
+		else
+			inEC[3] = 0;
+	}
+	else
+	{
+		inEC[2] += 1*fabs(delta);
+		inEC[3] = 0;
+	}
+
+	inEC[2] = gridDelta * fabs(deltaP*deltaPP) * 2 * __MUE0__ / inEC[2];
+	if (inEC[3]) inEC[3]=gridDelta*fabs(deltaP*deltaPP) * 2 / inEC[3];
+
+	return true;
+}
+
+bool CartOperator::Calc_EffMatPos(int n, unsigned int* pos, double* inMat)
+{
+	int nP = (n+1)%3;
+	int nPP = (n+2)%3;
+
+	unsigned int ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
+	double delta=MainOp->GetIndexDelta(n,pos[n]);
+	double deltaP=MainOp->GetIndexDelta(nP,pos[nP]);
+	double deltaPP=MainOp->GetIndexDelta(nPP,pos[nPP]);
+	
+	double delta_M=MainOp->GetIndexDelta(n,pos[n]-1);
+	double deltaP_M=MainOp->GetIndexDelta(nP,pos[nP]-1);
+	double deltaPP_M=MainOp->GetIndexDelta(nPP,pos[nPP]-1);
+
+	this->Calc_ECPos(n,pos,inMat);
+
+	inMat[0] *= (delta*delta)/MainOp->GetNodeVolume(ipos)/gridDelta;
+	inMat[1] *= (delta*delta)/MainOp->GetNodeVolume(ipos)/gridDelta;
+
+	inMat[2] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta;
+	inMat[3] *= 0.5*(fabs(delta_M) + fabs(delta)) / fabs(deltaP*deltaPP) / gridDelta;
+
+}
+
+
+bool CartOperator::Calc_EC()
+{
+	if (CSX==NULL) {cerr << "CartOperator::Calc_EC: CSX not given or invalid!!!" << endl; return false;}
+
+	unsigned int ipos;
+	unsigned int pos[3];
+	double inEC[4];
+	for (int n=0;n<3;++n)
+    {
+		//init x-cell-array
+		delete[] EC_C[n];
+		delete[] EC_G[n];
+		delete[] EC_L[n];
+		delete[] EC_R[n];
+		EC_C[n] = new double[MainOp->GetSize()];
+		EC_G[n] = new double[MainOp->GetSize()];
+		EC_L[n] = new double[MainOp->GetSize()];
+		EC_R[n] = new double[MainOp->GetSize()];
+		for (unsigned int i=0;i<MainOp->GetSize();i++) //init all
+		{
+			EC_C[n][i]=0;
+			EC_G[n][i]=0;
+			EC_L[n][i]=0;
+			EC_R[n][i]=0;
+		}
+		for (pos[2]=0;pos[2]<numLines[2];++pos[2])
+        {
+			for (pos[1]=0;pos[1]<numLines[1];++pos[1])
+            {
+				for (pos[0]=0;pos[0]<numLines[0];++pos[0])
+                {
+					Calc_ECPos(n,pos,inEC);
+					ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
+					EC_C[n][ipos]=inEC[0];
+					EC_G[n][ipos]=inEC[1];
+					EC_L[n][ipos]=inEC[2];
+					EC_R[n][ipos]=inEC[3];
+				}
+            }
+        }
+    }
+
+    return true;
+}
+
+double CartOperator::CalcTimestep()
+{
+	dT=1e200;
+	double newT;
+	unsigned int pos[3];
+	unsigned int ipos;
+	unsigned int ipos_PM;
+	unsigned int ipos_PPM;
+	MainOp->SetReflection2Cell();
+	for (int n=0;n<3;++n)
+	{
+		int nP = (n+1)%3;
+		int nPP = (n+2)%3;
+
+		for (pos[2]=0;pos[2]<numLines[2];++pos[2])
+		{
+			for (pos[1]=0;pos[1]<numLines[1];++pos[1])
+			{
+				for (pos[0]=0;pos[0]<numLines[0];++pos[0])
+				{
+					ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
+					ipos_PM = MainOp->Shift(nP,-1);
+					MainOp->ResetShift();
+					ipos_PPM= MainOp->Shift(nPP,-1);
+					MainOp->ResetShift();
+					newT = 2/sqrt( ( 4/EC_L[nP][ipos] + 4/EC_L[nP][ipos_PPM] + 4/EC_L[nPP][ipos] + 4/EC_L[nPP][ipos_PM]) / EC_C[n][ipos] );
+					if (newT<dT) dT=newT;
+				}
+			}
+		}
+	}
+//	cerr << "Timestep: " << dT << endl;
+	return 0;
+}
+
+unsigned int CartOperator::GetVoltageExcitation(unsigned int* &index, FDTD_FLOAT** &excit_amp, FDTD_FLOAT* &excit_delay)
+{
+	vector<unsigned int> vIndex;
+	vector<FDTD_FLOAT> vExcit[3];
+	vector<FDTD_FLOAT> vDelay;
+	unsigned int ipos;
+	unsigned int pos[3];
+	double coord[3];
+
+	for (pos[2]=0;pos[2]<numLines[2];++pos[2])
+	{
+		for (pos[1]=0;pos[1]<numLines[1];++pos[1])
+		{
+			for (pos[0]=0;pos[0]<numLines[0];++pos[0])
+			{
+				ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
+				coord[0] = discLines[0][pos[0]];
+				coord[1] = discLines[1][pos[1]];
+				coord[2] = discLines[2][pos[2]];
+				CSProperties* prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::ELECTRODE);
+				if (prop)
+				{
+					CSPropElectrode* elec = prop->ToElectrode();
+					if (elec->GetType()==0)
+					{
+						vIndex.push_back(ipos);
+						for (int n=0;n<3;++n)
+						{
+							double delta=MainOp->GetIndexDelta(n,pos[n]);
+							vExcit[n].push_back(elec->GetWeightedExcitation(n,coord)*delta);
+						}
+					}
+				}
+			}
+		}
+	}
+	cerr << "size:" << vIndex.size() << endl;
+}
diff --git a/FDTD/cartoperator.h b/FDTD/cartoperator.h
new file mode 100644
index 0000000..67d483b
--- /dev/null
+++ b/FDTD/cartoperator.h
@@ -0,0 +1,63 @@
+#ifndef CARTOPERATOR_H
+#define CARTOPERATOR_H
+
+#include "ContinuousStructure.h"
+#include "tools/AdrOp.h"
+#include "tools/constants.h"
+
+#define FDTD_FLOAT float
+
+class CartOperator
+{
+public:
+    CartOperator();
+    virtual ~CartOperator();
+
+    void SetGeometryCSX(ContinuousStructure* geo);
+
+	int CalcECOperator();
+
+	void ApplyElectricBC(bool* dirs); //applied by default to all boundaries
+	void ApplyMagneticBC(bool* dirs);
+
+	double GetTimestep() {return dT;};
+
+	/*!
+	  Get the voltage excitations. Returns number of excitations, listed position in index, with amplitude in all 3 directions and a possible time delay.
+	  */
+	unsigned int GetVoltageExcitation(unsigned int* &index, FDTD_FLOAT** &excit_amp, FDTD_FLOAT* &excit_delay);
+
+    void Reset();
+
+protected:
+    void Init();
+    ContinuousStructure* CSX;
+
+    AdrOp* MainOp;
+    AdrOp* DualOp;
+    double* discLines[3];
+    unsigned int numLines[3];
+	double gridDelta;
+
+	double dT; //FDTD timestep!
+
+	//EC operator
+	FDTD_FLOAT* vv[3]; //calc new voltage from old voltage
+	FDTD_FLOAT* vi[3]; //calc new voltage from old current
+	FDTD_FLOAT* ii[3]; //calc new current from old current
+	FDTD_FLOAT* iv[3]; //calc new current from old voltage
+
+	//Calc timestep only internal use
+	double CalcTimestep();
+
+	//EC elements, internal only!
+    bool Calc_EC();
+	bool Calc_ECPos(int n, unsigned int* pos, double* inEC);
+	bool Calc_EffMatPos(int n, unsigned int* pos, double* inMat);
+	double* EC_C[3];
+    double* EC_G[3];
+    double* EC_L[3];
+    double* EC_R[3];
+};
+
+#endif // CARTOPERATOR_H
diff --git a/main.cpp b/main.cpp
new file mode 100644
index 0000000..63bdfd7
--- /dev/null
+++ b/main.cpp
@@ -0,0 +1,34 @@
+#include <iostream>
+#include <time.h>
+#include "FDTD/cartoperator.h"
+#include "ContinuousStructure.h"
+
+using namespace std;
+
+int main(int argc, char *argv[])
+{
+    cout << "Hello World" << endl;
+	fprintf(stderr,"%e\n",1.4456);
+
+	time_t startTime=time(NULL);
+
+    ContinuousStructure CSX;
+	CSX.ReadFromXML("csx-files/1Mill.xml");
+
+	CartOperator cop;
+	cop.SetGeometryCSX(&CSX);
+
+	cop.CalcECOperator();
+
+	unsigned int* index = NULL;
+	FDTD_FLOAT** amp=NULL;
+	FDTD_FLOAT* delay=NULL;
+
+	unsigned int nEx = cop.GetVoltageExcitation(index,amp,delay);
+
+	time_t OpDoneTime=time(NULL);
+
+	cerr << "Time for operator: " << difftime(OpDoneTime,startTime) << endl;
+
+    return 0;
+}
diff --git a/openEMS.pro b/openEMS.pro
new file mode 100644
index 0000000..0158306
--- /dev/null
+++ b/openEMS.pro
@@ -0,0 +1,20 @@
+# -------------------------------------------------
+# Project created by QtCreator 2010-02-26T22:34:51
+# -------------------------------------------------
+QT -= gui
+TARGET = openEMS
+CONFIG += console
+CONFIG -= app_bundle
+TEMPLATE = app
+OBJECTS_DIR = obj
+INCLUDEPATH += ../CSXCAD
+LIBS += -L../CSXCAD \
+    -lCSXCAD
+SOURCES += main.cpp \
+    FDTD/cartoperator.cpp \
+    tools/ErrorMsg.cpp \
+    tools/AdrOp.cpp
+HEADERS += FDTD/cartoperator.h \
+    tools/ErrorMsg.h \
+    tools/AdrOp.h \
+    tools/constants.h
diff --git a/tools/AdrOp.cpp b/tools/AdrOp.cpp
new file mode 100644
index 0000000..d009ffd
--- /dev/null
+++ b/tools/AdrOp.cpp
@@ -0,0 +1,503 @@
+#include "AdrOp.h"
+
+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... \0"); exit(1); }
+	error->SetMsg(1,"Adress Operator: Memory allocation failed!! exiting... \0");
+	error->SetMsg(2,"Adress Operator: Invalid Adress requested!! exiting... \0");
+	error->SetMsg(3,"Adress Operator: Invalid Position set!! exiting...\0");
+	error->SetMsg(4,"Adress Operator: Invalid jump or passing end of iteration!! exiting...\0");
+	error->SetMsg(5,"Adress Operator: 4D not yet implemented!! exiting... \0");
+	error->SetMsg(6,"Adress Operator: Position not set!! exiting... \0");
+	error->SetMsg(7,"Adress Operator: Cells not added to Adress Operator!! exiting... \0");
+	error->SetMsg(8,"Adress Operator: Invalid Node!! exiting... \0");
+	error->SetMsg(9,"Adress Operator: Grid invalid!! exiting... \0");
+	
+	//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()
+{
+	return GetPos(iIshift,iJshift,iKshift);
+}
+
+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]={uiImax-1,uiJmax-1,uiKmax-1,uiLmax-1};
+	while (ny<0) ny+=uiDimension;
+	ny=ny%uiDimension;
+	iCellShift[ny]=step;
+	int uiPos[4]={uiIpos+iCellShift[0],uiJpos+iCellShift[1],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/tools/AdrOp.h b/tools/AdrOp.h
new file mode 100644
index 0000000..876fde0
--- /dev/null
+++ b/tools/AdrOp.h
@@ -0,0 +1,133 @@
+/*!
+\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"
+
+using namespace std;
+
+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"
+	unsigned int GetShiftedPos();
+	///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/tools/ErrorMsg.cpp b/tools/ErrorMsg.cpp
new file mode 100644
index 0000000..5fd67a1
--- /dev/null
+++ b/tools/ErrorMsg.cpp
@@ -0,0 +1,66 @@
+#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... \0"); 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, 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... \0"); 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/tools/ErrorMsg.h b/tools/ErrorMsg.h
new file mode 100644
index 0000000..d02c5b5
--- /dev/null
+++ b/tools/ErrorMsg.h
@@ -0,0 +1,33 @@
+/*!
+\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, 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/tools/ExpenseLog.cpp b/tools/ExpenseLog.cpp
new file mode 100644
index 0000000..1ddb7a7
--- /dev/null
+++ b/tools/ExpenseLog.cpp
@@ -0,0 +1,141 @@
+#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/tools/ExpenseLog.h b/tools/ExpenseLog.h
new file mode 100644
index 0000000..14613e9
--- /dev/null
+++ b/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/tools/constants.h b/tools/constants.h
new file mode 100644
index 0000000..d88fa6a
--- /dev/null
+++ b/tools/constants.h
@@ -0,0 +1,7 @@
+#ifndef CONSTANTS_H
+#define CONSTANTS_H
+
+#define __EPS0__ 8.85418781762e-12
+#define __MUE0__ 1.256637062e-6
+
+#endif // CONSTANTS_H