From 4d214c162f3c130e47cd338cc74fa635b63abb5a Mon Sep 17 00:00:00 2001
From: Thorsten Liebig <thorsten.liebig@gmx.de>
Date: Fri, 30 Jul 2010 17:02:21 +0200
Subject: [PATCH] sf_pml: user defined grading function using the fparser lib

---
 FDTD/operator_ext_pml_sf.cpp | 25 +++++++++++++++++++------
 FDTD/operator_ext_pml_sf.h   | 17 +++++++++++++++++
 2 files changed, 36 insertions(+), 6 deletions(-)

diff --git a/FDTD/operator_ext_pml_sf.cpp b/FDTD/operator_ext_pml_sf.cpp
index 766cf85..a3e760d 100644
--- a/FDTD/operator_ext_pml_sf.cpp
+++ b/FDTD/operator_ext_pml_sf.cpp
@@ -18,8 +18,8 @@
 #include "operator_ext_pml_sf.h"
 #include "engine_ext_pml_sf.h"
 #include "operator_cylinder.h"
-
 #include "tools/array_ops.h"
+#include "fparser.hh"
 
 bool Build_Split_Field_PML(Operator* op, int BC[6], int size[6])
 {
@@ -54,10 +54,16 @@ Operator_Ext_PML_SF::Operator_Ext_PML_SF(Operator* op) : Operator_Extension(op)
 
 	for (int n=0;n<6;++n)
 		m_BC[n]=0;
+
+	m_GradingFunction = new FunctionParser();
+	//default grading function
+	SetGradingFunction(" -log(1e-6)*log(2.5)/(2*dl*pow(2.5,W/dl)-1) * pow(2.5, D/dl) / Z ");
 }
 
 Operator_Ext_PML_SF::~Operator_Ext_PML_SF()
 {
+	delete m_GradingFunction;
+	m_GradingFunction = NULL;
 	DeleteOP();
 }
 
@@ -104,6 +110,16 @@ void Operator_Ext_PML_SF::DeleteOP()
 	iv[1] = NULL;
 }
 
+bool Operator_Ext_PML_SF::SetGradingFunction(string func)
+{
+	int res = m_GradingFunction->Parse(func.c_str(), "D,dl,W,Z,N");
+	if(res < 0) return true;
+
+	cerr << "Operator_Ext_PML_SF::SetGradingFunction: Warning, an error occured parsing the pml grading function (see below) ..." << endl;
+	cerr << func << "\n" << string(res, ' ') << "^\n" << m_GradingFunction->ErrorMsg() << "\n";
+	return false;
+}
+
 bool Operator_Ext_PML_SF::BuildExtension()
 {
 	if (!m_SetupDone)
@@ -240,11 +256,8 @@ double Operator_Ext_PML_SF_Plane::GetKappaGraded(double depth, double Zm) const
 	if (depth<0)
 		return 0.0;
 
-	//todo: use fparser to allow arbitrary, user-defined profiles and parameter
-	double g = 2.5;
-	double R0 = 1e-6;
-	double kappa0 = -log(R0)*log(g)/(2*m_pml_delta * pow(g,m_pml_width/m_pml_delta) -1);
-	return pow(g,depth/m_pml_delta)*kappa0 / Zm;
+	double vars[5] = {depth, m_pml_delta, m_pml_width, Zm, m_numLines[m_ny]};
+	return m_GradingFunction->Eval(vars);
 }
 
 bool Operator_Ext_PML_SF_Plane::Calc_ECPos(int nP, int n, unsigned int* pos, double* inEC) const
diff --git a/FDTD/operator_ext_pml_sf.h b/FDTD/operator_ext_pml_sf.h
index 2b537d7..a079cb6 100644
--- a/FDTD/operator_ext_pml_sf.h
+++ b/FDTD/operator_ext_pml_sf.h
@@ -21,6 +21,8 @@
 #include "operator.h"
 #include "operator_extension.h"
 
+class FunctionParser;
+
 //! Insert split field pml planes, edges and corner as necessary by the given boundary conditions
 bool Build_Split_Field_PML(Operator* op, int BC[6], int size[6]);
 
@@ -45,6 +47,19 @@ public:
 	//! This will resturn the pml parameter grading
 	virtual double GetKappaGraded(double depth, double Zm) const {UNUSED(depth);UNUSED(Zm);return 0.0;}
 
+	//! Set the grading function for the pml
+	/*!
+		Define the pml grading grading function.
+		Predefined variables in this grading function are:
+			D  = depth in the pml in meter
+			dl = mesh delta inside the pml in meter
+			W  = width (length) of the pml in meter
+			N  = number of cells for the pml
+			Z  = wave impedance at the current depth and position
+		example: SetGradingFunction("-log(1e-6)*log(2.5)/(2*dl*pow(2.5,W/dl)-1) * pow(2.5, D/dl) / Z");
+	*/
+	virtual bool SetGradingFunction(string func);
+
 	virtual bool BuildExtension();
 
 	virtual string GetExtensionName() const {return string("Split Field PML Extension");}
@@ -63,6 +78,8 @@ protected:
 
 	int m_BC[6];
 
+	FunctionParser* m_GradingFunction;
+
 	void InitOP();
 	void DeleteOP();
 	//split field EC operator