From bac2fc22f7f81d3272144e52ef24d018e05fb63d Mon Sep 17 00:00:00 2001
From: Thorsten Liebig <Thorsten.Liebig@gmx.de>
Date: Tue, 10 Aug 2010 07:50:53 +0200
Subject: [PATCH] Mur-ABC: set a phase velocity to handle dispersive waveguides

---
 FDTD/operator_ext_mur_abc.cpp | 21 +++++++++++++++++----
 FDTD/operator_ext_mur_abc.h   |  5 +++++
 matlab/SetBoundaryCond.m      |  6 ++++++
 openems.cpp                   |  3 +++
 4 files changed, 31 insertions(+), 4 deletions(-)

diff --git a/FDTD/operator_ext_mur_abc.cpp b/FDTD/operator_ext_mur_abc.cpp
index 083b823..d1f37f6 100644
--- a/FDTD/operator_ext_mur_abc.cpp
+++ b/FDTD/operator_ext_mur_abc.cpp
@@ -28,6 +28,8 @@ Operator_Ext_Mur_ABC::Operator_Ext_Mur_ABC(Operator* op) : Operator_Extension(op
 	m_LineNr = 0;
 	m_LineNr_Shift = 0;
 
+	m_v_phase = 0.0;
+
 	m_Mur_Coeff_nyP = NULL;
 	m_Mur_Coeff_nyPP = NULL;
 
@@ -111,19 +113,28 @@ bool Operator_Ext_Mur_ABC::BuildExtension()
 				//nP
 				eps = mat->GetEpsilonWeighted(m_nyP,coord);
 				mue = mat->GetMueWeighted(m_nyP,coord);
-				c0t = __C0__ * dT / sqrt(eps*mue);
+				if (m_v_phase>0.0)
+					c0t = m_v_phase * dT;
+				else
+					c0t = __C0__ * dT / sqrt(eps*mue);
 				m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] = (c0t - delta) / (c0t + delta);
 
 				//nPP
 				eps = mat->GetEpsilonWeighted(m_nyPP,coord);
 				mue = mat->GetMueWeighted(m_nyPP,coord);
-				c0t = __C0__ * dT / sqrt(eps*mue);
+				if (m_v_phase>0.0)
+					c0t = m_v_phase * dT;
+				else
+					c0t = __C0__ * dT / sqrt(eps*mue);
 				m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] = (c0t - delta) / (c0t + delta);
 
 			}
 			else
 			{
-				c0t = __C0__ * dT;
+				if (m_v_phase>0.0)
+					c0t = m_v_phase * dT;
+				else
+					c0t = __C0__ * dT;
 				m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]] = (c0t - delta) / (c0t + delta);
 				m_Mur_Coeff_nyPP[pos[m_nyP]][pos[m_nyPP]] = m_Mur_Coeff_nyP[pos[m_nyP]][pos[m_nyPP]];
 			}
@@ -145,5 +156,7 @@ void Operator_Ext_Mur_ABC::ShowStat(ostream &ostr)  const
 {
 	Operator_Extension::ShowStat(ostr);
 	string XYZ[3] = {"x","y","z"};
-	ostr << " Active direction: " << XYZ[m_ny] << " at line: " << m_LineNr << endl;
+	ostr << " Active direction\t: " << XYZ[m_ny] << " at line: " << m_LineNr << endl;
+	if (m_v_phase>0.0)
+		ostr << " Used phase velocity\t: " << m_v_phase << " (" << m_v_phase/__C0__ << " * c_0)" <<endl;
 }
diff --git a/FDTD/operator_ext_mur_abc.h b/FDTD/operator_ext_mur_abc.h
index df8291b..0d6b340 100644
--- a/FDTD/operator_ext_mur_abc.h
+++ b/FDTD/operator_ext_mur_abc.h
@@ -31,6 +31,9 @@ public:
 	//! Define the direction of this ABC: \a ny=0,1,2 -> x,y,z and if at bottom_ny -> e.g. x=0 or x=end
 	void SetDirection(int ny, bool top_ny);
 
+	//! Set (override) the expected phase velocity of the incoming wave
+	void SetPhaseVelocity(double c_phase) {m_v_phase=c_phase;};
+
 	virtual bool BuildExtension();
 
 	virtual Engine_Extension* CreateEngineExtention();
@@ -47,6 +50,8 @@ protected:
 	unsigned int m_LineNr;
 	int m_LineNr_Shift;
 
+	double m_v_phase;
+
 	unsigned int m_numLines[2];
 
 	FDTD_FLOAT** m_Mur_Coeff_nyP;
diff --git a/matlab/SetBoundaryCond.m b/matlab/SetBoundaryCond.m
index f4f3ec1..fbfbd4c 100644
--- a/matlab/SetBoundaryCond.m
+++ b/matlab/SetBoundaryCond.m
@@ -13,6 +13,12 @@ function FDTD = SetBoundaryCond(FDTD, BC, varargin)
 % BC = [ 1     1     0     0     2     3     ]  %using numbers or
 % BC = {'PMC' 'PMC' 'PEC' 'PEC' 'MUR' 'PML_8'}  %usign equivalent strings
 %
+% mur-abc definitions
+% define a phase-velocity to be used by the mur-abc
+% useful e.g. for dispersive waveguides
+% FDTD = SetBoundaryCond(FDTD,BC,'MUR_PhaseVelocity',299792457.93272);
+% 
+% 
 % pml definitions
 % 	arguments:  'PML_Grading','gradFunction'
 % 		Define the pml grading grading function.
diff --git a/openems.cpp b/openems.cpp
index bdfb118..e4453a1 100644
--- a/openems.cpp
+++ b/openems.cpp
@@ -204,6 +204,9 @@ bool openEMS::SetupBoundaryConditions(TiXmlElement* BC)
 		{
 			Operator_Ext_Mur_ABC* op_ext_mur = new Operator_Ext_Mur_ABC(FDTD_Op);
 			op_ext_mur->SetDirection(n/2,n%2);
+			double v_ph = 0;
+			if (BC->QueryDoubleAttribute("MUR_PhaseVelocity",&v_ph) == TIXML_SUCCESS)
+				op_ext_mur->SetPhaseVelocity(v_ph);
 			FDTD_Op->AddExtension(op_ext_mur);
 		}
 	}