From 796fd83f7b1239dafa6d1fe8a198bcd8c6621579 Mon Sep 17 00:00:00 2001
From: Thorsten Liebig <Thorsten.Liebig@gmx.de>
Date: Tue, 24 Apr 2012 11:05:19 +0200
Subject: [PATCH] dispersive material extension: prepared for higher order
 dispersion

---
 FDTD/extensions/engine_ext_dispersive.cpp     | 158 ++++++++++--------
 FDTD/extensions/engine_ext_dispersive.h       |   9 +-
 .../extensions/engine_ext_lorentzmaterial.cpp | 135 ++++++++-------
 FDTD/extensions/operator_ext_dispersive.cpp   |  42 +++--
 FDTD/extensions/operator_ext_dispersive.h     |  19 ++-
 .../operator_ext_lorentzmaterial.cpp          | 106 ++++++++----
 .../extensions/operator_ext_lorentzmaterial.h |  11 +-
 7 files changed, 288 insertions(+), 192 deletions(-)

diff --git a/FDTD/extensions/engine_ext_dispersive.cpp b/FDTD/extensions/engine_ext_dispersive.cpp
index bc6ec81..d9c431b 100644
--- a/FDTD/extensions/engine_ext_dispersive.cpp
+++ b/FDTD/extensions/engine_ext_dispersive.cpp
@@ -22,117 +22,141 @@
 Engine_Ext_Dispersive::Engine_Ext_Dispersive(Operator_Ext_Dispersive* op_ext_disp) : Engine_Extension(op_ext_disp)
 {
 	m_Op_Ext_Disp = op_ext_disp;
-	for (int n=0; n<3; ++n)
+	int order = m_Op_Ext_Disp->m_Order;
+	curr_ADE = new FDTD_FLOAT**[order];
+	volt_ADE = new FDTD_FLOAT**[order];
+	for (int o=0;o<order;++o)
 	{
-		if (m_Op_Ext_Disp->m_curr_ADE_On==true)
+		curr_ADE[o] = new FDTD_FLOAT*[3];
+		volt_ADE[o] = new FDTD_FLOAT*[3];
+		for (int n=0; n<3; ++n)
 		{
-			curr_ADE[n] = new FDTD_FLOAT[m_Op_Ext_Disp->m_LM_Count];
-			for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
-				curr_ADE[n][i]=0.0;
+			if (m_Op_Ext_Disp->m_curr_ADE_On[o]==true)
+			{
+				curr_ADE[o][n] = new FDTD_FLOAT[m_Op_Ext_Disp->m_LM_Count[o]];
+				for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count[o]; ++i)
+					curr_ADE[o][n][i]=0.0;
+			}
+			else
+				curr_ADE[o][n] = NULL;
+			if (m_Op_Ext_Disp->m_volt_ADE_On[o]==true)
+			{
+				volt_ADE[o][n] = new FDTD_FLOAT[m_Op_Ext_Disp->m_LM_Count[o]];
+				for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count[o]; ++i)
+					volt_ADE[o][n][i]=0.0;
+			}
+			else
+				volt_ADE[o][n] = NULL;
 		}
-		else
-			curr_ADE[n] = NULL;
-		if (m_Op_Ext_Disp->m_volt_ADE_On==true)
-		{
-			volt_ADE[n] = new FDTD_FLOAT[m_Op_Ext_Disp->m_LM_Count];
-			for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
-				volt_ADE[n][i]=0.0;
-		}
-		else
-			volt_ADE[n] = NULL;
 	}
 }
 
 Engine_Ext_Dispersive::~Engine_Ext_Dispersive()
 {
-	for (int n=0; n<3; ++n)
+	if (curr_ADE==NULL && volt_ADE==NULL)
+		return;
+
+	for (int o=0;o<m_Op_Ext_Disp->m_Order;++o)
 	{
-		delete[] curr_ADE[n];
-		curr_ADE[n] = NULL;
-		delete[] volt_ADE[n];
-		volt_ADE[n] = NULL;
+		for (int n=0; n<3; ++n)
+		{
+			delete[] curr_ADE[o][n];
+			delete[] volt_ADE[o][n];
+		}
+		delete[] curr_ADE[o];
+		delete[] volt_ADE[o];
 	}
+	delete[] curr_ADE;
+	curr_ADE=NULL;
+
+	delete[] volt_ADE;
+	volt_ADE=NULL;
 }
 
 void Engine_Ext_Dispersive::Apply2Voltages()
 {
-	if (m_Op_Ext_Disp->m_volt_ADE_On==false) return;
-
-	unsigned int **pos = m_Op_Ext_Disp->m_LM_pos;
-
-	//switch for different engine types to access faster inline engine functions
-	switch (m_Eng->GetType())
+	for (int o=0;o<m_Op_Ext_Disp->m_Order;++o)
 	{
-	case Engine::BASIC:
+		if (m_Op_Ext_Disp->m_volt_ADE_On[o]==false) continue;
+
+		unsigned int **pos = m_Op_Ext_Disp->m_LM_pos[o];
+
+		//switch for different engine types to access faster inline engine functions
+		switch (m_Eng->GetType())
 		{
-			for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
+		case Engine::BASIC:
+		{
+			for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count.at(o); ++i)
 			{
-				m_Eng->Engine::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]);
-				m_Eng->Engine::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]);
-				m_Eng->Engine::SetVolt(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[2][i]);
+				m_Eng->Engine::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][0][i]);
+				m_Eng->Engine::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][1][i]);
+				m_Eng->Engine::SetVolt(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][2][i]);
 			}
 			break;
 		}
-	case Engine::SSE:
+		case Engine::SSE:
 		{
 			Engine_sse* eng_sse = (Engine_sse*)m_Eng;
-			for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
+			for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count.at(o); ++i)
 			{
-				eng_sse->Engine_sse::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]);
-				eng_sse->Engine_sse::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]);
-				eng_sse->Engine_sse::SetVolt(2,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[2][i]);
+				eng_sse->Engine_sse::SetVolt(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][0][i]);
+				eng_sse->Engine_sse::SetVolt(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][1][i]);
+				eng_sse->Engine_sse::SetVolt(2,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][2][i]);
 			}
 			break;
 		}
-	default:
-		for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
-		{
-			m_Eng->SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[0][i]);
-			m_Eng->SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[1][i]);
-			m_Eng->SetVolt(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[2][i]);
+		default:
+			for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count.at(o); ++i)
+			{
+				m_Eng->SetVolt(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][0][i]);
+				m_Eng->SetVolt(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][1][i]);
+				m_Eng->SetVolt(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i]) - volt_ADE[o][2][i]);
+			}
+			break;
 		}
-		break;
 	}
 }
 
 void Engine_Ext_Dispersive::Apply2Current()
 {
-	if (m_Op_Ext_Disp->m_curr_ADE_On==false) return;
-
-	unsigned int **pos = m_Op_Ext_Disp->m_LM_pos;
-
-	//switch for different engine types to access faster inline engine functions
-	switch (m_Eng->GetType())
+	for (int o=0;o<m_Op_Ext_Disp->m_Order;++o)
 	{
-	case Engine::BASIC:
+		if (m_Op_Ext_Disp->m_curr_ADE_On[o]==false) continue;
+
+		unsigned int **pos = m_Op_Ext_Disp->m_LM_pos[o];
+
+		//switch for different engine types to access faster inline engine functions
+		switch (m_Eng->GetType())
 		{
-			for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
+		case Engine::BASIC:
+		{
+			for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count.at(o); ++i)
 			{
-				m_Eng->Engine::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]);
-				m_Eng->Engine::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]);
-				m_Eng->Engine::SetCurr(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[2][i]);
+				m_Eng->Engine::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][0][i]);
+				m_Eng->Engine::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][1][i]);
+				m_Eng->Engine::SetCurr(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][2][i]);
 			}
 			break;
 		}
-	case Engine::SSE:
+		case Engine::SSE:
 		{
 			Engine_sse* eng_sse = (Engine_sse*)m_Eng;
-			for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
+			for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count.at(o); ++i)
 			{
-				eng_sse->Engine_sse::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]);
-				eng_sse->Engine_sse::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]);
-				eng_sse->Engine_sse::SetCurr(2,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[2][i]);
+				eng_sse->Engine_sse::SetCurr(0,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][0][i]);
+				eng_sse->Engine_sse::SetCurr(1,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][1][i]);
+				eng_sse->Engine_sse::SetCurr(2,pos[0][i],pos[1][i],pos[2][i], eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][2][i]);
 			}
 			break;
 		}
-	default:
-		for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count; ++i)
-		{
-			m_Eng->SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[0][i]);
-			m_Eng->SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[1][i]);
-			m_Eng->SetCurr(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[2][i]);
+		default:
+			for (unsigned int i=0; i<m_Op_Ext_Disp->m_LM_Count.at(o); ++i)
+			{
+				m_Eng->SetCurr(0,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][0][i]);
+				m_Eng->SetCurr(1,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][1][i]);
+				m_Eng->SetCurr(2,pos[0][i],pos[1][i],pos[2][i], m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i]) - curr_ADE[o][2][i]);
+			}
+			break;
 		}
-		break;
 	}
-
 }
diff --git a/FDTD/extensions/engine_ext_dispersive.h b/FDTD/extensions/engine_ext_dispersive.h
index 536a5bf..c1b84ea 100644
--- a/FDTD/extensions/engine_ext_dispersive.h
+++ b/FDTD/extensions/engine_ext_dispersive.h
@@ -36,11 +36,16 @@ public:
 protected:
 	Operator_Ext_Dispersive* m_Op_Ext_Disp;
 
+	//! Dispersive order
+	int m_Order;
+
 	//! ADE currents
-	FDTD_FLOAT *curr_ADE[3];
+	// Array setup: curr_ADE[N_order][direction][mesh_pos]
+	FDTD_FLOAT ***curr_ADE;
 
 	//! ADE voltages
-	FDTD_FLOAT *volt_ADE[3];
+	// Array setup: volt_ADE[N_order][direction][mesh_pos]
+	FDTD_FLOAT ***volt_ADE;
 };
 
 #endif // ENGINE_EXT_DISPERSIVE_H
diff --git a/FDTD/extensions/engine_ext_lorentzmaterial.cpp b/FDTD/extensions/engine_ext_lorentzmaterial.cpp
index cedfdca..40e006f 100644
--- a/FDTD/extensions/engine_ext_lorentzmaterial.cpp
+++ b/FDTD/extensions/engine_ext_lorentzmaterial.cpp
@@ -22,6 +22,7 @@
 Engine_Ext_LorentzMaterial::Engine_Ext_LorentzMaterial(Operator_Ext_LorentzMaterial* op_ext_lorentz) : Engine_Ext_Dispersive(op_ext_lorentz)
 {
 	m_Op_Ext_Lor = op_ext_lorentz;
+	m_Order = m_Op_Ext_Lor->GetDispersionOrder();
 }
 
 Engine_Ext_LorentzMaterial::~Engine_Ext_LorentzMaterial()
@@ -31,113 +32,119 @@ Engine_Ext_LorentzMaterial::~Engine_Ext_LorentzMaterial()
 
 void Engine_Ext_LorentzMaterial::DoPreVoltageUpdates()
 {
-	if (m_Op_Ext_Lor->m_volt_ADE_On==false) return;
-
-	unsigned int **pos = m_Op_Ext_Lor->m_LM_pos;
-
-	//switch for different engine types to access faster inline engine functions
-	switch (m_Eng->GetType())
+	for (int o=0;o<m_Order;++o)
 	{
-	case Engine::BASIC:
+		if (m_Op_Ext_Lor->m_volt_ADE_On==false) continue;
+
+		unsigned int **pos = m_Op_Ext_Lor->m_LM_pos[o];
+
+		//switch for different engine types to access faster inline engine functions
+		switch (m_Eng->GetType())
 		{
-			for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count; ++i)
+		case Engine::BASIC:
+		{
+			for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i)
 			{
-				volt_ADE[0][i] *= m_Op_Ext_Lor->v_int_ADE[0][i];
-				volt_ADE[0][i] += m_Op_Ext_Lor->v_ext_ADE[0][i] * m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]);
+				volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i];
+				volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * m_Eng->Engine::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]);
 
-				volt_ADE[1][i] *= m_Op_Ext_Lor->v_int_ADE[1][i];
-				volt_ADE[1][i] += m_Op_Ext_Lor->v_ext_ADE[1][i] * m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]);
+				volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i];
+				volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * m_Eng->Engine::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]);
 
-				volt_ADE[2][i] *= m_Op_Ext_Lor->v_int_ADE[2][i];
-				volt_ADE[2][i] += m_Op_Ext_Lor->v_ext_ADE[2][i] * m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]);
+				volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i];
+				volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * m_Eng->Engine::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]);
 			}
 			break;
 		}
-	case Engine::SSE:
+		case Engine::SSE:
 		{
 			Engine_sse* eng_sse = (Engine_sse*)m_Eng;
-			for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count; ++i)
+			for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i)
 			{
-				volt_ADE[0][i] *= m_Op_Ext_Lor->v_int_ADE[0][i];
-				volt_ADE[0][i] += m_Op_Ext_Lor->v_ext_ADE[0][i] * eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]);
+				volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i];
+				volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * eng_sse->Engine_sse::GetVolt(0,pos[0][i],pos[1][i],pos[2][i]);
 
-				volt_ADE[1][i] *= m_Op_Ext_Lor->v_int_ADE[1][i];
-				volt_ADE[1][i] += m_Op_Ext_Lor->v_ext_ADE[1][i] * eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]);
+				volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i];
+				volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * eng_sse->Engine_sse::GetVolt(1,pos[0][i],pos[1][i],pos[2][i]);
 
-				volt_ADE[2][i] *= m_Op_Ext_Lor->v_int_ADE[2][i];
-				volt_ADE[2][i] += m_Op_Ext_Lor->v_ext_ADE[2][i] * eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]);
+				volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i];
+				volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * eng_sse->Engine_sse::GetVolt(2,pos[0][i],pos[1][i],pos[2][i]);
 			}
 			break;
 		}
-	default:
-		for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count; ++i)
-		{
-			volt_ADE[0][i] *= m_Op_Ext_Lor->v_int_ADE[0][i];
-			volt_ADE[0][i] += m_Op_Ext_Lor->v_ext_ADE[0][i] * m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]);
+		default:
+			for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i)
+			{
+				volt_ADE[o][0][i] *= m_Op_Ext_Lor->v_int_ADE[o][0][i];
+				volt_ADE[o][0][i] += m_Op_Ext_Lor->v_ext_ADE[o][0][i] * m_Eng->GetVolt(0,pos[0][i],pos[1][i],pos[2][i]);
 
-			volt_ADE[1][i] *= m_Op_Ext_Lor->v_int_ADE[1][i];
-			volt_ADE[1][i] += m_Op_Ext_Lor->v_ext_ADE[1][i] * m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]);
+				volt_ADE[o][1][i] *= m_Op_Ext_Lor->v_int_ADE[o][1][i];
+				volt_ADE[o][1][i] += m_Op_Ext_Lor->v_ext_ADE[o][1][i] * m_Eng->GetVolt(1,pos[0][i],pos[1][i],pos[2][i]);
 
-			volt_ADE[2][i] *= m_Op_Ext_Lor->v_int_ADE[2][i];
-			volt_ADE[2][i] += m_Op_Ext_Lor->v_ext_ADE[2][i] * m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i]);
+				volt_ADE[o][2][i] *= m_Op_Ext_Lor->v_int_ADE[o][2][i];
+				volt_ADE[o][2][i] += m_Op_Ext_Lor->v_ext_ADE[o][2][i] * m_Eng->GetVolt(2,pos[0][i],pos[1][i],pos[2][i]);
+			}
+			break;
 		}
-		break;
 	}
 }
 
 void Engine_Ext_LorentzMaterial::DoPreCurrentUpdates()
 {
-	if (m_Op_Ext_Lor->m_curr_ADE_On==false) return;
-
-	unsigned int **pos = m_Op_Ext_Lor->m_LM_pos;
-
-	//switch for different engine types to access faster inline engine functions
-	switch (m_Eng->GetType())
+	for (int o=0;o<m_Order;++o)
 	{
-	case Engine::BASIC:
+		if (m_Op_Ext_Lor->m_curr_ADE_On==false) continue;
+
+		unsigned int **pos = m_Op_Ext_Lor->m_LM_pos[o];
+
+		//switch for different engine types to access faster inline engine functions
+		switch (m_Eng->GetType())
 		{
-			for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count; ++i)
+		case Engine::BASIC:
+		{
+			for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i)
 			{
-				curr_ADE[0][i] *= m_Op_Ext_Lor->i_int_ADE[0][i];
-				curr_ADE[0][i] += m_Op_Ext_Lor->i_ext_ADE[0][i] * m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]);
+				curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i];
+				curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * m_Eng->Engine::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]);
 
-				curr_ADE[1][i] *= m_Op_Ext_Lor->i_int_ADE[1][i];
-				curr_ADE[1][i] += m_Op_Ext_Lor->i_ext_ADE[1][i] * m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]);
+				curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i];
+				curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * m_Eng->Engine::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]);
 
-				curr_ADE[2][i] *= m_Op_Ext_Lor->i_int_ADE[2][i];
-				curr_ADE[2][i] += m_Op_Ext_Lor->i_ext_ADE[2][i] * m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]);
+				curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i];
+				curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * m_Eng->Engine::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]);
 			}
 			break;
 		}
-	case Engine::SSE:
+		case Engine::SSE:
 		{
 			Engine_sse* eng_sse = (Engine_sse*)m_Eng;
-			for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count; ++i)
+			for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i)
 			{
-				curr_ADE[0][i] *= m_Op_Ext_Lor->i_int_ADE[0][i];
-				curr_ADE[0][i] += m_Op_Ext_Lor->i_ext_ADE[0][i] * eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]);
+				curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i];
+				curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * eng_sse->Engine_sse::GetCurr(0,pos[0][i],pos[1][i],pos[2][i]);
 
-				curr_ADE[1][i] *= m_Op_Ext_Lor->i_int_ADE[1][i];
-				curr_ADE[1][i] += m_Op_Ext_Lor->i_ext_ADE[1][i] * eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]);
+				curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i];
+				curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * eng_sse->Engine_sse::GetCurr(1,pos[0][i],pos[1][i],pos[2][i]);
 
-				curr_ADE[2][i] *= m_Op_Ext_Lor->i_int_ADE[2][i];
-				curr_ADE[2][i] += m_Op_Ext_Lor->i_ext_ADE[2][i] * eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]);
+				curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i];
+				curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * eng_sse->Engine_sse::GetCurr(2,pos[0][i],pos[1][i],pos[2][i]);
 			}
 			break;
 		}
-	default:
-		for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count; ++i)
-		{
-			curr_ADE[0][i] *= m_Op_Ext_Lor->i_int_ADE[0][i];
-			curr_ADE[0][i] += m_Op_Ext_Lor->i_ext_ADE[0][i] * m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]);
+		default:
+			for (unsigned int i=0; i<m_Op_Ext_Lor->m_LM_Count.at(o); ++i)
+			{
+				curr_ADE[o][0][i] *= m_Op_Ext_Lor->i_int_ADE[o][0][i];
+				curr_ADE[o][0][i] += m_Op_Ext_Lor->i_ext_ADE[o][0][i] * m_Eng->GetCurr(0,pos[0][i],pos[1][i],pos[2][i]);
 
-			curr_ADE[1][i] *= m_Op_Ext_Lor->i_int_ADE[1][i];
-			curr_ADE[1][i] += m_Op_Ext_Lor->i_ext_ADE[1][i] * m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]);
+				curr_ADE[o][1][i] *= m_Op_Ext_Lor->i_int_ADE[o][1][i];
+				curr_ADE[o][1][i] += m_Op_Ext_Lor->i_ext_ADE[o][1][i] * m_Eng->GetCurr(1,pos[0][i],pos[1][i],pos[2][i]);
 
-			curr_ADE[2][i] *= m_Op_Ext_Lor->i_int_ADE[2][i];
-			curr_ADE[2][i] += m_Op_Ext_Lor->i_ext_ADE[2][i] * m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i]);
+				curr_ADE[o][2][i] *= m_Op_Ext_Lor->i_int_ADE[o][2][i];
+				curr_ADE[o][2][i] += m_Op_Ext_Lor->i_ext_ADE[o][2][i] * m_Eng->GetCurr(2,pos[0][i],pos[1][i],pos[2][i]);
+			}
+			break;
 		}
-		break;
 	}
 }
 
diff --git a/FDTD/extensions/operator_ext_dispersive.cpp b/FDTD/extensions/operator_ext_dispersive.cpp
index 95b9525..56c333f 100644
--- a/FDTD/extensions/operator_ext_dispersive.cpp
+++ b/FDTD/extensions/operator_ext_dispersive.cpp
@@ -21,30 +21,44 @@
 
 Operator_Ext_Dispersive::Operator_Ext_Dispersive(Operator* op)  : Operator_Extension(op)
 {
-	m_curr_ADE_On = false;
-	m_volt_ADE_On = false;
+	m_curr_ADE_On = NULL;
+	m_volt_ADE_On = NULL;
 
-	m_LM_pos[0]=NULL;
-	m_LM_pos[1]=NULL;
-	m_LM_pos[2]=NULL;
+	m_LM_pos=NULL;
+	m_curr_ADE_On=NULL;
+	m_volt_ADE_On=NULL;
+
+	m_Order = 0;
 }
 
 Operator_Ext_Dispersive::~Operator_Ext_Dispersive()
 {
-	delete[] m_LM_pos[0];
-	delete[] m_LM_pos[1];
-	delete[] m_LM_pos[2];
+	delete[] m_curr_ADE_On;
+	delete[] m_volt_ADE_On;
+	m_curr_ADE_On=NULL;
+	m_volt_ADE_On=NULL;
 
-	m_LM_pos[0]=NULL;
-	m_LM_pos[1]=NULL;
-	m_LM_pos[2]=NULL;
+	for (int n=0;n<m_Order;++n)
+	{
+		delete[] m_LM_pos[n][0];
+		delete[] m_LM_pos[n][1];
+		delete[] m_LM_pos[n][2];
+	}
+	delete[] m_LM_pos;
+	m_LM_pos=NULL;
+	m_Order=0;
+	m_LM_Count.clear();
 }
 
 void Operator_Ext_Dispersive::ShowStat(ostream &ostr)  const
 {
 	Operator_Extension::ShowStat(ostr);
 	string On_Off[2] = {"Off", "On"};
-	ostr << " Active cells\t\t: " << 	m_LM_Count << endl;
-	ostr << " Voltage ADE is \t: " << On_Off[m_volt_ADE_On] << endl;
-	ostr << " Current ADE is \t: " << On_Off[m_curr_ADE_On] << endl;
+	ostr << " Max. Dispersion Order N = " << m_Order << endl;
+	for (int i=0;i<m_Order;++i)
+	{
+		ostr << " N=" << i << ":\t Active cells\t\t: " << 	m_LM_Count[i] << endl;
+		ostr << " N=" << i << ":\t Voltage ADE is \t: " << On_Off[m_volt_ADE_On[i]] << endl;
+		ostr << " N=" << i << ":\t Current ADE is \t: " << On_Off[m_curr_ADE_On[i]] << endl;
+	}
 }
diff --git a/FDTD/extensions/operator_ext_dispersive.h b/FDTD/extensions/operator_ext_dispersive.h
index e647632..02ca3fc 100644
--- a/FDTD/extensions/operator_ext_dispersive.h
+++ b/FDTD/extensions/operator_ext_dispersive.h
@@ -20,6 +20,7 @@
 
 //#include "operator.h"
 #include "operator_extension.h"
+#include "vector"
 
 //! Abstract base class for all dispersive material models, based on an ADE (additional differential equation)
 class Operator_Ext_Dispersive : public Operator_Extension
@@ -28,6 +29,8 @@ class Operator_Ext_Dispersive : public Operator_Extension
 public:
 	virtual ~Operator_Ext_Dispersive();
 
+	virtual int GetDispersionOrder() {return m_Order;}
+
 	virtual string GetExtensionName() const {return string("Dispersive Material Abstract Base class");}
 
 	virtual void ShowStat(ostream &ostr) const;
@@ -35,13 +38,17 @@ public:
 protected:
 	Operator_Ext_Dispersive(Operator* op);
 
-	//! Lorentz material count
-	unsigned int m_LM_Count;
-	//! Index with dispersive material
-	unsigned int *m_LM_pos[3];
+	//! Dispersive order
+	int m_Order;
 
-	bool m_curr_ADE_On;
-	bool m_volt_ADE_On;
+	//! Dispersive material count
+	vector<unsigned int> m_LM_Count;
+	//! Index with dispersive material
+	// Array setup: m_LM_pos[N_order][direction][mesh_pos]
+	unsigned int ***m_LM_pos;
+
+	bool *m_curr_ADE_On;
+	bool *m_volt_ADE_On;
 };
 
 #endif // OPERATOR_EXT_DISPERSIVE_H
diff --git a/FDTD/extensions/operator_ext_lorentzmaterial.cpp b/FDTD/extensions/operator_ext_lorentzmaterial.cpp
index c272204..b0645bb 100644
--- a/FDTD/extensions/operator_ext_lorentzmaterial.cpp
+++ b/FDTD/extensions/operator_ext_lorentzmaterial.cpp
@@ -20,28 +20,48 @@
 
 Operator_Ext_LorentzMaterial::Operator_Ext_LorentzMaterial(Operator* op) : Operator_Ext_Dispersive(op)
 {
-	for (int n=0; n<3; ++n)
-	{
-		v_int_ADE[n] = NULL;
-		v_ext_ADE[n] = NULL;
-		i_int_ADE[n] = NULL;
-		i_ext_ADE[n] = NULL;
-	}
+	v_int_ADE = NULL;
+	v_ext_ADE = NULL;
+	i_int_ADE = NULL;
+	i_ext_ADE = NULL;
 }
 
 Operator_Ext_LorentzMaterial::~Operator_Ext_LorentzMaterial()
 {
-	for (int n=0; n<3; ++n)
+	for (int i=0;i<m_Order;++i)
 	{
-		delete[] v_int_ADE[n];
-		v_int_ADE[n] = NULL;
-		delete[] v_ext_ADE[n];
-		v_ext_ADE[n] = NULL;
-		delete[] i_int_ADE[n];
-		i_int_ADE[n] = NULL;
-		delete[] i_ext_ADE[n];
-		i_ext_ADE[n] = NULL;
+		for (int n=0; n<3; ++n)
+		{
+			if (m_volt_ADE_On[i])
+			{
+				delete[] v_int_ADE[i][n];
+				delete[] v_ext_ADE[i][n];
+			}
+			if (m_curr_ADE_On[i])
+			{
+				delete[] i_int_ADE[i][n];
+				delete[] i_ext_ADE[i][n];
+			}
+		}
+		if (m_volt_ADE_On[i])
+		{
+			delete[] v_int_ADE[i];
+			delete[] v_ext_ADE[i];
+		}
+		if (m_curr_ADE_On[i])
+		{
+			delete[] i_int_ADE[i];
+			delete[] i_ext_ADE[i];
+		}
 	}
+	delete[] v_int_ADE;
+	delete[] v_ext_ADE;
+	delete[] i_int_ADE;
+	delete[] i_ext_ADE;
+	v_int_ADE = NULL;
+	v_ext_ADE = NULL;
+	i_int_ADE = NULL;
+	i_ext_ADE = NULL;
 }
 
 bool Operator_Ext_LorentzMaterial::BuildExtension()
@@ -61,6 +81,11 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
 	vector<double> i_int[3];
 	vector<double> i_ext[3];
 	vector<unsigned int> v_pos[3];
+	m_Order = 1;
+	m_volt_ADE_On = new bool[1];
+	m_volt_ADE_On[0]=false;
+	m_curr_ADE_On = new bool[1];
+	m_curr_ADE_On[0]=false;
 
 	for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
 	{
@@ -87,11 +112,11 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
 						if (w_plasma>0)
 						{
 							b_pos_on = true;
-							m_volt_ADE_On = true;
+							m_volt_ADE_On[0] = true;
 							L_D[n] = 1/(w_plasma*w_plasma*m_Op->EC_C[n][index]);
 						}
 						t_relax = mat->GetEpsRelaxTimeWeighted(n,coord);
-						if ((t_relax>0) && m_volt_ADE_On)
+						if ((t_relax>0) && m_volt_ADE_On[0])
 						{
 							R_D[n] = L_D[n]/t_relax;
 						}
@@ -114,11 +139,11 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
 						if (w_plasma>0)
 						{
 							b_pos_on = true;
-							m_curr_ADE_On = true;
+							m_curr_ADE_On[0] = true;
 							C_D[n] = 1/(w_plasma*w_plasma*m_Op->EC_L[n][index]);
 						}
 						t_relax = mat->GetMueRelaxTimeWeighted(n,coord);
-						if ((t_relax>0) && m_curr_ADE_On)
+						if ((t_relax>0) && m_curr_ADE_On[0])
 						{
 							G_D[n] = C_D[n]/t_relax;
 						}
@@ -158,33 +183,46 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
 	}
 
 	//copy all vectors into the array's
-	m_LM_Count = v_pos[0].size();
+	m_LM_Count.push_back(v_pos[0].size());
+
+	m_LM_pos = new unsigned int**[1];
+	m_LM_pos[0] = new unsigned int*[3];
+
+	v_int_ADE = new FDTD_FLOAT**[1];
+	v_ext_ADE = new FDTD_FLOAT**[1];
+	i_int_ADE = new FDTD_FLOAT**[1];
+	i_ext_ADE = new FDTD_FLOAT**[1];
+
+	v_int_ADE[0] = new FDTD_FLOAT*[3];
+	v_ext_ADE[0] = new FDTD_FLOAT*[3];
+	i_int_ADE[0] = new FDTD_FLOAT*[3];
+	i_ext_ADE[0] = new FDTD_FLOAT*[3];
 
 	for (int n=0; n<3; ++n)
 	{
-		m_LM_pos[n] = new unsigned int[m_LM_Count];
-		for (unsigned int i=0; i<m_LM_Count; ++i)
-			m_LM_pos[n][i] = v_pos[n].at(i);
+		m_LM_pos[0][n] = new unsigned int[m_LM_Count.at(0)];
+		for (unsigned int i=0; i<m_LM_Count.at(0); ++i)
+			m_LM_pos[0][n][i] = v_pos[n].at(i);
 		if (m_volt_ADE_On)
 		{
-			v_int_ADE[n]  = new FDTD_FLOAT[m_LM_Count];
-			v_ext_ADE[n]  = new FDTD_FLOAT[m_LM_Count];
+			v_int_ADE[0][n]  = new FDTD_FLOAT[m_LM_Count.at(0)];
+			v_ext_ADE[0][n]  = new FDTD_FLOAT[m_LM_Count.at(0)];
 
-			for (unsigned int i=0; i<m_LM_Count; ++i)
+			for (unsigned int i=0; i<m_LM_Count.at(0); ++i)
 			{
-				v_int_ADE[n][i] = v_int[n].at(i);
-				v_ext_ADE[n][i] = v_ext[n].at(i);
+				v_int_ADE[0][n][i] = v_int[n].at(i);
+				v_ext_ADE[0][n][i] = v_ext[n].at(i);
 			}
 		}
 		if (m_curr_ADE_On)
 		{
-			i_int_ADE[n]  = new FDTD_FLOAT[m_LM_Count];
-			i_ext_ADE[n]  = new FDTD_FLOAT[m_LM_Count];
+			i_int_ADE[0][n]  = new FDTD_FLOAT[m_LM_Count.at(0)];
+			i_ext_ADE[0][n]  = new FDTD_FLOAT[m_LM_Count.at(0)];
 
-			for (unsigned int i=0; i<m_LM_Count; ++i)
+			for (unsigned int i=0; i<m_LM_Count.at(0); ++i)
 			{
-				i_int_ADE[n][i] = i_int[n].at(i);
-				i_ext_ADE[n][i] = i_ext[n].at(i);
+				i_int_ADE[0][n][i] = i_int[n].at(i);
+				i_ext_ADE[0][n][i] = i_ext[n].at(i);
 			}
 		}
 
diff --git a/FDTD/extensions/operator_ext_lorentzmaterial.h b/FDTD/extensions/operator_ext_lorentzmaterial.h
index 9cea68f..66a10d8 100644
--- a/FDTD/extensions/operator_ext_lorentzmaterial.h
+++ b/FDTD/extensions/operator_ext_lorentzmaterial.h
@@ -34,16 +34,17 @@ public:
 
 	virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {return false;}
 
-	virtual std::string GetExtensionName() const {return std::string("Lorentz Dispersive Material Extension");}
+	virtual std::string GetExtensionName() const {return std::string("Drude/Lorentz Dispersive Material Extension");}
 
 	virtual void ShowStat(ostream &ostr) const;
 
 protected:
 
-	FDTD_FLOAT *v_int_ADE[3];
-	FDTD_FLOAT *v_ext_ADE[3];
-	FDTD_FLOAT *i_int_ADE[3];
-	FDTD_FLOAT *i_ext_ADE[3];
+	//ADE update coefficients, array setup: coeff[N_order][direction][mesh_pos_index]
+	FDTD_FLOAT ***v_int_ADE;
+	FDTD_FLOAT ***v_ext_ADE;
+	FDTD_FLOAT ***i_int_ADE;
+	FDTD_FLOAT ***i_ext_ADE;
 };
 
 #endif // OPERATOR_EXT_LORENTZMATERIAL_H