Imported Upstream version 0.0.34

1 files changed, 2131 insertions, 0 deletions

diff --git a/openEMS/FDTD/operator.cpp b/openEMS/FDTD/operator.cpp
new file mode 100644
index 0000000..a7582aa
--- /dev/null
+++ b/openEMS/FDTD/operator.cpp
@@ -0,0 +1,2131 @@
+/*
+*	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 <fstream>
+#include <algorithm>
+#include "operator.h"
+#include "engine.h"
+#include "extensions/operator_extension.h"
+#include "extensions/operator_ext_excitation.h"
+#include "Common/processfields.h"
+#include "tools/array_ops.h"
+#include "tools/vtk_file_writer.h"
+#include "fparser.hh"
+#include "extensions/operator_ext_excitation.h"
+
+#include "vtkPolyData.h"
+#include "vtkCellArray.h"
+#include "vtkPoints.h"
+#include "vtkXMLPolyDataWriter.h"
+#include "CSPrimBox.h"
+#include "CSPrimCurve.h"
+
+#include "CSPropMaterial.h"
+#include "CSPropLumpedElement.h"
+
+Operator* Operator::New()
+{
+	cout << "Create FDTD operator" << endl;
+	Operator* op = new Operator();
+	op->Init();
+	return op;
+}
+
+Operator::Operator() : Operator_Base()
+{
+	m_Exc = 0;
+	m_InvaildTimestep = false;
+	m_TimeStepVar = 3;
+}
+
+Operator::~Operator()
+{
+	for (size_t n=0; n<m_Op_exts.size(); ++n)
+		delete m_Op_exts.at(n);
+	m_Op_exts.clear();
+
+	Delete();
+}
+
+Engine* Operator::CreateEngine()
+{
+	m_Engine = Engine::New(this);
+	return m_Engine;
+}
+
+void Operator::Init()
+{
+	CSX = NULL;
+	m_Engine = NULL;
+
+	Operator_Base::Init();
+
+	vv=NULL;
+	vi=NULL;
+	iv=NULL;
+	ii=NULL;
+
+	m_epsR=NULL;
+	m_kappa=NULL;
+	m_mueR=NULL;
+	m_sigma=NULL;
+
+	MainOp=NULL;
+
+	for (int n=0; n<3; ++n)
+	{
+		EC_C[n]=NULL;
+		EC_G[n]=NULL;
+		EC_L[n]=NULL;
+		EC_R[n]=NULL;
+	}
+
+	m_Exc = 0;
+	m_TimeStepFactor = 1;
+	SetMaterialAvgMethod(QuarterCell);
+}
+
+void Operator::Delete()
+{
+	CSX = NULL;
+
+	Delete_N_3DArray(vv,numLines);
+	Delete_N_3DArray(vi,numLines);
+	Delete_N_3DArray(iv,numLines);
+	Delete_N_3DArray(ii,numLines);
+	vv=vi=iv=ii=0;
+	delete MainOp; MainOp=0;
+	for (int n=0; n<3; ++n)
+	{
+		delete[] EC_C[n];EC_C[n]=0;
+		delete[] EC_G[n];EC_G[n]=0;
+		delete[] EC_L[n];EC_L[n]=0;
+		delete[] EC_R[n];EC_R[n]=0;
+	}
+
+	Delete_N_3DArray(m_epsR,numLines);
+	m_epsR=0;
+	Delete_N_3DArray(m_kappa,numLines);
+	m_kappa=0;
+	Delete_N_3DArray(m_mueR,numLines);
+	m_mueR=0;
+	Delete_N_3DArray(m_sigma,numLines);
+	m_sigma=0;
+}
+
+void Operator::Reset()
+{
+	Delete();
+	Operator_Base::Reset();
+}
+
+double Operator::GetDiscLine(int n, unsigned int pos, bool dualMesh) const
+{
+	if ((n<0) || (n>2)) return 0.0;
+	if (pos>=numLines[n]) return 0.0;
+	if (dualMesh==false)
+		return discLines[n][pos];
+
+	// return dual mesh node
+	if (pos<numLines[n]-1)
+		return 0.5*(discLines[n][pos] + discLines[n][pos+1]);
+
+	// dual node for the last line (outside the field domain)
+	return discLines[n][pos] + 0.5*(discLines[n][pos] - discLines[n][pos-1]);
+}
+
+double Operator::GetDiscDelta(int n, unsigned int pos, bool dualMesh) const
+{
+	if ((n<0) || (n>2)) return 0.0;
+	if (pos>=numLines[n]) return 0.0;
+	double delta=0;
+	if (dualMesh==false)
+	{
+		if (pos<numLines[n]-1)
+			delta = GetDiscLine(n,pos+1,false) - GetDiscLine(n,pos,false);
+		else
+			delta = GetDiscLine(n,pos,false) - GetDiscLine(n,pos-1,false);
+		return delta;
+	}
+	else
+	{
+		if (pos>0)
+			delta = GetDiscLine(n,pos,true) - GetDiscLine(n,pos-1,true);
+		else
+			delta = GetDiscLine(n,1,false) - GetDiscLine(n,0,false);
+		return delta;
+	}
+}
+
+bool Operator::GetYeeCoords(int ny, unsigned int pos[3], double* coords, bool dualMesh) const
+{
+	for (int n=0;n<3;++n)
+		coords[n]=GetDiscLine(n,pos[n],dualMesh);
+	coords[ny]=GetDiscLine(ny,pos[ny],!dualMesh);
+
+	//check if position is inside the FDTD domain
+	if (dualMesh==false) //main grid
+	{
+		if (pos[ny]>=numLines[ny]-1)
+			return false;
+	}
+	else	//dual grid
+	{
+		int nP = (ny+1)%3;
+		int nPP = (ny+2)%3;
+		if ((pos[nP]>=numLines[nP]-1) || (pos[nPP]>=numLines[nPP]-1))
+			return false;
+	}
+	return true;
+}
+
+bool Operator::GetNodeCoords(const unsigned int pos[3], double* coords, bool dualMesh, CoordinateSystem c_system) const
+{
+	for (int n=0;n<3;++n)
+		coords[n]=GetDiscLine(n,pos[n],dualMesh);
+	TransformCoordSystem(coords,coords,m_MeshType,c_system);
+	return true;
+}
+
+double Operator::GetEdgeLength(int n, const unsigned int* pos, bool dualMesh) const
+{
+	return GetDiscDelta(n,pos[n],dualMesh)*gridDelta;
+}
+
+double Operator::GetCellVolume(const unsigned int pos[3], bool dualMesh) const
+{
+	double vol=1;
+	for (int n=0;n<3;++n)
+		vol*=GetEdgeLength(n,pos,dualMesh);
+	return vol;
+}
+
+double Operator::GetNodeWidth(int ny, const int pos[3], bool dualMesh) const
+{
+	if ( (pos[0]<0) || (pos[1]<0) || (pos[2]<0) )
+		return 0.0;
+
+	//call the unsigned int version of GetNodeWidth
+	unsigned int uiPos[]={(unsigned int)pos[0],(unsigned int)pos[1],(unsigned int)pos[2]};
+	return GetNodeWidth(ny, uiPos, dualMesh);
+}
+
+double Operator::GetNodeArea(int ny, const unsigned int pos[3], bool dualMesh) const
+{
+	int nyP = (ny+1)%3;
+	int nyPP = (ny+2)%3;
+	return GetNodeWidth(nyP,pos,dualMesh) * GetNodeWidth(nyPP,pos,dualMesh);
+}
+
+double Operator::GetNodeArea(int ny, const int pos[3], bool dualMesh) const
+{
+	if ( (pos[0]<0) || (pos[1]<0) || (pos[2]<0) )
+		return 0.0;
+
+	//call the unsigned int version of GetNodeArea
+	unsigned int uiPos[]={(unsigned int)pos[0],(unsigned int)pos[1],(unsigned int)pos[2]};
+	return GetNodeArea(ny, uiPos, dualMesh);
+}
+
+unsigned int Operator::SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh, bool fullMesh) const
+{
+	inside = false;
+	if ((ny<0) || (ny>2))
+		return 0;
+	if (coord<GetDiscLine(ny,0))
+		return 0;
+	unsigned int numLines = GetNumberOfLines(ny, fullMesh);
+	if (coord>GetDiscLine(ny,numLines-1))
+		return numLines-1;
+	inside=true;
+	if (dualMesh==false)
+	{
+		for (unsigned int n=0;n<numLines;++n)
+		{
+			if (coord<=GetDiscLine(ny,n,true))
+				return n;
+		}
+	}
+	else
+	{
+		for (unsigned int n=1;n<numLines;++n)
+		{
+			if (coord<=GetDiscLine(ny,n,false))
+				return n-1;
+		}
+	}
+	//should not happen
+	return 0;
+}
+
+bool Operator::SnapToMesh(const double* dcoord, unsigned int* uicoord, bool dualMesh, bool fullMesh, bool* inside) const
+{
+	bool meshInside=false;
+	bool ok=true;
+	for (int n=0; n<3; ++n)
+	{
+		uicoord[n] = SnapToMeshLine(n,dcoord[n],meshInside,dualMesh,fullMesh);
+		ok &= meshInside;
+		if (inside)
+			inside[n]=meshInside;
+	}
+	return ok;
+}
+
+int Operator::SnapBox2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh, bool fullMesh, int SnapMethod, bool* bStartIn, bool* bStopIn) const
+{
+	double l_start[3], l_stop[3];
+	for (int n=0;n<3;++n)
+	{
+		l_start[n] = fmin(start[n],stop[n]);
+		l_stop[n] = fmax(start[n], stop[n]);
+		double min = GetDiscLine(n,0);
+		double max = GetDiscLine(n,GetNumberOfLines(n, fullMesh)-1);
+		if ( ((l_start[n]<min) && (l_stop[n]<min)) || ((l_start[n]>max) && (l_stop[n]>max)) )
+		{
+			return -2;
+		}
+	}
+
+	SnapToMesh(l_start, uiStart, dualMesh, fullMesh, bStartIn);
+	SnapToMesh(l_stop, uiStop, dualMesh, fullMesh, bStopIn);
+	int iDim = 0;
+
+	if (SnapMethod==0)
+	{
+		for (int n=0;n<3;++n)
+			if (uiStop[n]>uiStart[n])
+				++iDim;
+		return iDim;
+	}
+	else if (SnapMethod==1)
+	{
+		for (int n=0;n<3;++n)
+		{
+			if (uiStop[n]>uiStart[n])
+			{
+				if ((GetDiscLine( n, uiStart[n], dualMesh ) > l_start[n]) && (uiStart[n]>0))
+					--uiStart[n];
+				if ((GetDiscLine( n, uiStop[n], dualMesh ) < l_stop[n]) && (uiStop[n]<GetNumberOfLines(n, fullMesh)-1))
+					++uiStop[n];
+			}
+			if (uiStop[n]>uiStart[n])
+				++iDim;
+		}
+		return iDim;
+	}
+	else if (SnapMethod==2)
+	{
+		for (int n=0;n<3;++n)
+		{
+			if (uiStop[n]>uiStart[n])
+			{
+				if ((GetDiscLine( n, uiStart[n], dualMesh ) < l_start[n]) && (uiStart[n]<GetNumberOfLines(n, fullMesh)-1))
+					++uiStart[n];
+				if ((GetDiscLine( n, uiStop[n], dualMesh ) > l_stop[n]) && (uiStop[n]>0))
+					--uiStop[n];
+			}
+			if (uiStop[n]>uiStart[n])
+				++iDim;
+		}
+		return iDim;
+	}
+	else
+		cerr << "Operator::SnapBox2Mesh: Unknown snapping method!" << endl;
+	return -1;
+}
+
+int Operator::SnapLine2Mesh(const double* start, const double* stop, unsigned int* uiStart, unsigned int* uiStop, bool dualMesh, bool fullMesh) const
+{
+	bool bStartIn[3];
+	bool bStopIn[3];
+	SnapToMesh(start, uiStart, dualMesh, fullMesh, bStartIn);
+	SnapToMesh(stop, uiStop, dualMesh, fullMesh, bStopIn);
+
+	for (int n=0;n<3;++n)
+	{
+		if ((start[n]<GetDiscLine(n,0)) && (stop[n]<GetDiscLine(n,0)))
+			return -1; //lower bound violation
+		if ((start[n]>GetDiscLine(n,GetNumberOfLines(n,true)-1)) && (stop[n]>GetDiscLine(n,GetNumberOfLines(n,true)-1)))
+			return -1; //upper bound violation
+	}
+
+	int ret = 0;
+	if (!(bStartIn[0] && bStartIn[1] && bStartIn[2]))
+		ret = ret + 1;
+	if (!(bStopIn[0] && bStopIn[1] && bStopIn[2]))
+		ret = ret + 2;
+	if (ret==0)
+		return ret;
+
+	//fixme, do we need to do something about start or stop being outside the field domain?
+	//maybe caclulate the intersection point and snap to that?
+	//it seems to work like this as well...
+
+	return ret;
+}
+
+
+Grid_Path Operator::FindPath(double start[], double stop[])
+{
+	Grid_Path path;
+	unsigned int uiStart[3],uiStop[3],currPos[3];
+
+	int ret = SnapLine2Mesh(start, stop, uiStart, uiStop, false, true);
+	if (ret<0)
+		return path;
+
+	currPos[0]=uiStart[0];
+	currPos[1]=uiStart[1];
+	currPos[2]=uiStart[2];
+	double meshStart[3] = {discLines[0][uiStart[0]], discLines[1][uiStart[1]], discLines[2][uiStart[2]]};
+	double meshStop[3] = {discLines[0][uiStop[0]], discLines[1][uiStop[1]], discLines[2][uiStop[2]]};
+	bool UpDir = false;
+	double foot=0,dist=0,minFoot=0,minDist=0;
+	int minDir=0;
+	unsigned int minPos[3];
+	double startFoot,stopFoot,currFoot;
+	Point_Line_Distance(meshStart,start,stop,startFoot,dist, m_MeshType);
+	Point_Line_Distance(meshStop,start,stop,stopFoot,dist, m_MeshType);
+	currFoot=startFoot;
+	minFoot=startFoot;
+	double P[3];
+
+	while (minFoot<stopFoot)
+	{
+		minDist=1e300;
+		for (int n=0; n<3; ++n) //check all 6 surrounding points
+		{
+			P[0] = discLines[0][currPos[0]];
+			P[1] = discLines[1][currPos[1]];
+			P[2] = discLines[2][currPos[2]];
+			if (((int)currPos[n]-1)>=0)
+			{
+				P[n] = discLines[n][currPos[n]-1];
+				Point_Line_Distance(P,start,stop,foot,dist, m_MeshType);
+				if ((foot>currFoot) && (dist<minDist))
+				{
+					minFoot=foot;
+					minDist=dist;
+					minDir = n;
+					UpDir = false;
+				}
+			}
+			if ((currPos[n]+1)<numLines[n])
+			{
+				P[n] = discLines[n][currPos[n]+1];
+				Point_Line_Distance(P,start,stop,foot,dist, m_MeshType);
+				if ((foot>currFoot) && (dist<minDist))
+				{
+					minFoot=foot;
+					minDist=dist;
+					minDir = n;
+					UpDir = true;
+				}
+			}
+		}
+		minPos[0]=currPos[0];
+		minPos[1]=currPos[1];
+		minPos[2]=currPos[2];
+		if (UpDir)
+		{
+			currPos[minDir]+=1;
+		}
+		else
+		{
+			currPos[minDir]+=-1;
+			minPos[minDir]-=1;
+		}
+		//check validity of current postion
+		for (int n=0;n<3;++n)
+			if (currPos[n]>=numLines[n])
+			{
+				cerr << __func__ << ": Error, path went out of simulation domain, skipping path!" << endl;
+				Grid_Path empty;
+				return empty;
+			}
+		path.posPath[0].push_back(minPos[0]);
+		path.posPath[1].push_back(minPos[1]);
+		path.posPath[2].push_back(minPos[2]);
+		currFoot=minFoot;
+		path.dir.push_back(minDir);
+	}
+
+	//close missing edges, if currPos is not equal to uiStopPos
+	for (int n=0; n<3; ++n)
+	{
+		if (currPos[n]>uiStop[n])
+		{
+			--currPos[n];
+			path.posPath[0].push_back(currPos[0]);
+			path.posPath[1].push_back(currPos[1]);
+			path.posPath[2].push_back(currPos[2]);
+			path.dir.push_back(n);
+		}
+		else if (currPos[n]<uiStop[n])
+		{
+			path.posPath[0].push_back(currPos[0]);
+			path.posPath[1].push_back(currPos[1]);
+			path.posPath[2].push_back(currPos[2]);
+			path.dir.push_back(n);
+		}
+	}
+
+	return path;
+}
+
+void Operator::SetMaterialAvgMethod(MatAverageMethods method)
+{
+	switch (method)
+	{
+	default:
+	case QuarterCell:
+		return SetQuarterCellMaterialAvg();
+	case CentralCell:
+		return SetCellConstantMaterial();
+	}
+}
+
+double Operator::GetNumberCells() const
+{
+	if (numLines)
+		return (numLines[0])*(numLines[1])*(numLines[2]); //it's more like number of nodes???
+	return 0;
+}
+
+void Operator::ShowStat() const
+{
+	unsigned int OpSize = 12*numLines[0]*numLines[1]*numLines[2]*sizeof(FDTD_FLOAT);
+	unsigned int FieldSize = 6*numLines[0]*numLines[1]*numLines[2]*sizeof(FDTD_FLOAT);
+	double MBdiff = 1024*1024;
+
+	cout << "------- Stat: FDTD Operator -------" << endl;
+	cout << "Dimensions\t\t: " << numLines[0] << "x" << numLines[1] << "x" << numLines[2] << " = " <<  numLines[0]*numLines[1]*numLines[2] << " Cells (" << numLines[0]*numLines[1]*numLines[2]/1e6 << " MCells)" << endl;
+	cout << "Size of Operator\t: " << OpSize << " Byte (" << (double)OpSize/MBdiff << " MiB) " << endl;
+	cout << "Size of Field-Data\t: " << FieldSize << " Byte (" << (double)FieldSize/MBdiff << " MiB) " << endl;
+	cout << "-----------------------------------" << endl;
+	cout << "Background materials (epsR/mueR/kappa/sigma): " << GetBackgroundEpsR() << "/" << GetBackgroundMueR() << "/" << GetBackgroundKappa() << "/" << GetBackgroundSigma() << endl;
+	cout << "-----------------------------------" << endl;
+	cout << "Number of PEC edges\t: " << m_Nr_PEC[0]+m_Nr_PEC[1]+m_Nr_PEC[2] << endl;
+	cout << "in " << GetDirName(0) << " direction\t\t: " << m_Nr_PEC[0] << endl;
+	cout << "in " << GetDirName(1) << " direction\t\t: " << m_Nr_PEC[1] << endl;
+	cout << "in " << GetDirName(2) << " direction\t\t: " << m_Nr_PEC[2] << endl;
+	cout << "-----------------------------------" << endl;
+	cout << "Timestep (s)\t\t: " << dT ;
+	if (opt_dT)
+		cout <<"\t(" << opt_dT << ")";
+	cout << endl;
+	cout << "Timestep method name\t: " << m_Used_TS_Name << endl;
+	cout << "Nyquist criteria (TS)\t: " << m_Exc->GetNyquistNum() << endl;
+	cout << "Nyquist criteria (s)\t: " << m_Exc->GetNyquistNum()*dT << endl;
+	cout << "-----------------------------------" << endl;
+}
+
+void Operator::ShowExtStat() const
+{
+	if (m_Op_exts.size()==0) return;
+	cout << "-----------------------------------" << endl;
+	for (size_t n=0; n<m_Op_exts.size(); ++n)
+		m_Op_exts.at(n)->ShowStat(cout);
+	cout << "-----------------------------------" << endl;
+}
+
+void Operator::DumpOperator2File(string filename)
+{
+#ifdef OUTPUT_IN_DRAWINGUNITS
+	double discLines_scaling = 1;
+#else
+	double discLines_scaling = GetGridDelta();
+#endif
+
+	cout << "Operator: Dumping FDTD operator information to vtk file: " << filename << " ..." << flush;
+
+	VTK_File_Writer* vtk_Writer = new VTK_File_Writer(filename.c_str(), m_MeshType);
+	vtk_Writer->SetMeshLines(discLines,numLines,discLines_scaling);
+	vtk_Writer->SetHeader("openEMS - Operator dump");
+
+	vtk_Writer->SetNativeDump(true);
+
+	//find excitation extension
+	Operator_Ext_Excitation* Op_Ext_Exc=GetExcitationExtension();
+
+	if (Op_Ext_Exc)
+	{
+		FDTD_FLOAT**** exc = NULL;
+		if (Op_Ext_Exc->Volt_Count>0)
+		{
+			exc = Create_N_3DArray<FDTD_FLOAT>(numLines);
+			for (unsigned int n=0; n<  Op_Ext_Exc->Volt_Count; ++n)
+				exc[  Op_Ext_Exc->Volt_dir[n]][  Op_Ext_Exc->Volt_index[0][n]][  Op_Ext_Exc->Volt_index[1][n]][  Op_Ext_Exc->Volt_index[2][n]] =   Op_Ext_Exc->Volt_amp[n];
+			vtk_Writer->AddVectorField("exc_volt",exc);
+			Delete_N_3DArray(exc,numLines);
+		}
+
+		if (  Op_Ext_Exc->Curr_Count>0)
+		{
+			exc = Create_N_3DArray<FDTD_FLOAT>(numLines);
+			for (unsigned int n=0; n<  Op_Ext_Exc->Curr_Count; ++n)
+				exc[  Op_Ext_Exc->Curr_dir[n]][  Op_Ext_Exc->Curr_index[0][n]][  Op_Ext_Exc->Curr_index[1][n]][  Op_Ext_Exc->Curr_index[2][n]] =   Op_Ext_Exc->Curr_amp[n];
+			vtk_Writer->AddVectorField("exc_curr",exc);
+			Delete_N_3DArray(exc,numLines);
+		}
+	}
+
+	FDTD_FLOAT**** vv_temp = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	FDTD_FLOAT**** vi_temp = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	FDTD_FLOAT**** iv_temp = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	FDTD_FLOAT**** ii_temp = Create_N_3DArray<FDTD_FLOAT>(numLines);
+
+	unsigned int pos[3], n;
+	for (n=0; n<3; n++)
+		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]++)
+				{
+					vv_temp[n][pos[0]][pos[1]][pos[2]] = GetVV(n,pos);
+					vi_temp[n][pos[0]][pos[1]][pos[2]] = GetVI(n,pos);
+					iv_temp[n][pos[0]][pos[1]][pos[2]] = GetIV(n,pos);
+					ii_temp[n][pos[0]][pos[1]][pos[2]] = GetII(n,pos);
+				}
+
+
+	vtk_Writer->AddVectorField("vv",vv_temp);
+	Delete_N_3DArray(vv_temp,numLines);
+	vtk_Writer->AddVectorField("vi",vi_temp);
+	Delete_N_3DArray(vi_temp,numLines);
+	vtk_Writer->AddVectorField("iv",iv_temp);
+	Delete_N_3DArray(iv_temp,numLines);
+	vtk_Writer->AddVectorField("ii",ii_temp);
+	Delete_N_3DArray(ii_temp,numLines);
+
+	if (vtk_Writer->Write()==false)
+		cerr << "Operator::DumpOperator2File: Error: Can't write file... skipping!" << endl;
+
+	delete vtk_Writer;
+}
+
+//! \brief dump PEC (perfect electric conductor) information (into VTK-file)
+//! visualization via paraview
+//! visualize only one component (x, y or z)
+void Operator::DumpPEC2File(string filename , unsigned int *range)
+{
+	cout << "Operator: Dumping PEC information to vtk file: " << filename << " ..." << flush;
+
+#ifdef OUTPUT_IN_DRAWINGUNITS
+	double scaling = 1.0;
+#else
+	double scaling = GetGridDelta();;
+#endif
+
+	unsigned int start[3] = {0, 0, 0};
+	unsigned int stop[3] = {numLines[0]-1,numLines[1]-1,numLines[2]-1};
+
+	if (range!=NULL)
+		for (int n=0;n<3;++n)
+		{
+			start[n] = range[2*n];
+			stop[n] = range[2*n+1];
+		}
+
+	vtkPolyData* polydata = vtkPolyData::New();
+	vtkCellArray *poly = vtkCellArray::New();
+	vtkPoints *points = vtkPoints::New();
+
+	int* pointIdx[2];
+	pointIdx[0] = new int[numLines[0]*numLines[1]];
+	pointIdx[1] = new int[numLines[0]*numLines[1]];
+	// init point idx
+	for (unsigned int n=0;n<numLines[0]*numLines[1];++n)
+	{
+		pointIdx[0][n]=-1;
+		pointIdx[1][n]=-1;
+	}
+
+	int nP,nPP;
+	double coord[3];
+	unsigned int pos[3],rpos[3];
+	unsigned int mesh_idx=0;
+	for (pos[2]=start[2];pos[2]<stop[2];++pos[2])
+	{ // each xy-plane
+		for (unsigned int n=0;n<numLines[0]*numLines[1];++n)
+		{
+			pointIdx[0][n]=pointIdx[1][n];
+			pointIdx[1][n]=-1;
+		}
+		for (pos[0]=start[0];pos[0]<stop[0];++pos[0])
+			for (pos[1]=start[1];pos[1]<stop[1];++pos[1])
+			{
+				for (int n=0;n<3;++n)
+				{
+					nP = (n+1)%3;
+					nPP = (n+2)%3;
+					if ((GetVV(n,pos) == 0) && (GetVI(n,pos) == 0) && (pos[nP]>0) && (pos[nPP]>0))
+					{
+						rpos[0]=pos[0];
+						rpos[1]=pos[1];
+						rpos[2]=pos[2];
+
+						poly->InsertNextCell(2);
+
+						mesh_idx = rpos[0] + rpos[1]*numLines[0];
+						if (pointIdx[0][mesh_idx]<0)
+						{
+							for (int m=0;m<3;++m)
+								coord[m] = discLines[m][rpos[m]];
+							TransformCoordSystem(coord, coord, m_MeshType, CARTESIAN);
+							for (int m=0;m<3;++m)
+								coord[m] *= scaling;
+							pointIdx[0][mesh_idx] = (int)points->InsertNextPoint(coord);
+						}
+						poly->InsertCellPoint(pointIdx[0][mesh_idx]);
+
+						++rpos[n];
+						mesh_idx = rpos[0] + rpos[1]*numLines[0];
+						if (pointIdx[n==2][mesh_idx]<0)
+						{
+							for (int m=0;m<3;++m)
+								coord[m] = discLines[m][rpos[m]];
+							TransformCoordSystem(coord, coord, m_MeshType, CARTESIAN);
+							for (int m=0;m<3;++m)
+								coord[m] *= scaling;
+							pointIdx[n==2][mesh_idx] = (int)points->InsertNextPoint(coord);
+						}
+						poly->InsertCellPoint(pointIdx[n==2][mesh_idx]);
+					}
+				}
+			}
+	}
+	delete[] pointIdx[0];
+	delete[] pointIdx[1];
+
+	polydata->SetPoints(points);
+	points->Delete();
+	polydata->SetLines(poly);
+	poly->Delete();
+
+	vtkXMLPolyDataWriter* writer  = vtkXMLPolyDataWriter::New();
+	filename += ".vtp";
+	writer->SetFileName(filename.c_str());
+
+#if VTK_MAJOR_VERSION>=6
+	writer->SetInputData(polydata);
+#else
+	writer->SetInput(polydata);
+#endif
+	writer->Write();
+
+	writer->Delete();
+	polydata->Delete();
+	cout << " done." << endl;
+}
+
+void Operator::DumpMaterial2File(string filename)
+{
+#ifdef OUTPUT_IN_DRAWINGUNITS
+	double discLines_scaling = 1;
+#else
+	double discLines_scaling = GetGridDelta();
+#endif
+
+	cout << "Operator: Dumping material information to vtk file: " << filename << " ..."  << flush;
+
+	FDTD_FLOAT**** epsilon = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	FDTD_FLOAT**** mue     = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	FDTD_FLOAT**** kappa   = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	FDTD_FLOAT**** sigma   = Create_N_3DArray<FDTD_FLOAT>(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])
+		{
+			vector<CSPrimitives*> vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL);
+			for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
+			{
+				for (int n=0; n<3; ++n)
+				{
+					double inMat[4];
+					Calc_EffMatPos(n, pos, inMat, vPrims);
+					epsilon[n][pos[0]][pos[1]][pos[2]] = inMat[0]/__EPS0__;
+					mue[n][pos[0]][pos[1]][pos[2]]     = inMat[2]/__MUE0__;
+					kappa[n][pos[0]][pos[1]][pos[2]]   = inMat[1];
+					sigma[n][pos[0]][pos[1]][pos[2]]   = inMat[3];
+				}
+			}
+		}
+	}
+
+	VTK_File_Writer* vtk_Writer = new VTK_File_Writer(filename.c_str(), m_MeshType);
+	vtk_Writer->SetMeshLines(discLines,numLines,discLines_scaling);
+	vtk_Writer->SetHeader("openEMS - material dump");
+
+	vtk_Writer->SetNativeDump(true);
+
+	vtk_Writer->AddVectorField("epsilon",epsilon);
+	Delete_N_3DArray(epsilon,numLines);
+	vtk_Writer->AddVectorField("mue",mue);
+	Delete_N_3DArray(mue,numLines);
+	vtk_Writer->AddVectorField("kappa",kappa);
+	Delete_N_3DArray(kappa,numLines);
+	vtk_Writer->AddVectorField("sigma",sigma);
+	Delete_N_3DArray(sigma,numLines);
+
+	if (vtk_Writer->Write()==false)
+		cerr << "Operator::DumpMaterial2File: Error: Can't write file... skipping!" << endl;
+
+	delete vtk_Writer;
+}
+
+ bool Operator::SetupCSXGrid(CSRectGrid* grid)
+ {
+	 for (int n=0; n<3; ++n)
+	 {
+		 discLines[n] = grid->GetLines(n,discLines[n],numLines[n],true);
+		 if (numLines[n]<3)
+		 {
+			 cerr << "CartOperator::SetupCSXGrid: you need at least 3 disc-lines in every direction (3D!)!!!" << endl;
+			 Reset();
+			 return false;
+		 }
+	 }
+	 MainOp = new AdrOp(numLines[0],numLines[1],numLines[2]);
+	 MainOp->SetGrid(discLines[0],discLines[1],discLines[2]);
+	 if (grid->GetDeltaUnit()<=0)
+	 {
+		 cerr << "CartOperator::SetupCSXGrid: grid delta unit must not be <=0 !!!" << endl;
+		 Reset();
+		 return false;
+	 }
+	 else gridDelta=grid->GetDeltaUnit();
+	 MainOp->SetGridDelta(1);
+	 MainOp->AddCellAdrOp();
+
+	 //delete the grid clone...
+	 delete grid;
+	 return true;
+ }
+
+bool Operator::SetGeometryCSX(ContinuousStructure* geo)
+{
+	if (geo==NULL) return false;
+
+	CSX = geo;
+
+	CSBackgroundMaterial* bg_mat=CSX->GetBackgroundMaterial();
+	SetBackgroundEpsR(bg_mat->GetEpsilon());
+	SetBackgroundMueR(bg_mat->GetMue());
+	SetBackgroundKappa(bg_mat->GetKappa());
+	SetBackgroundSigma(bg_mat->GetSigma());
+	SetBackgroundDensity(0);
+
+	CSRectGrid* grid=CSX->GetGrid();
+	return SetupCSXGrid(CSRectGrid::Clone(grid));
+}
+
+void Operator::InitOperator()
+{
+	Delete_N_3DArray(vv,numLines);
+	Delete_N_3DArray(vi,numLines);
+	Delete_N_3DArray(iv,numLines);
+	Delete_N_3DArray(ii,numLines);
+	vv = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	vi = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	iv = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	ii = Create_N_3DArray<FDTD_FLOAT>(numLines);
+}
+
+void Operator::InitDataStorage()
+{
+	if (m_StoreMaterial[0])
+	{
+		if (g_settings.GetVerboseLevel()>0)
+			cerr << "Operator::InitDataStorage(): Storing epsR material data..." << endl;
+		Delete_N_3DArray(m_epsR,numLines);
+		m_epsR = Create_N_3DArray<float>(numLines);
+	}
+	if (m_StoreMaterial[1])
+	{
+		if (g_settings.GetVerboseLevel()>0)
+			cerr << "Operator::InitDataStorage(): Storing kappa material data..." << endl;
+		Delete_N_3DArray(m_kappa,numLines);
+		m_kappa = Create_N_3DArray<float>(numLines);
+	}
+	if (m_StoreMaterial[2])
+	{
+		if (g_settings.GetVerboseLevel()>0)
+			cerr << "Operator::InitDataStorage(): Storing muR material data..." << endl;
+		Delete_N_3DArray(m_mueR,numLines);
+		m_mueR = Create_N_3DArray<float>(numLines);
+	}
+	if (m_StoreMaterial[3])
+	{
+		if (g_settings.GetVerboseLevel()>0)
+			cerr << "Operator::InitDataStorage(): Storing sigma material data..." << endl;
+		Delete_N_3DArray(m_sigma,numLines);
+		m_sigma = Create_N_3DArray<float>(numLines);
+	}
+}
+
+void Operator::CleanupMaterialStorage()
+{
+	if (!m_StoreMaterial[0] && m_epsR)
+	{
+		if (g_settings.GetVerboseLevel()>0)
+			cerr << "Operator::CleanupMaterialStorage(): Delete epsR material data..." << endl;
+		Delete_N_3DArray(m_epsR,numLines);
+		m_epsR = NULL;
+	}
+	if (!m_StoreMaterial[1] && m_kappa)
+	{
+		if (g_settings.GetVerboseLevel()>0)
+			cerr << "Operator::CleanupMaterialStorage(): Delete kappa material data..." << endl;
+		Delete_N_3DArray(m_kappa,numLines);
+		m_kappa = NULL;
+	}
+	if (!m_StoreMaterial[2] && m_mueR)
+	{
+		if (g_settings.GetVerboseLevel()>0)
+			cerr << "Operator::CleanupMaterialStorage(): Delete mueR material data..." << endl;
+		Delete_N_3DArray(m_mueR,numLines);
+		m_mueR = NULL;
+	}
+	if (!m_StoreMaterial[3] && m_sigma)
+	{
+		if (g_settings.GetVerboseLevel()>0)
+			cerr << "Operator::CleanupMaterialStorage(): Delete sigma material data..." << endl;
+		Delete_N_3DArray(m_sigma,numLines);
+		m_sigma = NULL;
+	}
+}
+
+double Operator::GetDiscMaterial(int type, int n, const unsigned int pos[3]) const
+{
+	switch (type)
+	{
+	case 0:
+		if (m_epsR==0)
+			return 0;
+		return m_epsR[n][pos[0]][pos[1]][pos[2]];
+	case 1:
+		if (m_kappa==0)
+			return 0;
+		return m_kappa[n][pos[0]][pos[1]][pos[2]];
+	case 2:
+		if (m_mueR==0)
+			return 0;
+		return m_mueR[n][pos[0]][pos[1]][pos[2]];
+	case 3:
+		if (m_sigma==0)
+			return 0;
+		return m_sigma[n][pos[0]][pos[1]][pos[2]];
+	}
+	return 0;
+}
+
+void Operator::SetExcitationSignal(Excitation* exc)
+{
+	m_Exc=exc;
+}
+
+void Operator::Calc_ECOperatorPos(int n, unsigned int* pos)
+{
+	unsigned int i = MainOp->SetPos(pos[0],pos[1],pos[2]);
+	double C = EC_C[n][i];
+	double G = EC_G[n][i];
+	if (C>0)
+	{
+		SetVV(n,pos[0],pos[1],pos[2], (1.0-dT*G/2.0/C)/(1.0+dT*G/2.0/C) );
+		SetVI(n,pos[0],pos[1],pos[2], (dT/C)/(1.0+dT*G/2.0/C) );
+	}
+	else
+	{
+		SetVV(n,pos[0],pos[1],pos[2], 0 );
+		SetVI(n,pos[0],pos[1],pos[2], 0 );
+	}
+
+	double L = EC_L[n][i];
+	double R = EC_R[n][i];
+	if (L>0)
+	{
+		SetII(n,pos[0],pos[1],pos[2], (1.0-dT*R/2.0/L)/(1.0+dT*R/2.0/L) );
+		SetIV(n,pos[0],pos[1],pos[2], (dT/L)/(1.0+dT*R/2.0/L) );
+	}
+	else
+	{
+		SetII(n,pos[0],pos[1],pos[2], 0 );
+		SetIV(n,pos[0],pos[1],pos[2], 0 );
+	}
+}
+
+int Operator::CalcECOperator( DebugFlags debugFlags )
+{
+	Init_EC();
+	InitDataStorage();
+
+	if (Calc_EC()==0)
+		return -1;
+
+	m_InvaildTimestep = false;
+	opt_dT = 0;
+	if (dT>0)
+	{
+		double save_dT = dT;
+		CalcTimestep();
+		opt_dT = dT;
+		if (dT<save_dT)
+		{
+			cerr << "Operator::CalcECOperator: Warning, forced timestep: " << save_dT << "s is larger than calculated timestep: " << dT << "s! It is not recommended using this timestep!! " << endl;
+			m_InvaildTimestep = true;
+		}
+
+		dT = save_dT;
+	}
+	else
+		CalcTimestep();
+
+	dT*=m_TimeStepFactor;
+
+	if (m_Exc->GetSignalPeriod()>0)
+	{
+		unsigned int TS = ceil(m_Exc->GetSignalPeriod()/dT);
+		double new_dT = m_Exc->GetSignalPeriod()/TS;
+		cout << "Operartor::CalcECOperator: Decreasing timestep by " << round((dT-new_dT)/dT*1000)/10.0 << "% to " << new_dT << " (" << dT << ") to match periodic signal" << endl;
+		dT = new_dT;
+	}
+
+	m_Exc->Reset(dT);
+
+	InitOperator();
+
+	unsigned int pos[3];
+
+	for (int n=0; n<3; ++n)
+	{
+		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])
+				{
+					Calc_ECOperatorPos(n,pos);
+				}
+			}
+		}
+	}
+
+	//Apply PEC to all boundary's
+	bool PEC[6]={1,1,1,1,1,1};
+	//make an exception for BC == -1
+	for (int n=0; n<6; ++n)
+		if ((m_BC[n]==-1))
+			PEC[n] = false;
+	ApplyElectricBC(PEC);
+
+	CalcPEC();
+
+	Calc_LumpedElements();
+
+	bool PMC[6];
+	for (int n=0; n<6; ++n)
+		PMC[n] = m_BC[n]==1;
+	ApplyMagneticBC(PMC);
+
+	//all information available for extension... create now...
+	for (size_t n=0; n<m_Op_exts.size(); ++n)
+		m_Op_exts.at(n)->BuildExtension();
+
+	//remove inactive extensions
+	vector<Operator_Extension*>::iterator it = m_Op_exts.begin();
+	while (it!=m_Op_exts.end())
+	{
+		if ( (*it)->IsActive() == false)
+		{
+			DeleteExtension((*it));
+			it = m_Op_exts.begin(); //restart search for inactive extension
+		}
+		else
+			++it;
+	}
+
+	if (debugFlags & debugMaterial)
+		DumpMaterial2File( "material_dump" );
+	if (debugFlags & debugOperator)
+		DumpOperator2File( "operator_dump" );
+	if (debugFlags & debugPEC)
+		DumpPEC2File( "PEC_dump" );
+
+	//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;
+	}
+
+	return 0;
+}
+
+void Operator::ApplyElectricBC(bool* dirs)
+{
+	if (!dirs)
+		return;
+
+	unsigned int pos[3];
+	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])
+			{
+				if (dirs[2*n])
+				{
+					// set to PEC
+					pos[n] = 0;
+					SetVV(nP, pos[0],pos[1],pos[2], 0 );
+					SetVI(nP, pos[0],pos[1],pos[2], 0 );
+					SetVV(nPP,pos[0],pos[1],pos[2], 0 );
+					SetVI(nPP,pos[0],pos[1],pos[2], 0 );
+				}
+
+				if (dirs[2*n+1])
+				{
+					// set to PEC
+					pos[n] = numLines[n]-1;
+					SetVV(n,  pos[0],pos[1],pos[2], 0 ); // these are outside the FDTD-domain as defined by the main disc
+					SetVI(n,  pos[0],pos[1],pos[2], 0 ); // these are outside the FDTD-domain as defined by the main disc
+
+					SetVV(nP, pos[0],pos[1],pos[2], 0 );
+					SetVI(nP, pos[0],pos[1],pos[2], 0 );
+					SetVV(nPP,pos[0],pos[1],pos[2], 0 );
+					SetVI(nPP,pos[0],pos[1],pos[2], 0 );
+				}
+			}
+		}
+	}
+}
+
+void Operator::ApplyMagneticBC(bool* dirs)
+{
+	if (!dirs)
+		return;
+
+	unsigned int pos[3];
+	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])
+			{
+				if (dirs[2*n])
+				{
+					// set to PMC
+					pos[n] = 0;
+					SetII(n,  pos[0],pos[1],pos[2], 0 );
+					SetIV(n,  pos[0],pos[1],pos[2], 0 );
+					SetII(nP, pos[0],pos[1],pos[2], 0 );
+					SetIV(nP, pos[0],pos[1],pos[2], 0 );
+					SetII(nPP,pos[0],pos[1],pos[2], 0 );
+					SetIV(nPP,pos[0],pos[1],pos[2], 0 );
+				}
+
+				if (dirs[2*n+1])
+				{
+					// set to PMC
+					pos[n] = numLines[n]-2;
+					SetII(nP, pos[0],pos[1],pos[2], 0 );
+					SetIV(nP, pos[0],pos[1],pos[2], 0 );
+					SetII(nPP,pos[0],pos[1],pos[2], 0 );
+					SetIV(nPP,pos[0],pos[1],pos[2], 0 );
+				}
+
+				// the last current lines are outside the FDTD domain and cannot be iterated by the FDTD engine
+				pos[n] = numLines[n]-1;
+				SetII(n,  pos[0],pos[1],pos[2], 0 );
+				SetIV(n,  pos[0],pos[1],pos[2], 0 );
+				SetII(nP, pos[0],pos[1],pos[2], 0 );
+				SetIV(nP, pos[0],pos[1],pos[2], 0 );
+				SetII(nPP,pos[0],pos[1],pos[2], 0 );
+				SetIV(nPP,pos[0],pos[1],pos[2], 0 );
+			}
+		}
+	}
+}
+
+bool Operator::Calc_ECPos(int ny, const unsigned int* pos, double* EC, vector<CSPrimitives*> vPrims) const
+{
+	double EffMat[4];
+	Calc_EffMatPos(ny,pos,EffMat, vPrims);
+
+	if (m_epsR)
+		m_epsR[ny][pos[0]][pos[1]][pos[2]] =  EffMat[0];
+	if (m_kappa)
+		m_kappa[ny][pos[0]][pos[1]][pos[2]] =  EffMat[1];
+	if (m_mueR)
+		m_mueR[ny][pos[0]][pos[1]][pos[2]] =  EffMat[2];
+	if (m_sigma)
+		m_sigma[ny][pos[0]][pos[1]][pos[2]] =  EffMat[3];
+
+	double delta = GetEdgeLength(ny,pos);
+	double area  = GetEdgeArea(ny,pos);
+
+//	if (isnan(EffMat[0]))
+//	{
+//		cerr << ny << " " << pos[0] << " " << pos[1] << " " << pos[2] << " : " << EffMat[0] << endl;
+//	}
+
+	if (delta)
+	{
+		EC[0] = EffMat[0] * area/delta;
+		EC[1] = EffMat[1] * area/delta;
+	}
+	else
+	{
+		EC[0] = 0;
+		EC[1] = 0;
+	}
+
+	delta = GetEdgeLength(ny,pos,true);
+	area  = GetEdgeArea(ny,pos,true);
+
+	if (delta)
+	{
+		EC[2] = EffMat[2] * area/delta;
+		EC[3] = EffMat[3] * area/delta;
+	}
+	else
+	{
+		EC[2] = 0;
+		EC[3] = 0;
+	}
+
+	return true;
+}
+
+double Operator::GetRawDiscDelta(int ny, const int pos) const
+{
+	//numLines[ny] is expected to be larger then 1 !
+
+	if (pos<0)
+		return (discLines[ny][0] - discLines[ny][1]);
+	if (pos>=(int)numLines[ny]-1)
+		return (discLines[ny][numLines[ny]-2] - discLines[ny][numLines[ny]-1]);
+
+	return (discLines[ny][pos+1] - discLines[ny][pos]);
+}
+
+bool Operator::GetCellCenterMaterialAvgCoord(const int pos[], double coord[3]) const
+{
+	unsigned int ui_pos[3];
+	for (int n=0;n<3;++n)
+	{
+		if ((pos[n]<0) || (pos[n]>=(int)numLines[n]))
+			return false;
+		ui_pos[n] = pos[n];
+	}
+	GetNodeCoords(ui_pos, coord, true);
+	return true;
+}
+
+double Operator::GetMaterial(int ny, const double* coords, int MatType, vector<CSPrimitives*> vPrims, bool markAsUsed) const
+{
+	CSProperties* prop = CSX->GetPropertyByCoordPriority(coords,vPrims,markAsUsed);
+//	CSProperties* old_prop = CSX->GetPropertyByCoordPriority(coords,CSProperties::MATERIAL,markAsUsed);
+//	if (old_prop!=prop)
+//	{
+//		cerr << "ERROR: Unequal properties!" << endl;
+//		exit(-1);
+//	}
+
+	CSPropMaterial* mat = dynamic_cast<CSPropMaterial*>(prop);
+	if (mat)
+	{
+		switch (MatType)
+		{
+		case 0:
+			return mat->GetEpsilonWeighted(ny,coords);
+		case 1:
+			return mat->GetKappaWeighted(ny,coords);
+		case 2:
+			return mat->GetMueWeighted(ny,coords);
+		case 3:
+			return mat->GetSigmaWeighted(ny,coords);
+		case 4:
+			return mat->GetDensityWeighted(coords);
+		default:
+			cerr << "Operator::GetMaterial: Error: unknown material type" << endl;
+			return 0;
+		}
+	}
+
+	switch (MatType)
+	{
+	case 0:
+		return GetBackgroundEpsR();
+	case 1:
+		return GetBackgroundKappa();
+	case 2:
+		return GetBackgroundMueR();
+	case 3:
+		return GetBackgroundSigma();
+	case 4:
+		return GetBackgroundDensity();
+	default:
+		cerr << "Operator::GetMaterial: Error: unknown material type" << endl;
+		return 0;
+	}
+}
+
+bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* EffMat, vector<CSPrimitives *> vPrims) const
+{
+	int n=ny;
+	double coord[3];
+	int nP = (n+1)%3;
+	int nPP = (n+2)%3;
+
+	int loc_pos[3] = {(int)pos[0],(int)pos[1],(int)pos[2]};
+	double A_n;
+	double area = 0;
+	EffMat[0] = 0;
+	EffMat[1] = 0;
+	EffMat[2] = 0;
+	EffMat[3] = 0;
+
+	//******************************* epsilon,kappa averaging *****************************//
+	//shift up-right
+	if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
+	{
+		A_n = GetNodeArea(ny,loc_pos,true);
+		EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n;
+		EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n;
+		area+=A_n;
+	}
+
+	//shift up-left
+	--loc_pos[nP];
+	if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
+	{
+		A_n = GetNodeArea(ny,loc_pos,true);
+		EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n;
+		EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n;
+		area+=A_n;
+	}
+
+	//shift down-right
+	++loc_pos[nP];
+	--loc_pos[nPP];
+	if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
+	{
+		A_n = GetNodeArea(ny,loc_pos,true);
+		EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n;
+		EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n;
+		area+=A_n;
+	}
+
+	//shift down-left
+	--loc_pos[nP];
+	if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
+	{
+		A_n = GetNodeArea(ny,loc_pos,true);
+		EffMat[0] += GetMaterial(n, coord, 0, vPrims)*A_n;
+		EffMat[1] += GetMaterial(n, coord, 1, vPrims)*A_n;
+		area+=A_n;
+	}
+
+	EffMat[0]*=__EPS0__/area;
+	EffMat[1]/=area;
+
+	//******************************* mu,sigma averaging *****************************//
+	loc_pos[0]=pos[0];
+	loc_pos[1]=pos[1];
+	loc_pos[2]=pos[2];
+	double length=0;
+	double delta_ny,sigma;
+	//shift down
+	--loc_pos[n];
+	if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
+	{
+		delta_ny = GetNodeWidth(n,loc_pos,true);
+		EffMat[2] += delta_ny / GetMaterial(n, coord, 2, vPrims);
+		sigma = GetMaterial(n, coord, 3, vPrims);
+		if (sigma)
+			EffMat[3] += delta_ny / sigma;
+		else
+			EffMat[3] = 0;
+		length+=delta_ny;
+	}
+
+	//shift up
+	++loc_pos[n];
+	if (GetCellCenterMaterialAvgCoord(loc_pos,coord))
+	{
+		delta_ny = GetNodeWidth(n,loc_pos,true);
+		EffMat[2] += delta_ny / GetMaterial(n, coord, 2, vPrims);
+		sigma = GetMaterial(n, coord, 3, vPrims);
+		if (sigma)
+			EffMat[3] += delta_ny / sigma;
+		else
+			EffMat[3] = 0;
+		length+=delta_ny;
+	}
+
+	EffMat[2] = length * __MUE0__ / EffMat[2];
+	if (EffMat[3]) EffMat[3]=length / EffMat[3];
+
+	for (int n=0; n<4; ++n)
+		if (isnan(EffMat[n]) || isinf(EffMat[n]))
+		{
+			cerr << "Operator::" << __func__ << ": Error, an effective material parameter is not a valid result, this should NOT have happend... exit..." << endl;
+			cerr << ny << "@" << n << " : " << pos[0] << "," << pos[1] << ","  << pos[2] << endl;
+			exit(0);
+		}
+	return true;
+}
+
+bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* EffMat, vector<CSPrimitives*> vPrims) const
+{
+	int n=ny;
+	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]];
+	double delta=GetRawDiscDelta(n,pos[n]);
+	double deltaP=GetRawDiscDelta(nP,pos[nP]);
+	double deltaPP=GetRawDiscDelta(nPP,pos[nPP]);
+	double delta_M=GetRawDiscDelta(n,pos[n]-1);
+	double deltaP_M=GetRawDiscDelta(nP,pos[nP]-1);
+	double deltaPP_M=GetRawDiscDelta(nPP,pos[nPP]-1);
+
+	int loc_pos[3] = {(int)pos[0],(int)pos[1],(int)pos[2]};
+	double A_n;
+	double area = 0;
+
+	//******************************* 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;
+	A_n = GetNodeArea(ny,loc_pos,true);
+	EffMat[0] = GetMaterial(n, shiftCoord, 0, vPrims)*A_n;
+	EffMat[1] = GetMaterial(n, shiftCoord, 1, vPrims)*A_n;
+	area+=A_n;
+
+	//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;
+
+	--loc_pos[nP];
+	A_n = GetNodeArea(ny,loc_pos,true);
+	EffMat[0] += GetMaterial(n, shiftCoord, 0, vPrims)*A_n;
+	EffMat[1] += GetMaterial(n, shiftCoord, 1, vPrims)*A_n;
+	area+=A_n;
+
+	//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;
+	++loc_pos[nP];
+	--loc_pos[nPP];
+	A_n = GetNodeArea(ny,loc_pos,true);
+	EffMat[0] += GetMaterial(n, shiftCoord, 0, vPrims)*A_n;
+	EffMat[1] += GetMaterial(n, shiftCoord, 1, vPrims)*A_n;
+	area+=A_n;
+
+	//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;
+	--loc_pos[nP];
+	A_n = GetNodeArea(ny,loc_pos,true);
+	EffMat[0] += GetMaterial(n, shiftCoord, 0, vPrims)*A_n;
+	EffMat[1] += GetMaterial(n, shiftCoord, 1, vPrims)*A_n;
+	area+=A_n;
+
+	EffMat[0]*=__EPS0__/area;
+	EffMat[1]/=area;
+
+	//******************************* mu,sigma averaging *****************************//
+	loc_pos[0]=pos[0];
+	loc_pos[1]=pos[1];
+	loc_pos[2]=pos[2];
+	double length=0;
+
+	//shift down
+	shiftCoord[n] = coord[n]-delta_M*0.25;
+	shiftCoord[nP] = coord[nP]+deltaP*0.5;
+	shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
+	--loc_pos[n];
+	double delta_ny = GetNodeWidth(n,loc_pos,true);
+	EffMat[2] = delta_ny / GetMaterial(n, shiftCoord, 2, vPrims);
+	double sigma = GetMaterial(n, shiftCoord, 3, vPrims);
+	if (sigma)
+		EffMat[3] = delta_ny / sigma;
+	else
+		EffMat[3] = 0;
+	length=delta_ny;
+
+	//shift up
+	shiftCoord[n] = coord[n]+delta*0.25;
+	shiftCoord[nP] = coord[nP]+deltaP*0.5;
+	shiftCoord[nPP] = coord[nPP]+deltaPP*0.5;
+	++loc_pos[n];
+	delta_ny = GetNodeWidth(n,loc_pos,true);
+	EffMat[2] += delta_ny / GetMaterial(n, shiftCoord, 2, vPrims);
+	sigma = GetMaterial(n, shiftCoord, 3, vPrims);
+	if (sigma)
+		EffMat[3] += delta_ny / sigma;
+	else
+		EffMat[3] = 0;
+	length+=delta_ny;
+
+	EffMat[2] = length * __MUE0__ / EffMat[2];
+	if (EffMat[3]) EffMat[3]=length / EffMat[3];
+
+	for (int n=0; n<4; ++n)
+		if (isnan(EffMat[n]) || isinf(EffMat[n]))
+		{
+			cerr << "Operator::" << __func__ << ": Error, An effective material parameter is not a valid result, this should NOT have happend... exit..." << endl;
+			cerr << ny << "@" << n << " : " << pos[0] << "," << pos[1] << ","  << pos[2] << endl;
+			exit(0);
+		}
+
+	return true;
+}
+
+bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat, vector<CSPrimitives *> vPrims) const
+{
+	switch (m_MatAverageMethod)
+	{
+	case QuarterCell:
+		return AverageMatQuarterCell(ny, pos, EffMat, vPrims);
+	case CentralCell:
+		return AverageMatCellCenter(ny, pos, EffMat, vPrims);
+	default:
+		cerr << "Operator:: " << __func__ << ":  Error, unknown material averaging method... exit" << endl;
+		exit(1);
+	}
+	return false;
+}
+
+bool Operator::Calc_LumpedElements()
+{
+	vector<CSProperties*> props = CSX->GetPropertyByType(CSProperties::LUMPED_ELEMENT);
+	for (size_t i=0;i<props.size();++i)
+	{
+		CSPropLumpedElement* PLE = dynamic_cast<CSPropLumpedElement*>(props.at(i));
+		if (PLE==NULL)
+			return false; //sanity check: this should never happen!
+		vector<CSPrimitives*> prims = PLE->GetAllPrimitives();
+		for (size_t bn=0;bn<prims.size();++bn)
+		{
+			CSPrimBox* box = dynamic_cast<CSPrimBox*>(prims.at(bn));
+			if (box)
+			{	//calculate lumped element parameter
+
+				double C = PLE->GetCapacity();
+				if (C<=0)
+					C = NAN;
+				double R = PLE->GetResistance();
+				if (R<0)
+					R = NAN;
+
+				if ((isnan(R)) && (isnan(C)))
+				{
+					cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element R or C not specified! skipping. "
+							<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
+					continue;
+				}
+
+				int ny = PLE->GetDirection();
+				if ((ny<0) || (ny>2))
+				{
+					cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element direction is invalid! skipping. "
+							<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
+					continue;
+				}
+				int nyP = (ny+1)%3;
+				int nyPP = (ny+2)%3;
+
+				unsigned int uiStart[3];
+				unsigned int uiStop[3];
+				// snap to the native coordinate system
+				int Snap_Dimension = Operator::SnapBox2Mesh(box->GetStartCoord()->GetCoords(m_MeshType), box->GetStopCoord()->GetCoords(m_MeshType), uiStart, uiStop, false, true);
+				if (Snap_Dimension<=0)
+				{
+					if (Snap_Dimension>=-1)
+						cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element snapping failed! Dimension is: " << Snap_Dimension << " skipping. "
+								<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
+					// Snap_Dimension == -2 means outside the simulation domain --> no special warning, but box probably marked as unused!
+					continue;
+				}
+
+				if (uiStart[ny]==uiStop[ny])
+				{
+					cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element with zero (snapped) length is invalid! skipping. "
+							<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
+					continue;
+				}
+
+				//calculate geometric property for this lumped element
+				unsigned int pos[3];
+				double unitGC=0;
+				int ipos=0;
+				for (pos[ny]=uiStart[ny];pos[ny]<uiStop[ny];++pos[ny])
+				{
+					double unitGC_Plane=0;
+					for (pos[nyP]=uiStart[nyP];pos[nyP]<=uiStop[nyP];++pos[nyP])
+					{
+						for (pos[nyPP]=uiStart[nyPP];pos[nyPP]<=uiStop[nyPP];++pos[nyPP])
+						{
+							// capacity/conductivity in parallel: add values
+							unitGC_Plane += GetEdgeArea(ny,pos)/GetEdgeLength(ny,pos);
+						}
+					}
+
+					//capacity/conductivity in series: add reciprocal values
+					unitGC += 1/unitGC_Plane;
+				}
+				unitGC = 1/unitGC;
+
+				bool caps = PLE->GetCaps();
+				double kappa = 0;
+				double epsilon = 0;
+				if (R>0)
+					kappa = 1 / R / unitGC;
+				if (C>0)
+				{
+					epsilon =  C / unitGC;
+
+					if (epsilon< __EPS0__)
+					{
+						cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element capacity is too small for its size! skipping. "
+								<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
+						C = 0;
+					}
+				}
+
+				for (pos[ny]=uiStart[ny];pos[ny]<uiStop[ny];++pos[ny])
+				{
+					for (pos[nyP]=uiStart[nyP];pos[nyP]<=uiStop[nyP];++pos[nyP])
+					{
+						for (pos[nyPP]=uiStart[nyPP];pos[nyPP]<=uiStop[nyPP];++pos[nyPP])
+						{
+							ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
+							if (C>0)
+								EC_C[ny][ipos] = epsilon * GetEdgeArea(ny,pos)/GetEdgeLength(ny,pos);
+							if (R>0)
+								EC_G[ny][ipos] = kappa * GetEdgeArea(ny,pos)/GetEdgeLength(ny,pos);
+
+							if (R==0) //make lumped element a PEC if resistance is zero
+							{
+								SetVV(ny,pos[0],pos[1],pos[2], 0 );
+								SetVI(ny,pos[0],pos[1],pos[2], 0 );
+							}
+							else //recalculate operator inside the lumped element
+								Calc_ECOperatorPos(ny,pos);
+						}
+					}
+				}
+
+				// setup metal caps
+				if (caps)
+				{
+					for (pos[nyP]=uiStart[nyP];pos[nyP]<=uiStop[nyP];++pos[nyP])
+					{
+						for (pos[nyPP]=uiStart[nyPP];pos[nyPP]<=uiStop[nyPP];++pos[nyPP])
+						{
+							pos[ny]=uiStart[ny];
+							if (pos[nyP]<uiStop[nyP])
+							{
+								SetVV(nyP,pos[0],pos[1],pos[2], 0 );
+								SetVI(nyP,pos[0],pos[1],pos[2], 0 );
+								++m_Nr_PEC[nyP];
+							}
+
+							if (pos[nyPP]<uiStop[nyPP])
+							{
+								SetVV(nyPP,pos[0],pos[1],pos[2], 0 );
+								SetVI(nyPP,pos[0],pos[1],pos[2], 0 );
+								++m_Nr_PEC[nyPP];
+							}
+
+							pos[ny]=uiStop[ny];
+							if (pos[nyP]<uiStop[nyP])
+							{
+								SetVV(nyP,pos[0],pos[1],pos[2], 0 );
+								SetVI(nyP,pos[0],pos[1],pos[2], 0 );
+								++m_Nr_PEC[nyP];
+							}
+
+							if (pos[nyPP]<uiStop[nyPP])
+							{
+								SetVV(nyPP,pos[0],pos[1],pos[2], 0 );
+								SetVI(nyPP,pos[0],pos[1],pos[2], 0 );
+								++m_Nr_PEC[nyPP];
+							}
+						}
+					}
+				}
+				box->SetPrimitiveUsed(true);
+
+			}
+			else
+				cerr << "Operator::Calc_LumpedElements(): Warning: Primitves other than boxes are not supported for lumped elements! skipping "
+						<< prims.at(bn)->GetTypeName() << " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
+		}
+	}
+	return true;
+}
+
+void Operator::Init_EC()
+{
+	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 FDTD_FLOAT[MainOp->GetSize()];
+		EC_G[n] = new FDTD_FLOAT[MainOp->GetSize()];
+		EC_L[n] = new FDTD_FLOAT[MainOp->GetSize()];
+		EC_R[n] = new FDTD_FLOAT[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;
+		}
+	}
+}
+
+bool Operator::Calc_EC()
+{
+	if (CSX==NULL)
+	{
+		cerr << "CartOperator::Calc_EC: CSX not given or invalid!!!" << endl;
+		return false;
+	}
+	
+	MainOp->SetPos(0,0,0);
+	Calc_EC_Range(0,numLines[0]-1);
+	return true;
+}
+
+vector<CSPrimitives*> Operator::GetPrimitivesBoundBox(int posX, int posY, int posZ, CSProperties::PropertyType type) const
+{
+	double boundBox[6];
+	int BBpos[3] = {posX, posY, posZ};
+	for (int n=0;n<3;++n)
+	{
+		if (BBpos[n]<0)
+		{
+			boundBox[2*n]   = this->GetDiscLine(n,0);
+			boundBox[2*n+1] = this->GetDiscLine(n,numLines[n]-1);
+		}
+		else
+		{
+			boundBox[2*n]   = this->GetDiscLine(n, max(0, BBpos[n]-1));
+			boundBox[2*n+1] = this->GetDiscLine(n, min(int(numLines[n])-1, BBpos[n]+1));
+		}
+	}
+
+	vector<CSPrimitives*> vPrim = this->CSX->GetPrimitivesByBoundBox(boundBox, true, type);
+	return vPrim;
+}
+
+void Operator::Calc_EC_Range(unsigned int xStart, unsigned int xStop)
+{
+//	vector<CSPrimitives*> vPrims = this->CSX->GetAllPrimitives(true, CSProperties::MATERIAL);
+	unsigned int ipos;
+	unsigned int pos[3];
+	double inEC[4];
+	for (pos[0]=xStart; pos[0]<=xStop; ++pos[0])
+	{
+		for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
+		{
+			vector<CSPrimitives*> vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, CSProperties::MATERIAL);
+			for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
+			{
+				ipos = MainOp->GetPos(pos[0],pos[1],pos[2]);
+				for (int n=0; n<3; ++n)
+				{
+					Calc_ECPos(n,pos,inEC,vPrims);
+					EC_C[n][ipos]=inEC[0];
+					EC_G[n][ipos]=inEC[1];
+					EC_L[n][ipos]=inEC[2];
+					EC_R[n][ipos]=inEC[3];
+				}
+			}
+		}
+	}
+}
+
+void Operator::SetTimestepFactor(double factor)
+{
+	if ((factor<=0) || (factor>1))
+	{
+		cerr << "Operator::SetTimestepFactor: Warning, invalid timestep factor, skipping!" << endl;
+		return;
+	}
+
+	cout << "Operator::SetTimestepFactor: Setting timestep factor to " << factor << endl;
+	m_TimeStepFactor=factor;
+}
+
+double Operator::CalcTimestep()
+{
+	if (m_TimeStepVar==3)
+		return CalcTimestep_Var3(); //the biggest one for cartesian meshes
+
+	//variant 1 is default
+	return CalcTimestep_Var1();
+}
+
+////Berechnung nach Andreas Rennings Dissertation 2008, Seite 66, Formel 4.52
+double Operator::CalcTimestep_Var1()
+{
+	m_Used_TS_Name = string("Rennings_1");
+//	cout << "Operator::CalcTimestep(): Using timestep algorithm by Andreas Rennings, Dissertation @ University Duisburg-Essen, 2008, pp. 66, eq. 4.52" << endl;
+	dT=1e200;
+	double newT;
+	unsigned int pos[3];
+	unsigned int smallest_pos[3] = {0, 0, 0};
+	unsigned int smallest_n = 0;
+	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) && (newT>0.0))
+					{
+						dT=newT;
+						smallest_pos[0]=pos[0];smallest_pos[1]=pos[1];smallest_pos[2]=pos[2];
+						smallest_n = n;
+					}
+				}
+			}
+		}
+	}
+	if (dT==0)
+	{
+		cerr << "Operator::CalcTimestep: Timestep is zero... this is not supposed to happen!!! exit!" << endl;
+		exit(3);
+	}
+	if (g_settings.GetVerboseLevel()>1)
+	{
+		cout << "Operator::CalcTimestep_Var1: Smallest timestep (" << dT << "s) found at position: " <<  smallest_n << " : " << smallest_pos[0] << ";" <<  smallest_pos[1] << ";" <<  smallest_pos[2] << endl;
+	}
+	return 0;
+}
+
+double min(double* val, unsigned int count)
+{
+	if (count==0)
+		return 0.0;
+	double min = val[0];
+	for (unsigned int n=1; n<count; ++n)
+		if (val[n]<min)
+			min = val[n];
+	return min;
+}
+
+//Berechnung nach Andreas Rennings Dissertation 2008, Seite 76 ff, Formel 4.77 ff
+double Operator::CalcTimestep_Var3()
+{
+	dT=1e200;
+	m_Used_TS_Name = string("Rennings_2");
+//	cout << "Operator::CalcTimestep(): Using timestep algorithm by Andreas Rennings, Dissertation @ University Duisburg-Essen, 2008, pp. 76, eq. 4.77 ff." << endl;
+	double newT;
+	unsigned int pos[3];
+	unsigned int smallest_pos[3] = {0, 0, 0};
+	unsigned int smallest_n = 0;
+	unsigned int ipos;
+	double w_total=0;
+	double wqp=0,wt1=0,wt2=0;
+	double wt_4[4]={0,0,0,0};
+	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])
+				{
+					MainOp->ResetShift();
+					ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
+					wqp  = 1/(EC_L[nPP][ipos]*EC_C[n][MainOp->GetShiftedPos(nP ,1)]) + 1/(EC_L[nPP][ipos]*EC_C[n][ipos]);
+					wqp += 1/(EC_L[nP ][ipos]*EC_C[n][MainOp->GetShiftedPos(nPP,1)]) + 1/(EC_L[nP ][ipos]*EC_C[n][ipos]);
+					ipos = MainOp->Shift(nP,-1);
+					wqp += 1/(EC_L[nPP][ipos]*EC_C[n][MainOp->GetShiftedPos(nP ,1)]) + 1/(EC_L[nPP][ipos]*EC_C[n][ipos]);
+					ipos = MainOp->Shift(nPP,-1);
+					wqp += 1/(EC_L[nP ][ipos]*EC_C[n][MainOp->GetShiftedPos(nPP,1)]) + 1/(EC_L[nP ][ipos]*EC_C[n][ipos]);
+
+					MainOp->ResetShift();
+					ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
+					wt_4[0] = 1/(EC_L[nPP][ipos]						  *EC_C[nP ][ipos]);
+					wt_4[1] = 1/(EC_L[nPP][MainOp->GetShiftedPos(nP ,-1)] *EC_C[nP ][ipos]);
+					wt_4[2] = 1/(EC_L[nP ][ipos]						  *EC_C[nPP][ipos]);
+					wt_4[3] = 1/(EC_L[nP ][MainOp->GetShiftedPos(nPP,-1)] *EC_C[nPP][ipos]);
+
+					wt1 = wt_4[0]+wt_4[1]+wt_4[2]+wt_4[3] - 2*min(wt_4,4);
+
+					MainOp->ResetShift();
+					ipos = MainOp->SetPos(pos[0],pos[1],pos[2]);
+					wt_4[0] = 1/(EC_L[nPP][ipos]						  *EC_C[nP ][MainOp->GetShiftedPos(n,1)]);
+					wt_4[1] = 1/(EC_L[nPP][MainOp->GetShiftedPos(nP ,-1)] *EC_C[nP ][MainOp->GetShiftedPos(n,1)]);
+					wt_4[2] = 1/(EC_L[nP ][ipos]						  *EC_C[nPP][MainOp->GetShiftedPos(n,1)]);
+					wt_4[3] = 1/(EC_L[nP ][MainOp->GetShiftedPos(nPP,-1)] *EC_C[nPP][MainOp->GetShiftedPos(n,1)]);
+
+					wt2 = wt_4[0]+wt_4[1]+wt_4[2]+wt_4[3] - 2*min(wt_4,4);
+
+					w_total = wqp + wt1 + wt2;
+					newT = 2/sqrt( w_total );
+					if ((newT<dT) && (newT>0.0))
+					{
+						dT=newT;
+						smallest_pos[0]=pos[0];smallest_pos[1]=pos[1];smallest_pos[2]=pos[2];
+						smallest_n = n;
+					}
+				}
+			}
+		}
+	}
+	if (dT==0)
+	{
+		cerr << "Operator::CalcTimestep: Timestep is zero... this is not supposed to happen!!! exit!" << endl;
+		exit(3);
+	}
+	if (g_settings.GetVerboseLevel()>1)
+	{
+		cout << "Operator::CalcTimestep_Var3: Smallest timestep (" << dT << "s) found at position: " <<  smallest_n << " : " << smallest_pos[0] << ";" <<  smallest_pos[1] << ";" <<  smallest_pos[2] << endl;
+	}
+	return 0;
+}
+
+bool Operator::CalcPEC()
+{
+	m_Nr_PEC[0]=0;
+	m_Nr_PEC[1]=0;
+	m_Nr_PEC[2]=0;
+
+	CalcPEC_Range(0,numLines[0]-1,m_Nr_PEC);
+
+	CalcPEC_Curves();
+
+	return true;
+}
+
+void Operator::CalcPEC_Range(unsigned int startX, unsigned int stopX, unsigned int* counter)
+{
+	double coord[3];
+	unsigned int pos[3];
+	for (pos[0]=startX; pos[0]<=stopX; ++pos[0])
+	{
+		for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
+		{
+			vector<CSPrimitives*> vPrims = this->GetPrimitivesBoundBox(pos[0], pos[1], -1, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL));
+			for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
+			{
+				for (int n=0; n<3; ++n)
+				{
+					GetYeeCoords(n,pos,coord,false);
+					CSProperties* prop = CSX->GetPropertyByCoordPriority(coord, vPrims, true);
+//					CSProperties* old_prop = CSX->GetPropertyByCoordPriority(coord, (CSProperties::PropertyType)(CSProperties::MATERIAL | CSProperties::METAL), true);
+//					if (old_prop!=prop)
+//					{
+//						cerr << "CalcPEC_Range: " << old_prop << " vs " << prop << endl;
+//						exit(-1);
+//					}
+					if (prop)
+					{
+						if (prop->GetType()==CSProperties::METAL) //set to PEC
+						{
+							SetVV(n,pos[0],pos[1],pos[2], 0 );
+							SetVI(n,pos[0],pos[1],pos[2], 0 );
+							++counter[n];
+						}
+					}
+				}
+			}
+		}
+	}
+}
+
+void Operator::CalcPEC_Curves()
+{
+	//special treatment for primitives of type curve (treated as wires)
+	double p1[3];
+	double p2[3];
+	Grid_Path path;
+	vector<CSProperties*> vec_prop = CSX->GetPropertyByType(CSProperties::METAL);
+	for (size_t p=0; p<vec_prop.size(); ++p)
+	{
+		CSProperties* prop = vec_prop.at(p);
+		for (size_t n=0; n<prop->GetQtyPrimitives(); ++n)
+		{
+			CSPrimitives* prim = prop->GetPrimitive(n);
+			CSPrimCurve* curv = prim->ToCurve();
+			if (curv)
+			{
+				for (size_t i=1; i<curv->GetNumberOfPoints(); ++i)
+				{
+					curv->GetPoint(i-1,p1,m_MeshType);
+					curv->GetPoint(i,p2,m_MeshType);
+					path = FindPath(p1,p2);
+					if (path.dir.size()>0)
+						prim->SetPrimitiveUsed(true);
+					for (size_t t=0; t<path.dir.size(); ++t)
+					{
+						SetVV(path.dir.at(t),path.posPath[0].at(t),path.posPath[1].at(t),path.posPath[2].at(t), 0 );
+						SetVI(path.dir.at(t),path.posPath[0].at(t),path.posPath[1].at(t),path.posPath[2].at(t), 0 );
+						++m_Nr_PEC[path.dir.at(t)];
+					}
+				}
+			}
+		}
+	}
+}
+
+Operator_Ext_Excitation* Operator::GetExcitationExtension() const
+{
+	//search for excitation extension
+	Operator_Ext_Excitation* Op_Ext_Exc=0;
+	for (size_t n=0; n<m_Op_exts.size(); ++n)
+	{
+		Op_Ext_Exc = dynamic_cast<Operator_Ext_Excitation*>(m_Op_exts.at(n));
+		if (Op_Ext_Exc)
+			break;
+	}
+	return Op_Ext_Exc;
+}
+
+void Operator::AddExtension(Operator_Extension* op_ext)
+{
+	m_Op_exts.push_back(op_ext);
+}
+
+void Operator::DeleteExtension(Operator_Extension* op_ext)
+{
+	for (size_t n=0;n<m_Op_exts.size();++n)
+	{
+		if (m_Op_exts.at(n)==op_ext)
+		{
+			m_Op_exts.erase(m_Op_exts.begin()+n);
+			return;
+		}
+	}
+}
+
+double Operator::CalcNumericPhaseVelocity(unsigned int start[3], unsigned int stop[3], double propDir[3], float freq) const
+{
+	double average_mesh_disc[3];
+	double c0 = __C0__/sqrt(GetBackgroundEpsR()*GetBackgroundMueR());
+
+	//calculate average mesh deltas
+	for (int n=0;n<3;++n)
+	{
+		average_mesh_disc[n] = fabs(GetDiscLine(n,start[n])-GetDiscLine(n,stop[n]))*GetGridDelta() / (fabs(stop[n]-start[n]));
+	}
+
+	// if propagation is in a Cartesian direction, return analytic solution
+	for (int n=0;n<3;++n)
+	{
+		int nP = (n+1)%3;
+		int nPP = (n+2)%3;
+		if ((fabs(propDir[n])==1) && (propDir[nP]==0) && (propDir[nPP]==0))
+		{
+			double kx = asin(average_mesh_disc[0]/c0/dT*sin(2*PI*freq*dT/2))*2/average_mesh_disc[0];
+			return 2*PI*freq/kx;
+		}
+	}
+
+	// else, do an newton iterative estimation
+	double k0=2*PI*freq/c0;
+	double k=k0;
+	double RHS = pow(sin(2*PI*freq*dT/2)/c0/dT,2);
+	double fk=1,fdk=0;
+	double old_phv=0;
+	double phv=c0;
+	double err_est = 1e-6;
+	int it_count=0;
+	while (fabs(old_phv-phv)>err_est)
+	{
+		++it_count;
+		old_phv=phv;
+		fk=0;
+		fdk=0;
+		for (int n=0;n<3;++n)
+		{
+			fk+= pow(sin(propDir[n]*k*average_mesh_disc[n]/2)/average_mesh_disc[n],2);
+			fdk+= propDir[n]*sin(propDir[n]*k*average_mesh_disc[n]/2)*cos(propDir[n]*k*average_mesh_disc[n]/2)/average_mesh_disc[n];
+		}
+		fk -= RHS;
+		k-=fk/fdk;
+
+		// do not allow a speed greater than c0 due to a numerical inaccuracy
+		if (k<k0)
+			k=k0;
+
+		phv=2*PI*freq/k;
+
+		//abort if iteration count is getting to high!
+		if (it_count>99)
+		{
+			cerr << "Operator::CalcNumericPhaseVelocity: Error, newton iteration estimation can't find a solution!!" << endl;
+			break;
+		}
+	}
+
+	if (g_settings.GetVerboseLevel()>1)
+		cerr << "Operator::CalcNumericPhaseVelocity: Newton iteration estimated solution: " << phv/__C0__ << "*c0 in " <<  it_count << " iterations." << endl;
+
+	return phv;
+}