new allow multiple ProcessIntegral and ProcessModeMatch using this to dump mode purity additionally

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();
+	CalcMultipleIntegrals();
-	integral*=m_weight;
+	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;
+}

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

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;
 			}
 		}
 	}
-
+	if (purity!=0)
-	return value;
+		m_Results[1] = value*value/purity;
+	else
+		m_Results[1] = 0;
+	m_Results[0] = value;
+	return m_Results;
 }

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

matlab/ReadUI.m

@@ -43,6 +43,10 @@ for n=1:numel(filenames)
     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 );
     else