From ab701c4a7fe0c581af4d5d8f48c02a15da815a4e Mon Sep 17 00:00:00 2001
From: Thorsten Liebig <Thorsten.Liebig@gmx.de>
Date: Thu, 9 Sep 2010 10:02:33 +0200
Subject: [PATCH] new allow multiple ProcessIntegral and ProcessModeMatch using
 this to dump mode purity additionally

---
 FDTD/processintegral.cpp  | 20 +++++++++++++++++---
 FDTD/processintegral.h    | 15 +++++++++++++++
 FDTD/processmodematch.cpp | 19 +++++++++++++++----
 FDTD/processmodematch.h   |  4 +++-
 matlab/ReadUI.m           |  6 +++++-
 5 files changed, 55 insertions(+), 9 deletions(-)

diff --git a/FDTD/processintegral.cpp b/FDTD/processintegral.cpp
index ddf8b6d..c4478b5 100644
--- a/FDTD/processintegral.cpp
+++ b/FDTD/processintegral.cpp
@@ -22,10 +22,13 @@
 ProcessIntegral::ProcessIntegral(Operator* op, Engine* eng)  : Processing(op, eng)
 {
 	m_TimeShift = 0.0;
+	m_Results=NULL;
 }
 
 ProcessIntegral::~ProcessIntegral()
 {
+	delete[] m_Results;
+	m_Results = NULL;
 	ProcessIntegral::FlushData();
 }
 
@@ -50,8 +53,9 @@ int ProcessIntegral::Process()
 	if (Enabled==false) return -1;
 	if (CheckTimestep()==false) return GetNextInterval();
 
-	FDTD_FLOAT integral=CalcIntegral();
-	integral*=m_weight;
+	CalcMultipleIntegrals();
+	int NrInt = GetNumberOfIntegrals();
+	double integral = m_Results[0] * m_weight;
 
 	double time = (double)Eng->GetNumberOfTimesteps()*Op->GetTimestep() + m_TimeShift;
 
@@ -60,7 +64,10 @@ int ProcessIntegral::Process()
 		if (Eng->GetNumberOfTimesteps()%ProcessInterval==0)
 		{
 			TD_Values.push_back(integral);
-			file << setprecision(m_precision) << time << "\t" << integral << endl;
+			file << setprecision(m_precision) << time;
+			for (int n=0;n<NrInt;++n)
+				file << "\t" << m_Results[n] * m_weight;
+			file << endl;
 		}
 	}
 
@@ -83,3 +90,10 @@ int ProcessIntegral::Process()
 	return GetNextInterval();
 }
 
+double* ProcessIntegral::CalcMultipleIntegrals()
+{
+	if (m_Results==NULL)
+		m_Results = new double[1];
+	m_Results[0] = CalcIntegral();
+	return m_Results;
+}
diff --git a/FDTD/processintegral.h b/FDTD/processintegral.h
index 5c99407..a5d8a35 100644
--- a/FDTD/processintegral.h
+++ b/FDTD/processintegral.h
@@ -21,6 +21,9 @@
 #include "processing.h"
 
 //! Abstract base class for integral parameter processing
+/*!
+  \todo Weighting is applied equally to all integral parameter --> todo: weighting for each result individually
+  */
 class ProcessIntegral : public Processing
 {
 public:
@@ -31,6 +34,16 @@ public:
 	//! Flush FD data to disk
 	virtual void FlushData();
 
+	//! This method can calculate multiple integral parameter and must be overloaded for each derived class. \sa GetNumberOfIntegrals
+	/*!
+	  This method will store its integral results internally with a size given by GetNumberOfIntegrals()
+	  It will return the result for the CalcIntegral() as default.
+	  */
+	virtual double* CalcMultipleIntegrals();
+
+	//! Number of calculated results produced by this integral processing. \sa CalcMultipleIntegrals
+	virtual int GetNumberOfIntegrals() const {return 1;}
+
 	//! This method should calculate the integral parameter and must be overloaded for each derived class
 	virtual double CalcIntegral() {return 0;}
 
@@ -45,6 +58,8 @@ protected:
 
 	vector<FDTD_FLOAT> TD_Values;
 	vector<_Complex double> FD_Values;
+
+	double *m_Results;
 };
 
 #endif // PROCESSINTEGRAL_H
diff --git a/FDTD/processmodematch.cpp b/FDTD/processmodematch.cpp
index 546fe40..0239f97 100644
--- a/FDTD/processmodematch.cpp
+++ b/FDTD/processmodematch.cpp
@@ -27,6 +27,9 @@ ProcessModeMatch::ProcessModeMatch(Operator* op, Engine* eng) : ProcessIntegral(
 		m_ModeDist[n] = NULL;
 	}
 	m_dualMesh = false;
+
+	delete[] m_Results;
+	m_Results = new double[2];
 }
 
 ProcessModeMatch::~ProcessModeMatch()
@@ -232,9 +235,11 @@ double ProcessModeMatch::GetHField(int ny, unsigned int pos[3])
 }
 
 
-double ProcessModeMatch::CalcIntegral()
+double* ProcessModeMatch::CalcMultipleIntegrals()
 {
 	double value = 0;
+	double field = 0;
+	double purity = 0;
 	double area = 0;
 
 	int nP = (m_ny+1)%3;
@@ -253,10 +258,16 @@ double ProcessModeMatch::CalcIntegral()
 
 			for (int n=0;n<2;++n)
 			{
-				value += GetField((m_ny+n+1)%3,pos) * m_ModeDist[n][posP][posPP] * area;
+				field = GetField((m_ny+n+1)%3,pos);
+				value += field * m_ModeDist[n][posP][posPP] * area;
+				purity += field*field * area;
 			}
 		}
 	}
-
-	return value;
+	if (purity!=0)
+		m_Results[1] = value*value/purity;
+	else
+		m_Results[1] = 0;
+	m_Results[0] = value;
+	return m_Results;
 }
diff --git a/FDTD/processmodematch.h b/FDTD/processmodematch.h
index d35e7be..2a36e70 100644
--- a/FDTD/processmodematch.h
+++ b/FDTD/processmodematch.h
@@ -33,7 +33,9 @@ public:
 
 	void SetFieldType(int type);
 	void SetModeFunction(int ny, string function);
-	virtual double CalcIntegral();
+
+	virtual int GetNumberOfIntegrals() const {return 2;}
+	virtual double* CalcMultipleIntegrals();
 
 protected:
 	//normal direction of the mode plane
diff --git a/matlab/ReadUI.m b/matlab/ReadUI.m
index e519acb..fed2826 100644
--- a/matlab/ReadUI.m
+++ b/matlab/ReadUI.m
@@ -39,9 +39,13 @@ for n=1:numel(filenames)
     tmp = load( fullfile(path,filenames{n}) );
     t = tmp(:,1)';
     val = tmp(:,2)';
-
+    
     UI.TD{n}.t = t;
     UI.TD{n}.val = val;
+
+    if (numel(tmp(1,:))>2)
+        UI.TD{n}.additional = tmp(:,3:end)';
+    end    
     
     if (nargin<3)
         [UI.FD{n}.f,UI.FD{n}.val] = FFT_time2freq( t,val );