new excitation: Custom excite and some excitation related code-cleanup

FDTD/operator.cpp

@@ -19,6 +19,7 @@
 #include "operator.h"
 #include "processfields.h"
 #include "tools/array_ops.h"
+#include "fparser.hh"
 
 Operator::Operator()
 {
@@ -289,6 +290,35 @@ unsigned int Operator::CalcStepExcitation()
 	return 1;
 }
 
+unsigned int Operator::CalcCustomExcitation(double f0, int nTS, string signal)
+{
+	if (dT==0) return 0;
+	if (nTS<=0) return 0;
+
+	ExciteLength = (unsigned int)(nTS);
+//	cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
+	delete[] ExciteSignal;
+	ExciteSignal = new FDTD_FLOAT[ExciteLength+1];
+	ExciteSignal[0]=0.0;
+	FunctionParser fParse;
+	fParse.AddConstant("pi", 3.14159265358979323846);
+	fParse.AddConstant("e", 2.71828182845904523536);
+	fParse.Parse(signal,"t");
+	if (fParse.GetParseErrorType()!=FunctionParser::FP_NO_ERROR)
+	{
+		cerr << "Operator::CalcCustomExcitation: Function Parser error: " << fParse.ErrorMsg() << endl;
+		exit(1);
+	}
+	double vars[1];
+	for (unsigned int n=1;n<ExciteLength+1;++n)
+	{
+		vars[0] = (n-1)*GetTimestep();
+		ExciteSignal[n] = fParse.Eval(vars);
+//		cerr << ExciteSignal[n] << endl;
+	}
+	return CalcNyquistNum(f0);
+}
+
 unsigned int Operator::CalcSinusExcitation(double f0, int nTS)
 {
 	if (dT==0) return 0;

FDTD/operator.h

@@ -35,6 +35,8 @@ public:
 
 	virtual int CalcECOperator();
 
+	//! Calculate a custom signal \return number of Nyquist timesteps defined by f0
+	virtual unsigned int CalcCustomExcitation(double f0, int nTS, string signal);
 	//! Calculate an excitation with center of f0 and the half bandwidth fc \return number of Nyquist timesteps
 	virtual unsigned int CalcGaussianPulsExcitation(double f0, double fc);
 	//! Calculate a sinusoidal excitation with frequency f0 and a duration of nTS number of timesteps \return number of Nyquist timesteps

matlab/SetCustomExcite.m

@@ -0,0 +1,5 @@
+function FDTD = SetCustomExcite(FDTD,f0,funcStr);
+
+FDTD.Excitation.ATTRIBUTE.Type=10;
+FDTD.Excitation.ATTRIBUTE.f0=f0;
+FDTD.Excitation.ATTRIBUTE.Function=funcStr;

openEMS.pro

@@ -8,6 +8,7 @@ CONFIG -= app_bundle
 TEMPLATE = app
 OBJECTS_DIR = obj
 INCLUDEPATH += ../CSXCAD \
+	../fparser \
     ../tinyxml
 LIBS += -L../CSXCAD \
     -lCSXCAD \

openems.cpp

@@ -108,13 +108,55 @@ bool openEMS::parseCommandLineArgument( const char *argv )
 	return false;
 }
 
+void openEMS::SetupExcitation(TiXmlElement* Excite)
+{
+	if (Excite==NULL)
+	{
+		cerr << "Can't read openEMS Excitation Settings... " << endl;
+		exit(-2);
+	}
+
+	int Excit_Type=0;
+	double f0=0;
+	double fc=0;
+	Excite->QueryIntAttribute("Type",&Excit_Type);
+
+	unsigned int Nyquist = 0;
+	switch (Excit_Type)
+	{
+		case 0:
+			Excite->QueryDoubleAttribute("f0",&f0);
+			Excite->QueryDoubleAttribute("fc",&fc);
+			Nyquist = FDTD_Op->CalcGaussianPulsExcitation(f0,fc);
+			break;
+		case 1:
+			Excite->QueryDoubleAttribute("f0",&f0);
+			Nyquist = FDTD_Op->CalcSinusExcitation(f0,NrTS);
+			break;
+		case 2:
+			Nyquist = FDTD_Op->CalcDiracPulsExcitation();
+			break;
+		case 3:
+			Nyquist = FDTD_Op->CalcStepExcitation();
+			break;
+		case 10:
+			Excite->QueryDoubleAttribute("f0",&f0);
+			Nyquist = FDTD_Op->CalcCustomExcitation(f0,NrTS,Excite->Attribute("Function"));
+			break;
+	}
+
+	if (!Nyquist)
+	{
+		cerr << "openEMS: excitation setup failed!!" << endl;
+		exit(2);
+	}
+	FDTD_Op->SetNyquistNum(Nyquist);
+}
+
 int openEMS::SetupFDTD(const char* file)
 {
 	if (file==NULL) return -1;
 	Reset();
-	double f0=0;
-	double fc=0;
-	int Excit_Type=0;
 	int bounds[6];
 
 	time_t startTime=time(NULL);
@@ -135,6 +177,7 @@ int openEMS::SetupFDTD(const char* file)
 	}
 
 	TiXmlElement* FDTD_Opts = openEMSxml->FirstChildElement("FDTD");
+
 	if (FDTD_Opts==NULL)
 	{
 		cerr << "Can't read openEMS FDTD Settings... " << endl;
@@ -146,6 +189,7 @@ int openEMS::SetupFDTD(const char* file)
 		NrTS=0;
 	else
 		NrTS = help;
+
 	FDTD_Opts->QueryDoubleAttribute("endCriteria",&endCrit);
 	if (endCrit==0)
 		endCrit=1e-6;
@@ -154,29 +198,6 @@ int openEMS::SetupFDTD(const char* file)
 	if (m_OverSampling<2)
 		m_OverSampling=2;
 
-	TiXmlElement* Excite = FDTD_Opts->FirstChildElement("Excitation");
-	if (Excite==NULL)
-	{
-		cerr << "Can't read openEMS Excitation Settings... " << endl;
-		exit(-2);
-	}
-	Excite->QueryIntAttribute("Type",&Excit_Type);
-	if (Excit_Type==0)
-	{
-		Excite->QueryDoubleAttribute("f0",&f0);
-		Excite->QueryDoubleAttribute("fc",&fc);
-	}
-	else if (Excit_Type==1)
-	{
-		Excite->QueryDoubleAttribute("f0",&f0);
-		fc = 0;
-	}
-	else if (Excit_Type==2)
-	{
-		Excite->QueryDoubleAttribute("f0",&f0);
-		fc = 0;
-	}
-
 	TiXmlElement* BC = FDTD_Opts->FirstChildElement("BoundaryCond");
 	if (BC==NULL)
 	{
@@ -209,50 +230,8 @@ int openEMS::SetupFDTD(const char* file)
 	if (FDTD_Op->SetGeometryCSX(&CSX)==false) return(-1);
 
 	FDTD_Op->CalcECOperator();
-	unsigned int Nyquist = 0;
-	if (Excit_Type==0)
-	{
-		Nyquist = FDTD_Op->CalcGaussianPulsExcitation(f0,fc);
-		if (!Nyquist)
-		{
-			cerr << "openEMS: excitation setup failed!!" << endl;
-			exit(2);
-		}
-	}
-	else if (Excit_Type==1)
-	{
-		Nyquist = FDTD_Op->CalcSinusExcitation(f0,NrTS);
-		if (!Nyquist)
-		{
-			cerr << "openEMS: excitation setup failed!!" << endl;
-			exit(2);
-		}
-	}
-	else if (Excit_Type==2)
-	{
-		Nyquist = FDTD_Op->CalcDiracPulsExcitation();
-		if (!Nyquist)
-		{
-			cerr << "openEMS: excitation setup failed!!" << endl;
-			exit(2);
-		}
-	}
-	else if (Excit_Type==3)
-	{
-		Nyquist = FDTD_Op->CalcStepExcitation();
-		if (!Nyquist)
-		{
-			cerr << "openEMS: excitation setup failed!!" << endl;
-			exit(2);
-		}
-	}
-	else
-	{
-		cerr << "openEMS: Excitation type is unknown" << endl;
-		exit(-1);
-	}
 
-	FDTD_Op->SetNyquistNum(Nyquist);
+	SetupExcitation(FDTD_Opts->FirstChildElement("Excitation"));
 
 	if (DebugMat)
 	{
@@ -281,11 +260,11 @@ int openEMS::SetupFDTD(const char* file)
 		break;
 	}
 
-
 	time_t currTime = time(NULL);
 
 	//*************** setup processing ************//
 	cout << "Setting up processing..." << endl;
+	unsigned int Nyquist = FDTD_Op->GetNyquistNum();
 	PA = new ProcessingArray(Nyquist);
 
 	double start[3];

openems.h

@@ -21,6 +21,7 @@
 class Operator;
 class Engine;
 class ProcessingArray;
+class TiXmlElement;
 
 class openEMS
 {
@@ -43,6 +44,7 @@ public:
 	void DebugOperator() {DebugOp=true;}
 
 protected:
+	void SetupExcitation(TiXmlElement* Excite);
 	//! Number of Timesteps
 	unsigned int NrTS;
 	bool Enable_Dumps;