lumpedRLC: replace all "uint" to "unsigned int".

The original contributor used "uint" in all code, which is a
non-standard language extension that doesn't exist in ISO C
or ISO C++, and causes build failures on macOS as reported by
PR #144. Replace all "uint" with "unsigned int" for standard
conformance to maximize portability.

Signed-off-by: Yifeng Li <tomli@tomli.me>

.gitignore

@@ -26,6 +26,8 @@ localConfig.cmake
 *.pyc
 *.pyo
 python/**/*.cpp
+python/openEMS.egg-info/*
+python/dist
 !python/doc
 python/doc/_build
 python/doc/Tutorials/__*

CMakeLists.txt

@@ -6,13 +6,13 @@ ELSE()
   SET( CMAKE_BUILD_TYPE Release CACHE STRING "Set to either \"Release\" or \"Debug\"" )
 ENDIF()
 
-PROJECT(openEMS CXX)
-cmake_minimum_required(VERSION 2.8)
+PROJECT(openEMS CXX C)
+cmake_minimum_required(VERSION 3.0)
 
 # default
 set(LIB_VERSION_MAJOR 0)
 set(LIB_VERSION_MINOR 0)
-set(LIB_VERSION_PATCH 34)
+set(LIB_VERSION_PATCH 36)
 set(LIB_VERSION_STRING ${LIB_VERSION_MAJOR}.${LIB_VERSION_MINOR}.${LIB_VERSION_PATCH})
 
 set(VERSION "v${LIB_VERSION_STRING}")
@@ -21,7 +21,9 @@ IF(EXISTS ${PROJECT_SOURCE_DIR}/localConfig.cmake)
    include(${PROJECT_SOURCE_DIR}/localConfig.cmake)
 ENDIF()
 
-set(VERSION "v0.0.34")
+#ADD_DEFINITIONS( -D__SSE2__ )
+
+set(VERSION "v0.0.36")
 
 # add git revision
 IF(EXISTS ${PROJECT_SOURCE_DIR}/.git )
@@ -59,7 +61,6 @@ if(UNIX AND ENABLE_RPATH)
   set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
 endif()
 
-
 if (WITH_MPI)
 	ADD_DEFINITIONS(-DMPI_SUPPORT)
 	# Require MPI for this project:
@@ -106,11 +107,18 @@ INCLUDE_DIRECTORIES( ${CSXCAD_INCLUDE_DIR} )
 # TinyXML module from https://github.com/ros/cmake_modules
 find_package(TinyXML REQUIRED)
 ADD_DEFINITIONS( -DTIXML_USE_STL )
+message(STATUS "TinyXML_INCLUDE_DIR: ${TinyXML_INCLUDE_DIR}")
+message(STATUS "TinyXML_LIBRARY: ${TinyXML_LIBRARY}")
+
+INCLUDE_DIRECTORIES( ${TinyXML_INCLUDE_DIR} )
 
 # hdf5
 find_package(HDF5 1.8 COMPONENTS C HL REQUIRED)
 INCLUDE_DIRECTORIES (${HDF5_INCLUDE_DIRS})
 link_directories(${HDF5_LIBRARIES})
+# hdf5 compat
+#ADD_DEFINITIONS( -DH5_USE_16_API )
+#ADD_DEFINITIONS( -DH5_BUILT_AS_DYNAMIC_LIB )
 
 # boost
 find_package(Boost 1.46 COMPONENTS
@@ -119,34 +127,53 @@ find_package(Boost 1.46 COMPONENTS
   date_time
   serialization
   chrono
+  REQUIRED
 )
+message(STATUS "Boost_INCLUDE_DIRS: ${Boost_INCLUDE_DIRS}")
+message(STATUS "Boost_LIBRARIES: ${Boost_LIBRARIES}")
+INCLUDE_DIRECTORIES (${Boost_INCLUDE_DIRS})
 
 # vtk
-if (WIN32)
-    find_package(VTK 6.1 REQUIRED)
-else()
-    # prefer >=6.1, fallback to >=5.4
-    find_package(VTK 6.1 COMPONENTS vtkIOXML vtkIOGeometry vtkIOLegacy vtkIOPLY NO_MODULE)
-    IF (NOT ${VTK_FOUND})
-        find_package(VTK REQUIRED)
-    endif()
-endif()
+find_package(VTK COMPONENTS vtkIOXML vtkIOGeometry vtkIOLegacy vtkIOPLY NO_MODULE REQUIRED)
 
 message(STATUS "Found package VTK. Using version " ${VTK_VERSION})
-if("${VTK_MAJOR_VERSION}" GREATER 5)
-    set( vtk_LIBS ${VTK_LIBRARIES} )
-else()
-    set( vtk_LIBS
-        vtkCommon
-        vtkIO
-    )    
-endif()
+set( vtk_LIBS ${VTK_LIBRARIES} )
 message(STATUS "vtk libraries " ${vtk_LIBS})
 
 include(${VTK_USE_FILE})
 INCLUDE_DIRECTORIES (${VTK_INCLUDE_DIR})
 
-#set(CMAKE_CXX_FLAGS "-msse -march=native")
+if(${CMAKE_SYSTEM_PROCESSOR} STREQUAL "x86_64")
+  set(ARCH "x86_64")
+elseif(${CMAKE_SYSTEM_PROCESSOR} STREQUAL "amd64")
+  set(ARCH "x86_64")
+elseif(${CMAKE_SYSTEM_PROCESSOR} STREQUAL "AMD64")
+  set(ARCH "x86_64")
+elseif(CMAKE_SYSTEM_PROCESSOR MATCHES "^(ppc64.*|PPC64.*|powerpc64.*)")
+  set(ARCH "ppc64")
+elseif(${CMAKE_SYSTEM_PROCESSOR} MATCHES "^(aarch64|arm64)")
+  set(ARCH "aarch64")
+elseif(${CMAKE_SYSTEM_PROCESSOR} STREQUAL "unknown")
+  set(ARCH "unknown")
+  message(FATAL_ERROR "Unable to determine target architecture!  Try setting CMAKE_SYSTEM_PROCESSOR to a valid value.")
+else()
+  set(ARCH "unsupported")
+endif()
+
+if(${ARCH} STREQUAL "x86_64")
+  message(STATUS "Detected 64-bit x86 target")
+  #set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -msse -march=native")
+elseif(${ARCH} STREQUAL "ppc64")
+  message(STATUS "Detected 64-bit POWER target")
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNO_WARN_X86_INTRINSICS -DSSE_CORRECT_DENORMALS")
+elseif(${ARCH} STREQUAL "aarch64")
+  message(STATUS "Detected 64-bit ARM target")
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNO_WARN_X86_INTRINSICS -DSSE_CORRECT_DENORMALS")
+elseif(${ARCH} STREQUAL "unsupported")
+  message(FATAL_ERROR "Unsupported target architecture!  Try porting openEMS to your architecture...")
+else()
+  message(FATAL_ERROR "Unable to determine target architecture!  Aborting.")
+endif()
 
 # independent tool
 ADD_SUBDIRECTORY( nf2ff )
@@ -155,7 +182,7 @@ set(SOURCES
   openems.cpp
 )
 
-set(PUB_HEADERS openems.h)
+set(PUB_HEADERS openems.h openems_global.h)
 
 # libs
 ADD_SUBDIRECTORY( tools )
@@ -174,24 +201,35 @@ if (${MPI_CXX_FOUND})
 endif()
 
 add_library( openEMS SHARED ${SOURCES})
-set_target_properties(openEMS PROPERTIES VERSION ${LIB_VERSION_STRING} SOVERSION ${LIB_VERSION_MAJOR})
+#ADD_EXECUTABLE( openEMS main.cpp ${SOURCES})
+set_target_properties(openEMS PROPERTIES VERSION ${LIB_VERSION_STRING} SOVERSION ${LIB_VERSION_MAJOR} )
 TARGET_LINK_LIBRARIES( openEMS
   ${CSXCAD_LIBRARIES}
   ${fparser_LIBRARIES}
-  tinyxml
+  ${TinyXML_LIBRARY}
   ${HDF5_LIBRARIES}
+  ${HDF5_HL_LIBRARIES}
   ${Boost_LIBRARIES}
   ${vtk_LIBS}
   ${MPI_LIBRARIES}
 )
 
+if (WIN32)
+    # make sure "BUILD_OPENEMS_LIB" is only set for openEMS (dll) not openEMS_bin (exe)
+    target_compile_definitions(openEMS PRIVATE -DBUILD_OPENEMS_LIB )
+endif (WIN32)
+
 # main program
 ADD_EXECUTABLE( openEMS_bin main.cpp )
 SET_TARGET_PROPERTIES(openEMS_bin PROPERTIES OUTPUT_NAME openEMS)
 TARGET_LINK_LIBRARIES(openEMS_bin openEMS)
  
-INSTALL(TARGETS openEMS DESTINATION lib${LIB_SUFFIX})
-INSTALL(TARGETS openEMS_bin DESTINATION bin )
+if (WIN32)
+    INSTALL(TARGETS openEMS DESTINATION bin)
+else()
+    INSTALL(TARGETS openEMS DESTINATION lib${LIB_SUFFIX})
+endif()
+INSTALL(TARGETS openEMS_bin DESTINATION bin)
 
 if (UNIX)
     INSTALL( FILES openEMS.sh 

Common/engine_interface_base.h

@@ -22,9 +22,9 @@
 
 class Operator_Base;
 
-//! This is the abstact base for all Engine Interface classes.
+//! This is the abstract base for all Engine Interface classes.
 /*!
-	 This is the abstact base for all Engine Interface classes. It will provide unified access to the field information of the corresponding engine.
+	 This is the abstract base for all Engine Interface classes. It will provide unified access to the field information of the corresponding engine.
 	 All processing methods should only access this base class.
 */
 class Engine_Interface_Base
@@ -56,6 +56,10 @@ public:
 	virtual double* GetJField(const unsigned int* pos, double* out) const =0;
 	//! Get the total current density field by rot(H) at \p pos. \sa SetInterpolationType
 	virtual double* GetRotHField(const unsigned int* pos, double* out) const =0;
+	//! Get the (interpolated) electric flux density field at \p pos. \sa SetInterpolationType
+	virtual double* GetDField(const unsigned int* pos, double* out) const =0;
+	//! Get the (interpolated) magnetic flux density field at \p pos. \sa SetInterpolationType
+	virtual double* GetBField(const unsigned int* pos, double* out) const =0;
 
 	//! Calculate the electric field integral along a given line
 	virtual double CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const =0;

Common/operator_base.h

@@ -71,7 +71,7 @@ public:
 	//! Get the length of an FDTD edge (unit is meter).
 	virtual double GetEdgeLength(int ny, const unsigned int pos[3], bool dualMesh = false) const =0;
 
-	//! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos
+	//! Get the area around an edge for a given direction \a n and a given mesh position \a pos
 	/*!
 		This will return the area around an edge with a given direction, measured at the middle of the edge.
 		In a cartesian mesh this is equal to the NodeArea, may be different in other coordinate systems.
@@ -81,7 +81,7 @@ public:
 	//! Get the volume of an FDTD cell
 	virtual double GetCellVolume(const unsigned int pos[3], bool dualMesh = false) const =0;
 
-	//! Snap the given coodinates to mesh indices, return box dimension
+	//! Snap the given coordinates to mesh indices, return box dimension
 	virtual bool SnapToMesh(const double* coord, unsigned int* uicoord, bool dualMesh=false, bool fullMesh=false, bool* inside=NULL) const =0;
 
 	//! Snap a given box to the operator mesh, uiStart will be always <= uiStop

Common/processfields.cpp

@@ -71,6 +71,10 @@ string ProcessFields::GetFieldNameByType(DumpType type)
 		return "J-Field";
 	case ROTH_FIELD_DUMP:
 		return "RotH-Field";
+	case D_FIELD_DUMP:
+		return "D-Field";
+	case B_FIELD_DUMP:
+		return "B-Field";
 	case SAR_LOCAL_DUMP:
 		return "SAR-local";
 	case SAR_1G_DUMP:
@@ -95,6 +99,30 @@ bool ProcessFields::NeedConductivity() const
 	return false;
 }
 
+bool ProcessFields::NeedPermittivity() const
+{
+	switch (m_DumpType)
+	{
+	case D_FIELD_DUMP:
+		return true;
+	default:
+		return false;
+	}
+	return false;
+}
+
+bool ProcessFields::NeedPermeability() const
+{
+	switch (m_DumpType)
+	{
+	case B_FIELD_DUMP:
+		return true;
+	default:
+		return false;
+	}
+	return false;
+}
+
 void ProcessFields::InitProcess()
 {
 	if (Enabled==false) return;
@@ -333,6 +361,44 @@ FDTD_FLOAT**** ProcessFields::CalcField()
 			}
 		}
 		return field;
+	case D_FIELD_DUMP:
+		for (unsigned int i=0; i<numLines[0]; ++i)
+		{
+			pos[0]=posLines[0][i];
+			for (unsigned int j=0; j<numLines[1]; ++j)
+			{
+				pos[1]=posLines[1][j];
+				for (unsigned int k=0; k<numLines[2]; ++k)
+				{
+					pos[2]=posLines[2][k];
+
+					m_Eng_Interface->GetDField(pos,out);
+					field[0][i][j][k] = out[0];
+					field[1][i][j][k] = out[1];
+					field[2][i][j][k] = out[2];
+				}
+			}
+		}
+		return field;
+	case B_FIELD_DUMP:
+		for (unsigned int i=0; i<numLines[0]; ++i)
+		{
+			pos[0]=posLines[0][i];
+			for (unsigned int j=0; j<numLines[1]; ++j)
+			{
+				pos[1]=posLines[1][j];
+				for (unsigned int k=0; k<numLines[2]; ++k)
+				{
+					pos[2]=posLines[2][k];
+
+					m_Eng_Interface->GetBField(pos,out);
+					field[0][i][j][k] = out[0];
+					field[1][i][j][k] = out[1];
+					field[2][i][j][k] = out[2];
+				}
+			}
+		}
+		return field;
 	default:
 		cerr << "ProcessFields::CalcField(): Error, unknown dump type..." << endl;
 		return field;

Common/processfields.h

@@ -40,7 +40,7 @@ public:
 	  Current dump types are electric field (E_FIELD_DUMP), magnetic field (H_FIELD_DUMP),
 	  (conduction) electric current density (kappa*E) (J_FIELD_DUMP) and total current density (rotH)
 	  */
-	enum DumpType { E_FIELD_DUMP=0, H_FIELD_DUMP=1, J_FIELD_DUMP=2, ROTH_FIELD_DUMP=3, SAR_LOCAL_DUMP=20, SAR_1G_DUMP=21, SAR_10G_DUMP=22, SAR_RAW_DATA=29};
+	enum DumpType { E_FIELD_DUMP=0, H_FIELD_DUMP=1, J_FIELD_DUMP=2, ROTH_FIELD_DUMP=3, D_FIELD_DUMP=4, B_FIELD_DUMP=5, SAR_LOCAL_DUMP=20, SAR_1G_DUMP=21, SAR_10G_DUMP=22, SAR_RAW_DATA=29};
 
 	virtual std::string GetProcessingName() const {return "common field processing";}
 
@@ -60,9 +60,9 @@ public:
 
 	//! Define the Dump-Mode
 	void SetDumpMode(Engine_Interface_Base::InterpolationType mode);
-	//! This methode will dump all fields on a main cell node using 2 E-field and 4 H-fields per direction.
+	//! This method will dump all fields on a main cell node using 2 E-field and 4 H-fields per direction.
 	void SetDumpMode2Node() {SetDumpMode(Engine_Interface_Base::NODE_INTERPOLATE);}
-	//! This methode will dump all fields in the center of a main cell (dual-node) using 4 E-field and 2 H-fields per direction.
+	//! This method will dump all fields in the center of a main cell (dual-node) using 4 E-field and 2 H-fields per direction.
 	void SetDumpMode2Cell() {SetDumpMode(Engine_Interface_Base::CELL_INTERPOLATE);}
 
 	//! Set dump type: 0 for E-fields, 1 for H-fields, 2 for D-fields, 3 for B-fields, 4 for J-fields, etc...
@@ -75,6 +75,8 @@ public:
 	static std::string GetFieldNameByType(DumpType type);
 
 	virtual bool NeedConductivity() const;
+	virtual bool NeedPermittivity() const;
+	virtual bool NeedPermeability() const;
 
 protected:
 	DumpType m_DumpType;

Common/processing.cpp

@@ -203,6 +203,11 @@ void Processing::ShowSnappedCoords()
 		 << stop[0] << "," << stop[1] << "," << stop[2] << "]" << endl;
 }
 
+void Processing::SetProcessInterval(unsigned int interval)
+{
+	ProcessInterval=std::max((unsigned int)1,interval);
+}
+
 void Processing::SetProcessStartStopTime(double start, double stop)
 {
 	double dT = Op->GetTimestep();

Common/processing.h

@@ -18,6 +18,10 @@
 #ifndef PROCESSING_H
 #define PROCESSING_H
 
+#ifndef __GNUC__ // not GCC
+#include <emmintrin.h>
+#endif
+
 #include <complex>
 typedef std::complex<double> double_complex;
 #define _I double_complex(0.0,1.0)
@@ -31,6 +35,8 @@ typedef std::complex<double> double_complex;
 #include <string>
 #include <vector>
 
+#define _USE_MATH_DEFINES
+
 #include "Common/engine_interface_base.h"
 
 class Operator_Base;
@@ -59,7 +65,7 @@ public:
 
 	virtual void ShowSnappedCoords();
 
-	void SetProcessInterval(unsigned int interval) {ProcessInterval=std::max((unsigned int)1,interval);}
+	void SetProcessInterval(unsigned int interval);
 	void SetProcessStartStopTime(double start, double stop);
 
 	void AddStep(unsigned int step);
@@ -96,7 +102,7 @@ public:
 	//! Set the dump precision
 	void SetPrecision(unsigned int val) {m_precision = val;}
 
-	//! Dump probe geometry to file (will obay main or dual mesh property)
+	//! Dump probe geometry to file (will obey main or dual mesh property)
 	virtual void DumpBox2File(std::string vtkfilenameprefix) const {DumpBox2File(vtkfilenameprefix,m_dualMesh);}
 
 	//! Dump probe geometry to file
@@ -178,7 +184,7 @@ public:
 
 	void Reset();
 
-	//! Deletes all given processing's, can be helpful, but use carefull!!!
+	//! Deletes all given processing's, can be helpful, but use carefully!!!
 	void DeleteAll();
 
 	//! Invoke PreProcess() on all Processings.

Common/processmodematch.cpp

@@ -74,7 +74,7 @@ void ProcessModeMatch::InitProcess()
 
 	if (m_Eng_Interface==NULL)
 	{
-		cerr << "ProcessModeMatch::InitProcess: Error, Engine_Interface is NULL, abort mode mathcing..." << endl;
+		cerr << "ProcessModeMatch::InitProcess: Error, Engine_Interface is NULL, abort mode matching..." << endl;
 		Enabled=false;
 		return;
 	}
@@ -123,7 +123,7 @@ void ProcessModeMatch::InitProcess()
 		int res = m_ModeParser[n]->Parse(m_ModeFunction[ny], "x,y,z,rho,a,r,t");
 		if (res >= 0)
 		{
-			cerr << "ProcessModeMatch::InitProcess(): Warning, an error occured parsing the mode matching function (see below) ..." << endl;
+			cerr << "ProcessModeMatch::InitProcess(): Warning, an error occurred parsing the mode matching function (see below) ..." << endl;
 			cerr << m_ModeFunction[ny] << "\n" << string(res, ' ') << "^\n" << m_ModeParser[n]->ErrorMsg() << "\n";
 			SetEnable(false);
 			Reset();
@@ -138,7 +138,7 @@ void ProcessModeMatch::InitProcess()
 	bool dualMesh = m_ModeFieldType==1;
 	unsigned int pos[3] = {0,0,0};
 	double discLine[3] = {0,0,0};
-	double gridDelta = 1; // 1 -> mode-matching function is definied in drawing units...
+	double gridDelta = 1; // 1 -> mode-matching function is defined in drawing units...
 	double var[7];
 	pos[m_ny] = start[m_ny];
 	discLine[m_ny] = Op->GetDiscLine(m_ny,pos[m_ny],dualMesh);
@@ -174,7 +174,7 @@ void ProcessModeMatch::InitProcess()
 			for (int n=0; n<2; ++n)
 			{
 				m_ModeDist[n][posP][posPP] = m_ModeParser[n]->Eval(var); //calc mode template
-				if ((isnan(m_ModeDist[n][posP][posPP])) || (isinf(m_ModeDist[n][posP][posPP])))
+				if ((std::isnan(m_ModeDist[n][posP][posPP])) || (std::isinf(m_ModeDist[n][posP][posPP])))
 					m_ModeDist[n][posP][posPP] = 0.0;
 				norm += pow(m_ModeDist[n][posP][posPP],2) * area;
 			}

Common/processmodematch.h

@@ -24,7 +24,7 @@ class CSFunctionParser;
 
 //! Processing class to match a mode to a given analytic function and return the integral value.
 /*!
-  The analytric function has to be definied in drawing units.
+  The analytric function has to be defined in drawing units.
   It will return the integral value and the mode purity as a secondary value.
   */
 class ProcessModeMatch : public ProcessIntegral

Common/processvoltage.cpp

@@ -16,6 +16,7 @@
 */
 
 #include "processvoltage.h"
+#include "FDTD/engine_interface_fdtd.h"
 #include <iomanip>
 
 ProcessVoltage::ProcessVoltage(Engine_Interface_Base* eng_if) : ProcessIntegral(eng_if)
@@ -33,6 +34,12 @@ std::string ProcessVoltage::GetIntegralName(int row) const
 	return "unknown";
 }
 
+void ProcessVoltage::DefineStartStopCoord(double* dstart, double* dstop)
+{
+	Op->SnapToMesh(dstart, start, m_dualMesh, false, m_start_inside);
+	Op->SnapToMesh(dstop, stop, m_dualMesh, false, m_stop_inside);
+}
+
 double ProcessVoltage::CalcIntegral()
 {
 	//integrate voltages from start to stop on a line

Common/processvoltage.h

@@ -31,6 +31,8 @@ public:
 
 	virtual std::string GetIntegralName(int row) const;
 
+	virtual void DefineStartStopCoord(double* dstart, double* dstop);
+
 	virtual double CalcIntegral();
 
 protected:

FDTD/engine.cpp

@@ -123,16 +123,37 @@ void Engine::UpdateVoltages(unsigned int startX, unsigned int numX)
 				shift[2]=pos[2];
 				//do the updates here
 				//for x
-				volt[0][pos[0]][pos[1]][pos[2]] *= Op->vv[0][pos[0]][pos[1]][pos[2]];
-				volt[0][pos[0]][pos[1]][pos[2]] += Op->vi[0][pos[0]][pos[1]][pos[2]] * ( curr[2][pos[0]][pos[1]][pos[2]] - curr[2][pos[0]][pos[1]-shift[1]][pos[2]] - curr[1][pos[0]][pos[1]][pos[2]] + curr[1][pos[0]][pos[1]][pos[2]-shift[2]]);
+				volt[0][pos[0]][pos[1]][pos[2]] *=
+				    Op->vv[0][pos[0]][pos[1]][pos[2]];
+				volt[0][pos[0]][pos[1]][pos[2]] +=
+				    Op->vi[0][pos[0]][pos[1]][pos[2]] * (
+				        curr[2][pos[0]][pos[1]         ][pos[2]         ] -
+				        curr[2][pos[0]][pos[1]-shift[1]][pos[2]         ] -
+				        curr[1][pos[0]][pos[1]         ][pos[2]         ] +
+				        curr[1][pos[0]][pos[1]         ][pos[2]-shift[2]]
+				    );
 
 				//for y
-				volt[1][pos[0]][pos[1]][pos[2]] *= Op->vv[1][pos[0]][pos[1]][pos[2]];
-				volt[1][pos[0]][pos[1]][pos[2]] += Op->vi[1][pos[0]][pos[1]][pos[2]] * ( curr[0][pos[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]-shift[2]] - curr[2][pos[0]][pos[1]][pos[2]] + curr[2][pos[0]-shift[0]][pos[1]][pos[2]]);
+				volt[1][pos[0]][pos[1]][pos[2]] *=
+				    Op->vv[1][pos[0]][pos[1]][pos[2]];
+				volt[1][pos[0]][pos[1]][pos[2]] +=
+				    Op->vi[1][pos[0]][pos[1]][pos[2]] * (
+				        curr[0][pos[0]         ][pos[1]][pos[2]         ] -
+				        curr[0][pos[0]         ][pos[1]][pos[2]-shift[2]] -
+				        curr[2][pos[0]         ][pos[1]][pos[2]         ] +
+				        curr[2][pos[0]-shift[0]][pos[1]][pos[2]         ]
+				    );
 
 				//for z
-				volt[2][pos[0]][pos[1]][pos[2]] *= Op->vv[2][pos[0]][pos[1]][pos[2]];
-				volt[2][pos[0]][pos[1]][pos[2]] += Op->vi[2][pos[0]][pos[1]][pos[2]] * ( curr[1][pos[0]][pos[1]][pos[2]] - curr[1][pos[0]-shift[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]] + curr[0][pos[0]][pos[1]-shift[1]][pos[2]]);
+				volt[2][pos[0]][pos[1]][pos[2]] *=
+				    Op->vv[2][pos[0]][pos[1]][pos[2]];
+				volt[2][pos[0]][pos[1]][pos[2]] +=
+				    Op->vi[2][pos[0]][pos[1]][pos[2]] * (
+				        curr[1][pos[0]         ][pos[1]         ][pos[2]] -
+				        curr[1][pos[0]-shift[0]][pos[1]         ][pos[2]] -
+				        curr[0][pos[0]         ][pos[1]         ][pos[2]] +
+				        curr[0][pos[0]         ][pos[1]-shift[1]][pos[2]]
+				    );
 			}
 		}
 		++pos[0];
@@ -151,16 +172,37 @@ void Engine::UpdateCurrents(unsigned int startX, unsigned int numX)
 			{
 				//do the updates here
 				//for x
-				curr[0][pos[0]][pos[1]][pos[2]] *= Op->ii[0][pos[0]][pos[1]][pos[2]];
-				curr[0][pos[0]][pos[1]][pos[2]] += Op->iv[0][pos[0]][pos[1]][pos[2]] * ( volt[2][pos[0]][pos[1]][pos[2]] - volt[2][pos[0]][pos[1]+1][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]] + volt[1][pos[0]][pos[1]][pos[2]+1]);
+				curr[0][pos[0]][pos[1]][pos[2]] *=
+				    Op->ii[0][pos[0]][pos[1]][pos[2]];
+				curr[0][pos[0]][pos[1]][pos[2]] +=
+				    Op->iv[0][pos[0]][pos[1]][pos[2]] * (
+				        volt[2][pos[0]][pos[1]  ][pos[2]  ] -
+				        volt[2][pos[0]][pos[1]+1][pos[2]  ] -
+				        volt[1][pos[0]][pos[1]  ][pos[2]  ] +
+				        volt[1][pos[0]][pos[1]  ][pos[2]+1]
+				    );
 
 				//for y
-				curr[1][pos[0]][pos[1]][pos[2]] *= Op->ii[1][pos[0]][pos[1]][pos[2]];
-				curr[1][pos[0]][pos[1]][pos[2]] += Op->iv[1][pos[0]][pos[1]][pos[2]] * ( volt[0][pos[0]][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]+1] - volt[2][pos[0]][pos[1]][pos[2]] + volt[2][pos[0]+1][pos[1]][pos[2]]);
+				curr[1][pos[0]][pos[1]][pos[2]] *=
+				    Op->ii[1][pos[0]][pos[1]][pos[2]];
+				curr[1][pos[0]][pos[1]][pos[2]] +=
+				    Op->iv[1][pos[0]][pos[1]][pos[2]] * (
+				        volt[0][pos[0]  ][pos[1]][pos[2]  ] -
+				        volt[0][pos[0]  ][pos[1]][pos[2]+1] -
+				        volt[2][pos[0]  ][pos[1]][pos[2]  ] +
+				        volt[2][pos[0]+1][pos[1]][pos[2]  ]
+				    );
 
 				//for z
-				curr[2][pos[0]][pos[1]][pos[2]] *= Op->ii[2][pos[0]][pos[1]][pos[2]];
-				curr[2][pos[0]][pos[1]][pos[2]] += Op->iv[2][pos[0]][pos[1]][pos[2]] * ( volt[1][pos[0]][pos[1]][pos[2]] - volt[1][pos[0]+1][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]] + volt[0][pos[0]][pos[1]+1][pos[2]]);
+				curr[2][pos[0]][pos[1]][pos[2]] *=
+				    Op->ii[2][pos[0]][pos[1]][pos[2]];
+				curr[2][pos[0]][pos[1]][pos[2]] +=
+				    Op->iv[2][pos[0]][pos[1]][pos[2]] * (
+				        volt[1][pos[0]  ][pos[1]  ][pos[2]] -
+				        volt[1][pos[0]+1][pos[1]  ][pos[2]] -
+				        volt[0][pos[0]  ][pos[1]  ][pos[2]] +
+				        volt[0][pos[0]  ][pos[1]+1][pos[2]]
+				    );
 			}
 		}
 		++pos[0];

FDTD/engine.h

@@ -47,6 +47,8 @@ public:
 
 	virtual unsigned int GetNumberOfTimesteps() {return numTS;}
 
+	virtual void NextInterval(float curr_speed) {};
+
 	//this access functions muss be overloaded by any new engine using a different storage model
 	inline virtual FDTD_FLOAT GetVolt( unsigned int n, unsigned int x, unsigned int y, unsigned int z )		const { return volt[n][x][y][z]; }
 	inline virtual FDTD_FLOAT GetVolt( unsigned int n, const unsigned int pos[3] )							const { return volt[n][pos[0]][pos[1]][pos[2]]; }

FDTD/engine_interface_fdtd.cpp

@@ -47,6 +47,11 @@ double* Engine_Interface_FDTD::GetJField(const unsigned int* pos, double* out) c
 	return GetRawInterpolatedField(pos, out, 1);
 }
 
+double* Engine_Interface_FDTD::GetDField(const unsigned int* pos, double* out) const
+{
+	return GetRawInterpolatedField(pos, out, 3);
+}
+
 double* Engine_Interface_FDTD::GetRotHField(const unsigned int* pos, double* out) const
 {
 	return GetRawInterpolatedField(pos, out, 2);
@@ -116,6 +121,27 @@ double* Engine_Interface_FDTD::GetRawInterpolatedField(const unsigned int* pos,
 }
 
 double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) const
+{
+		return GetRawInterpolatedDualField(pos, out, 0);
+}
+
+double* Engine_Interface_FDTD::GetBField(const unsigned int* pos, double* out) const
+{
+		return GetRawInterpolatedDualField(pos, out, 1);
+}
+
+double Engine_Interface_FDTD::GetRawDualField(unsigned int n, const unsigned int* pos, int type) const
+{
+	double value = m_Eng->GetCurr(n,pos[0],pos[1],pos[2]);
+	double delta = m_Op->GetEdgeLength(n,pos,true);
+	if ((type==0) && (delta))
+		return value/delta;
+	if ((type==1) && (m_Op->m_mueR) && (delta))
+		return value*m_Op->m_mueR[n][pos[0]][pos[1]][pos[2]]/delta;
+	return 0.0;
+}
+
+double* Engine_Interface_FDTD::GetRawInterpolatedDualField(const unsigned int* pos, double* out, int type) const
 {
 	unsigned int iPos[] = {pos[0],pos[1],pos[2]};
 	int nP,nPP;
@@ -124,9 +150,9 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
 	{
 	default:
 	case NO_INTERPOLATION:
-		out[0] = m_Eng->GetCurr(0,pos) / m_Op->GetEdgeLength(0,pos,true);
-		out[1] = m_Eng->GetCurr(1,pos) / m_Op->GetEdgeLength(1,pos,true);
-		out[2] = m_Eng->GetCurr(2,pos) / m_Op->GetEdgeLength(2,pos,true);
+		out[0] = GetRawDualField(0, pos, type);
+		out[1] = GetRawDualField(1, pos, type);
+		out[2] = GetRawDualField(2, pos, type);
 		break;
 	case NODE_INTERPOLATE:
 		for (int n=0; n<3; ++n)
@@ -138,13 +164,13 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
 				out[n] = 0;
 				continue;
 			}
-			out[n]=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true);
+			out[n] = GetRawDualField(n, iPos, type);
 			--iPos[nP];
-			out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true);
+			out[n]+= GetRawDualField(n, iPos, type);
 			--iPos[nPP];
-			out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true);
+			out[n]+= GetRawDualField(n, iPos, type);
 			++iPos[nP];
-			out[n]+=m_Eng->GetCurr(n,iPos)/m_Op->GetEdgeLength(n,iPos,true);
+			out[n]+= GetRawDualField(n, iPos, type);
 			++iPos[nPP];
 			out[n]/=4;
 		}
@@ -153,7 +179,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
 		for (int n=0; n<3; ++n)
 		{
 			delta = m_Op->GetEdgeLength(n,iPos,true);
-			out[n] = m_Eng->GetCurr(n,iPos);
+			out[n] = GetRawDualField(n, iPos, type);
 			if ((pos[n]>=m_Op->GetNumberOfLines(n,true)-1))
 			{
 				out[n] = 0; //magnetic field on the outer boundaries is always zero
@@ -162,7 +188,7 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
 			++iPos[n];
 			double deltaUp = m_Op->GetEdgeLength(n,iPos,true);
 			double deltaRel = delta / (delta+deltaUp);
-			out[n] = out[n]*(1.0-deltaRel)/delta + (double)m_Eng->GetCurr(n,iPos)/deltaUp*deltaRel;
+			out[n] = out[n]*(1.0-deltaRel) + (double)GetRawDualField(n, iPos, type)*deltaRel;
 			--iPos[n];
 		}
 		break;
@@ -173,6 +199,12 @@ double* Engine_Interface_FDTD::GetHField(const unsigned int* pos, double* out) c
 
 double Engine_Interface_FDTD::CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const
 {
+	if (((start[0]!=stop[0]) + (start[1]!=stop[1]) + (start[2]!=stop[2]))!=1)
+	{
+		cerr << "Engine_Interface_FDTD::CalcVoltageIntegral: Error, only a 1D/line integration is allowed" << endl;
+		return 0;
+	}
+	//cerr << "CalcVoltageIntegral" << start[0] << ", " << start[1] << ", " << start[2] << " -> " << stop[0] << ", " << stop[1] << ", " << stop[2] << ", " << endl;
 	double result=0;
 	for (int n=0; n<3; ++n)
 	{
@@ -181,6 +213,7 @@ double Engine_Interface_FDTD::CalcVoltageIntegral(const unsigned int* start, con
 			unsigned int pos[3]={start[0],start[1],start[2]};
 			for (; pos[n]<stop[n]; ++pos[n])
 				result += m_Eng->GetVolt(n,pos[0],pos[1],pos[2]);
+
 		}
 		else
 		{
@@ -192,6 +225,34 @@ double Engine_Interface_FDTD::CalcVoltageIntegral(const unsigned int* start, con
 	return result;
 }
 
+//double Engine_Interface_FDTD::CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const
+//{
+//	//cerr << "CalcVoltageIntegral" << start[0] << ", " << start[1] << ", " << start[2] << " -> " << stop[0] << ", " << stop[1] << ", " << stop[2] << ", " << endl;
+//	double result=0;
+//	//unsigned int pos[3]={min(start[0],stop[0]),min(start[1],stop[1]),min(start[2],stop[2])};
+//	unsigned int pos[3]={start[0],start[1],start[2]};
+//	for (int n=0; n<3; ++n)
+//	{
+//		if (start[n]<stop[n])
+//		{
+//			for (; pos[n]<stop[n]; ++pos[n])
+//			{
+//				//cerr << "at pos: " << n << ": " << pos[0] << ", " << pos[1] << ", " << pos[2] << endl;
+//				result += m_Eng->GetVolt(n,pos[0],pos[1],pos[2]);
+//			}
+//		}
+//		else
+//		{
+//			for (--pos[n]; pos[n]>=stop[n]; --pos[n])
+//			{
+//				//cerr << "at neg: " << n << ": " << pos[0] << ", " << pos[1] << ", " << pos[2] << endl;
+//				result -= m_Eng->GetVolt(n,pos[0],pos[1],pos[2]);
+//			}
+//		}
+//		pos[n] = stop[n];
+//	}
+//	return result;
+//}
 
 double Engine_Interface_FDTD::GetRawField(unsigned int n, const unsigned int* pos, int type) const
 {
@@ -201,6 +262,8 @@ double Engine_Interface_FDTD::GetRawField(unsigned int n, const unsigned int* po
 		return value/delta;
 	if ((type==1) && (m_Op->m_kappa) && (delta))
 		return value*m_Op->m_kappa[n][pos[0]][pos[1]][pos[2]]/delta;
+	if ((type==3) && (m_Op->m_epsR) && (delta))
+		return value*m_Op->m_epsR[n][pos[0]][pos[1]][pos[2]]/delta;
 	if (type==2) //calc rot(H)
 	{
 		int nP = (n+1)%3;

FDTD/engine_interface_fdtd.h

@@ -44,6 +44,8 @@ public:
 	virtual double* GetHField(const unsigned int* pos, double* out) const;
 	virtual double* GetJField(const unsigned int* pos, double* out) const;
 	virtual double* GetRotHField(const unsigned int* pos, double* out) const;
+	virtual double* GetDField(const unsigned int* pos, double* out) const;
+	virtual double* GetBField(const unsigned int* pos, double* out) const;
 
 	virtual double CalcVoltageIntegral(const unsigned int* start, const unsigned int* stop) const;
 
@@ -56,10 +58,15 @@ protected:
 	Operator* m_Op;
 	Engine* m_Eng;
 
-	//! Internal method to get an interpolated field of a given type. (0: E, 1: J, 2: rotH)
+	//! Internal method to get an interpolated field of a given type. (0: E, 1: J, 2: rotH, 3: D)
 	virtual double* GetRawInterpolatedField(const unsigned int* pos, double* out, int type) const;
-	//! Internal method to get a raw field of a given type. (0: E, 1: J, 2: rotH)
+	//! Internal method to get a raw field of a given type. (0: E, 1: J, 2: rotH, 3: D)
 	virtual double GetRawField(unsigned int n, const unsigned int* pos, int type) const;
+
+	//! Internal method to get an interpolated dual field of a given type. (0: H, 1: B)
+	virtual double* GetRawInterpolatedDualField(const unsigned int* pos, double* out, int type) const;
+	//! Internal method to get a raw dual field of a given type. (0: H, 1: B)
+	virtual double GetRawDualField(unsigned int n, const unsigned int* pos, int type) const;
 };
 
 #endif // ENGINE_INTERFACE_FDTD_H

FDTD/engine_multithread.cpp

@@ -32,7 +32,10 @@
 #include "boost/date_time/posix_time/posix_time.hpp"
 #include "boost/date_time/gregorian/gregorian.hpp"
 #include <iomanip>
+
+#ifndef SSE_CORRECT_DENORMALS
 #include <xmmintrin.h>
+#endif
 
 //! \brief construct an Engine_Multithread instance
 //! it's the responsibility of the caller to free the returned pointer
@@ -52,6 +55,12 @@ Engine_Multithread::Engine_Multithread(const Operator_Multithread* op) : ENGINE_
 	m_IterateBarrier = 0;
 	m_startBarrier = 0;
 	m_stopBarrier = 0;
+	m_thread_group = 0;
+	m_max_numThreads = boost::thread::hardware_concurrency();
+	m_numThreads = 0;
+	m_last_speed = 0;
+	m_opt_speed = false;
+	m_stopThreads = true;
 
 #ifdef ENABLE_DEBUG_TIME
 	m_MPI_Barrier = 0;
@@ -89,27 +98,85 @@ void Engine_Multithread::setNumThreads( unsigned int numThreads )
 void Engine_Multithread::Init()
 {
 	m_stopThreads = true;
+	m_opt_speed = false;
 	ENGINE_MULTITHREAD_BASE::Init();
 
 	// initialize threads
 	m_stopThreads = false;
 	if (m_numThreads == 0)
-		m_numThreads = boost::thread::hardware_concurrency();
+	{
+		m_opt_speed = true;
+		m_numThreads = 1;
+	}
+	else if (m_numThreads > m_max_numThreads)
+		m_numThreads = m_max_numThreads;
+
+#ifdef MPI_SUPPORT
+	m_MPI_Barrier = 0;
+#endif
+	this->changeNumThreads(m_numThreads);
+}
+
+void Engine_Multithread::Reset()
+{
+	if (m_thread_group!=0) // prevent multiple invocations
+	{
+		ClearExtensions(); //prevent extensions from interfering with thread reset...
+
+		// stop the threads
+		//NS_Engine_Multithread::DBG().cout() << "stopping all threads" << endl;
+		m_thread_group->interrupt_all();
+		m_thread_group->join_all(); // wait for termination
+		delete m_IterateBarrier;
+		m_IterateBarrier = 0;
+		delete m_startBarrier;
+		m_startBarrier = 0;
+		delete m_stopBarrier;
+		m_stopBarrier = 0;
+		delete m_thread_group;
+		m_thread_group = 0;
+	}
+
+	ENGINE_MULTITHREAD_BASE::Reset();
+}
+
+void Engine_Multithread::changeNumThreads(unsigned int numThreads)
+{
+	if (m_thread_group!=0)
+	{
+		m_thread_group->interrupt_all();
+		m_thread_group->join_all(); // wait for termination
+		delete m_thread_group;
+		m_thread_group = 0;
+		//m_stopThreads = false;
+	}
+
+	m_numThreads = numThreads;
+
+	if (g_settings.GetVerboseLevel()>0)
+		cout << "Multithreaded engine using " << m_numThreads << " threads. Utilization: (";
 
 	vector<unsigned int> m_Start_Lines;
 	vector<unsigned int> m_Stop_Lines;
 	m_Op_MT->CalcStartStopLines( m_numThreads, m_Start_Lines, m_Stop_Lines );
 
-	if (g_settings.GetVerboseLevel()>0)
-		cout << "Multithreaded engine using " << m_numThreads << " threads. Utilization: (";
+	if (m_IterateBarrier!=0)
+		delete m_IterateBarrier;
+	// make sure all threads are waiting
 	m_IterateBarrier = new boost::barrier(m_numThreads); // numThread workers
 
+	if (m_startBarrier!=0)
+		delete m_startBarrier;
 	m_startBarrier = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller
+
+	if (m_stopBarrier!=0)
+		delete m_stopBarrier;
 	m_stopBarrier = new boost::barrier(m_numThreads+1); // numThread workers + 1 controller
 #ifdef MPI_SUPPORT
 	m_MPI_Barrier = 0;
 #endif
 
+	m_thread_group = new boost::thread_group();
 	for (unsigned int n=0; n<m_numThreads; n++)
 	{
 		unsigned int start = m_Start_Lines.at(n);
@@ -127,50 +194,44 @@ void Engine_Multithread::Init()
 				cout << stop-start+1 << ";";
 //		NS_Engine_Multithread::DBG().cout() << "###DEBUG## Thread " << n << ": start=" << start << " stop=" << stop  << " stop_h=" << stop_h << std::endl;
 		boost::thread *t = new boost::thread( NS_Engine_Multithread::thread(this,start,stop,stop_h,n) );
-		m_thread_group.add_thread( t );
+		m_thread_group->add_thread( t );
 	}
 
 	for (size_t n=0; n<m_Eng_exts.size(); ++n)
 		m_Eng_exts.at(n)->SetNumberOfThreads(m_numThreads);
 }
 
-void Engine_Multithread::Reset()
-{
-	if (!m_stopThreads) // prevent multiple invocations
-	{
-		ClearExtensions(); //prevent extensions from interfering with thread reset...
-
-		// stop the threads
-		//NS_Engine_Multithread::DBG().cout() << "stopping all threads" << endl;
-		m_iterTS = 1;
-		m_startBarrier->wait(); // start the threads
-		m_stopThreads = true;
-		m_stopBarrier->wait(); // wait for the threads to finish
-		m_thread_group.join_all(); // wait for termination
-		delete m_IterateBarrier;
-		m_IterateBarrier = 0;
-		delete m_startBarrier;
-		m_startBarrier = 0;
-		delete m_stopBarrier;
-		m_stopBarrier = 0;
-	}
-
-	ENGINE_MULTITHREAD_BASE::Reset();
-}
-
 bool Engine_Multithread::IterateTS(unsigned int iterTS)
 {
 	m_iterTS = iterTS;
 
-	//cout << "bool Engine_Multithread::IterateTS(): starting threads ...";
+	//cerr << "bool Engine_Multithread::IterateTS(): starting threads ...";
 	m_startBarrier->wait(); // start the threads
 
-	//cout << "... threads started";
+	//cerr << "... threads started" << endl;
 
 	m_stopBarrier->wait(); // wait for the threads to finish <iterTS> time steps
+
 	return true;
 }
 
+void Engine_Multithread::NextInterval(float curr_speed)
+{
+	ENGINE_MULTITHREAD_BASE::NextInterval(curr_speed);
+	if (!m_opt_speed) return;
+	if (curr_speed<m_last_speed)
+	{
+		this->changeNumThreads(m_numThreads-1);
+		cout << "Multithreaded Engine: Best performance found using " << m_numThreads << " threads." << std::endl;
+		m_opt_speed = false;
+	}
+	else if (m_numThreads<m_max_numThreads)
+	{
+		m_last_speed = curr_speed;
+		this->changeNumThreads(m_numThreads+1);
+	}
+}
+
 void Engine_Multithread::DoPreVoltageUpdates(int threadID)
 {
 	//execute extensions in reverse order -> highest priority gets access to the voltages last
@@ -266,6 +327,12 @@ void thread::operator()()
 		m_enginePtr->m_startBarrier->wait();
 		//cout << "Thread " << boost::this_thread::get_id() << " waiting... started." << endl;
 
+		if (m_enginePtr->m_stopThreads)
+		{
+			//DBG().cout() << "Thread " << m_threadID << " (" << boost::this_thread::get_id() << ") stop!." << endl;
+			return;
+		}
+
 		DEBUG_TIME( Timer timer1 );
 
 		for (unsigned int iter=0; iter<m_enginePtr->m_iterTS; ++iter)

FDTD/engine_multithread.h

@@ -26,11 +26,13 @@
 #include <boost/fusion/container/list/list_fwd.hpp>
 #include <boost/fusion/include/list_fwd.hpp>
 
-//#ifdef WIN32
-//#include <Winsock2.h> // for struct timeval
-//#endif
-
+#include "tools/useful.h"
+#ifndef __GNUC__
+#include <Winsock2.h> // for struct timeval
+#else
 #include <sys/time.h>
+#endif
+
 
 #ifdef MPI_SUPPORT
 	#define ENGINE_MULTITHREAD_BASE Engine_MPI
@@ -88,6 +90,7 @@ public:
 	virtual void setNumThreads( unsigned int numThreads );
 	virtual void Init();
 	virtual void Reset();
+	virtual void NextInterval(float curr_speed);
 
 	//! Iterate \a iterTS number of timesteps
 	virtual bool IterateTS(unsigned int iterTS);
@@ -102,13 +105,17 @@ public:
 
 protected:
 	Engine_Multithread(const Operator_Multithread* op);
+	void changeNumThreads(unsigned int numThreads);
 	const Operator_Multithread* m_Op_MT;
-	boost::thread_group m_thread_group;
+	boost::thread_group *m_thread_group;
 	boost::barrier *m_startBarrier, *m_stopBarrier;
 	boost::barrier *m_IterateBarrier;
 	volatile unsigned int m_iterTS;
 	unsigned int m_numThreads; //!< number of worker threads
+	unsigned int m_max_numThreads; //!< max. number of worker threads
 	volatile bool m_stopThreads;
+	bool m_opt_speed;
+	float m_last_speed;
 
 #ifdef MPI_SUPPORT
 	/*! Workaround needed for subgridding scheme... (see Engine_CylinderMultiGrid)

FDTD/engine_sse.cpp

@@ -15,7 +15,10 @@
 *	along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#ifndef SSE_CORRECT_DENORMALS
 #include <xmmintrin.h>
+#endif
+
 #include "engine_sse.h"
 
 //! \brief construct an Engine_sse instance
@@ -88,16 +91,37 @@ void Engine_sse::UpdateVoltages(unsigned int startX, unsigned int numX)
 			for (pos[2]=1; pos[2]<numVectors; ++pos[2])
 			{
 				// x-polarization
-				f4_volt[0][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv[0][pos[0]][pos[1]][pos[2]].v;
-				f4_volt[0][pos[0]][pos[1]][pos[2]].v += Op->f4_vi[0][pos[0]][pos[1]][pos[2]].v * ( f4_curr[2][pos[0]][pos[1]][pos[2]].v - f4_curr[2][pos[0]][pos[1]-shift[1]][pos[2]].v - f4_curr[1][pos[0]][pos[1]][pos[2]].v + f4_curr[1][pos[0]][pos[1]][pos[2]-1].v );
+				f4_volt[0][pos[0]][pos[1]][pos[2]].v *=
+				    Op->f4_vv[0][pos[0]][pos[1]][pos[2]].v;
+				f4_volt[0][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_vi[0][pos[0]][pos[1]][pos[2]].v * (
+				        f4_curr[2][pos[0]][pos[1]         ][pos[2]].v -
+				        f4_curr[2][pos[0]][pos[1]-shift[1]][pos[2]].v -
+				        f4_curr[1][pos[0]][pos[1]         ][pos[2]].v +
+				        f4_curr[1][pos[0]][pos[1]         ][pos[2]-1].v
+				    );
 
 				// y-polarization
-				f4_volt[1][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv[1][pos[0]][pos[1]][pos[2]].v;
-				f4_volt[1][pos[0]][pos[1]][pos[2]].v += Op->f4_vi[1][pos[0]][pos[1]][pos[2]].v * ( f4_curr[0][pos[0]][pos[1]][pos[2]].v - f4_curr[0][pos[0]][pos[1]][pos[2]-1].v - f4_curr[2][pos[0]][pos[1]][pos[2]].v + f4_curr[2][pos[0]-shift[0]][pos[1]][pos[2]].v);
+				f4_volt[1][pos[0]][pos[1]][pos[2]].v *=
+				    Op->f4_vv[1][pos[0]][pos[1]][pos[2]].v;
+				f4_volt[1][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_vi[1][pos[0]][pos[1]][pos[2]].v * (
+				        f4_curr[0][pos[0]         ][pos[1]][pos[2]  ].v -
+				        f4_curr[0][pos[0]         ][pos[1]][pos[2]-1].v -
+				        f4_curr[2][pos[0]         ][pos[1]][pos[2]  ].v +
+				        f4_curr[2][pos[0]-shift[0]][pos[1]][pos[2]  ].v
+				    );
 
 				// z-polarization
-				f4_volt[2][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv[2][pos[0]][pos[1]][pos[2]].v;
-				f4_volt[2][pos[0]][pos[1]][pos[2]].v += Op->f4_vi[2][pos[0]][pos[1]][pos[2]].v * ( f4_curr[1][pos[0]][pos[1]][pos[2]].v - f4_curr[1][pos[0]-shift[0]][pos[1]][pos[2]].v - f4_curr[0][pos[0]][pos[1]][pos[2]].v + f4_curr[0][pos[0]][pos[1]-shift[1]][pos[2]].v);
+				f4_volt[2][pos[0]][pos[1]][pos[2]].v *=
+				    Op->f4_vv[2][pos[0]][pos[1]][pos[2]].v;
+				f4_volt[2][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_vi[2][pos[0]][pos[1]][pos[2]].v * (
+				        f4_curr[1][pos[0]         ][pos[1]         ][pos[2]].v -
+				        f4_curr[1][pos[0]-shift[0]][pos[1]         ][pos[2]].v -
+				        f4_curr[0][pos[0]         ][pos[1]         ][pos[2]].v +
+				        f4_curr[0][pos[0]         ][pos[1]-shift[1]][pos[2]].v
+				    );
 			}
 
 			// for pos[2] = 0
@@ -106,20 +130,41 @@ void Engine_sse::UpdateVoltages(unsigned int startX, unsigned int numX)
 			temp.f[1] = f4_curr[1][pos[0]][pos[1]][numVectors-1].f[0];
 			temp.f[2] = f4_curr[1][pos[0]][pos[1]][numVectors-1].f[1];
 			temp.f[3] = f4_curr[1][pos[0]][pos[1]][numVectors-1].f[2];
-			f4_volt[0][pos[0]][pos[1]][0].v *= Op->f4_vv[0][pos[0]][pos[1]][0].v;
-			f4_volt[0][pos[0]][pos[1]][0].v += Op->f4_vi[0][pos[0]][pos[1]][0].v * ( f4_curr[2][pos[0]][pos[1]][0].v - f4_curr[2][pos[0]][pos[1]-shift[1]][0].v - f4_curr[1][pos[0]][pos[1]][0].v + temp.v );
+			f4_volt[0][pos[0]][pos[1]][0].v *=
+			    Op->f4_vv[0][pos[0]][pos[1]][0].v;
+			f4_volt[0][pos[0]][pos[1]][0].v +=
+			    Op->f4_vi[0][pos[0]][pos[1]][0].v * (
+			        f4_curr[2][pos[0]][pos[1]         ][0].v -
+			        f4_curr[2][pos[0]][pos[1]-shift[1]][0].v -
+			        f4_curr[1][pos[0]][pos[1]         ][0].v +
+			        temp.v
+			    );
 
 			// y-polarization
 			temp.f[0] = 0;
 			temp.f[1] = f4_curr[0][pos[0]][pos[1]][numVectors-1].f[0];
 			temp.f[2] = f4_curr[0][pos[0]][pos[1]][numVectors-1].f[1];
 			temp.f[3] = f4_curr[0][pos[0]][pos[1]][numVectors-1].f[2];
-			f4_volt[1][pos[0]][pos[1]][0].v *= Op->f4_vv[1][pos[0]][pos[1]][0].v;
-			f4_volt[1][pos[0]][pos[1]][0].v += Op->f4_vi[1][pos[0]][pos[1]][0].v * ( f4_curr[0][pos[0]][pos[1]][0].v - temp.v - f4_curr[2][pos[0]][pos[1]][0].v + f4_curr[2][pos[0]-shift[0]][pos[1]][0].v);
+			f4_volt[1][pos[0]][pos[1]][0].v *=
+			    Op->f4_vv[1][pos[0]][pos[1]][0].v;
+			f4_volt[1][pos[0]][pos[1]][0].v +=
+			    Op->f4_vi[1][pos[0]][pos[1]][0].v * (
+			        f4_curr[0][pos[0]         ][pos[1]][0].v -
+			        temp.v -
+			        f4_curr[2][pos[0]         ][pos[1]][0].v +
+			        f4_curr[2][pos[0]-shift[0]][pos[1]][0].v
+			    );
 
 			// z-polarization
-			f4_volt[2][pos[0]][pos[1]][0].v *= Op->f4_vv[2][pos[0]][pos[1]][0].v;
-			f4_volt[2][pos[0]][pos[1]][0].v += Op->f4_vi[2][pos[0]][pos[1]][0].v * ( f4_curr[1][pos[0]][pos[1]][0].v - f4_curr[1][pos[0]-shift[0]][pos[1]][0].v - f4_curr[0][pos[0]][pos[1]][0].v + f4_curr[0][pos[0]][pos[1]-shift[1]][0].v);
+			f4_volt[2][pos[0]][pos[1]][0].v *=
+			    Op->f4_vv[2][pos[0]][pos[1]][0].v;
+			f4_volt[2][pos[0]][pos[1]][0].v +=
+			    Op->f4_vi[2][pos[0]][pos[1]][0].v * (
+			        f4_curr[1][pos[0]         ][pos[1]         ][0].v -
+			        f4_curr[1][pos[0]-shift[0]][pos[1]         ][0].v -
+			        f4_curr[0][pos[0]         ][pos[1]         ][0].v +
+			        f4_curr[0][pos[0]         ][pos[1]-shift[1]][0].v
+			    );
 		}
 		++pos[0];
 	}
@@ -138,16 +183,37 @@ void Engine_sse::UpdateCurrents(unsigned int startX, unsigned int numX)
 			for (pos[2]=0; pos[2]<numVectors-1; ++pos[2])
 			{
 				// x-pol
-				f4_curr[0][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii[0][pos[0]][pos[1]][pos[2]].v;
-				f4_curr[0][pos[0]][pos[1]][pos[2]].v += Op->f4_iv[0][pos[0]][pos[1]][pos[2]].v * ( f4_volt[2][pos[0]][pos[1]][pos[2]].v - f4_volt[2][pos[0]][pos[1]+1][pos[2]].v - f4_volt[1][pos[0]][pos[1]][pos[2]].v + f4_volt[1][pos[0]][pos[1]][pos[2]+1].v);
+				f4_curr[0][pos[0]][pos[1]][pos[2]].v *=
+				    Op->f4_ii[0][pos[0]][pos[1]][pos[2]].v;
+				f4_curr[0][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_iv[0][pos[0]][pos[1]][pos[2]].v * (
+				        f4_volt[2][pos[0]][pos[1]  ][pos[2]  ].v -
+				        f4_volt[2][pos[0]][pos[1]+1][pos[2]  ].v -
+				        f4_volt[1][pos[0]][pos[1]  ][pos[2]  ].v +
+				        f4_volt[1][pos[0]][pos[1]  ][pos[2]+1].v
+				    );
 
 				// y-pol
-				f4_curr[1][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii[1][pos[0]][pos[1]][pos[2]].v;
-				f4_curr[1][pos[0]][pos[1]][pos[2]].v += Op->f4_iv[1][pos[0]][pos[1]][pos[2]].v * ( f4_volt[0][pos[0]][pos[1]][pos[2]].v - f4_volt[0][pos[0]][pos[1]][pos[2]+1].v - f4_volt[2][pos[0]][pos[1]][pos[2]].v + f4_volt[2][pos[0]+1][pos[1]][pos[2]].v);
+				f4_curr[1][pos[0]][pos[1]][pos[2]].v *=
+				    Op->f4_ii[1][pos[0]][pos[1]][pos[2]].v;
+				f4_curr[1][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_iv[1][pos[0]][pos[1]][pos[2]].v * (
+				        f4_volt[0][pos[0]  ][pos[1]][pos[2]  ].v -
+				        f4_volt[0][pos[0]  ][pos[1]][pos[2]+1].v -
+				        f4_volt[2][pos[0]  ][pos[1]][pos[2]  ].v +
+				        f4_volt[2][pos[0]+1][pos[1]][pos[2]  ].v
+				    );
 
 				// z-pol
-				f4_curr[2][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii[2][pos[0]][pos[1]][pos[2]].v;
-				f4_curr[2][pos[0]][pos[1]][pos[2]].v += Op->f4_iv[2][pos[0]][pos[1]][pos[2]].v * ( f4_volt[1][pos[0]][pos[1]][pos[2]].v - f4_volt[1][pos[0]+1][pos[1]][pos[2]].v - f4_volt[0][pos[0]][pos[1]][pos[2]].v + f4_volt[0][pos[0]][pos[1]+1][pos[2]].v);
+				f4_curr[2][pos[0]][pos[1]][pos[2]].v *=
+				    Op->f4_ii[2][pos[0]][pos[1]][pos[2]].v;
+				f4_curr[2][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_iv[2][pos[0]][pos[1]][pos[2]].v * (
+				        f4_volt[1][pos[0]  ][pos[1]  ][pos[2]].v -
+				        f4_volt[1][pos[0]+1][pos[1]  ][pos[2]].v -
+				        f4_volt[0][pos[0]  ][pos[1]  ][pos[2]].v +
+				        f4_volt[0][pos[0]  ][pos[1]+1][pos[2]].v
+				    );
 			}
 
 			// for pos[2] = numVectors-1
@@ -156,20 +222,41 @@ void Engine_sse::UpdateCurrents(unsigned int startX, unsigned int numX)
 			temp.f[1] = f4_volt[1][pos[0]][pos[1]][0].f[2];
 			temp.f[2] = f4_volt[1][pos[0]][pos[1]][0].f[3];
 			temp.f[3] = 0;
-			f4_curr[0][pos[0]][pos[1]][numVectors-1].v *= Op->f4_ii[0][pos[0]][pos[1]][numVectors-1].v;
-			f4_curr[0][pos[0]][pos[1]][numVectors-1].v += Op->f4_iv[0][pos[0]][pos[1]][numVectors-1].v * ( f4_volt[2][pos[0]][pos[1]][numVectors-1].v - f4_volt[2][pos[0]][pos[1]+1][numVectors-1].v - f4_volt[1][pos[0]][pos[1]][numVectors-1].v + temp.v);
+			f4_curr[0][pos[0]][pos[1]][numVectors-1].v *=
+			    Op->f4_ii[0][pos[0]][pos[1]][numVectors-1].v;
+			f4_curr[0][pos[0]][pos[1]][numVectors-1].v +=
+			    Op->f4_iv[0][pos[0]][pos[1]][numVectors-1].v * (
+			        f4_volt[2][pos[0]][pos[1]  ][numVectors-1].v -
+			        f4_volt[2][pos[0]][pos[1]+1][numVectors-1].v -
+			        f4_volt[1][pos[0]][pos[1]  ][numVectors-1].v +
+			        temp.v
+			    );
 
 			// y-pol
 			temp.f[0] = f4_volt[0][pos[0]][pos[1]][0].f[1];
 			temp.f[1] = f4_volt[0][pos[0]][pos[1]][0].f[2];
 			temp.f[2] = f4_volt[0][pos[0]][pos[1]][0].f[3];
 			temp.f[3] = 0;
-			f4_curr[1][pos[0]][pos[1]][numVectors-1].v *= Op->f4_ii[1][pos[0]][pos[1]][numVectors-1].v;
-			f4_curr[1][pos[0]][pos[1]][numVectors-1].v += Op->f4_iv[1][pos[0]][pos[1]][numVectors-1].v * ( f4_volt[0][pos[0]][pos[1]][numVectors-1].v - temp.v - f4_volt[2][pos[0]][pos[1]][numVectors-1].v + f4_volt[2][pos[0]+1][pos[1]][numVectors-1].v);
+			f4_curr[1][pos[0]][pos[1]][numVectors-1].v *=
+			    Op->f4_ii[1][pos[0]][pos[1]][numVectors-1].v;
+			f4_curr[1][pos[0]][pos[1]][numVectors-1].v +=
+			    Op->f4_iv[1][pos[0]][pos[1]][numVectors-1].v * (
+			        f4_volt[0][pos[0]  ][pos[1]][numVectors-1].v -
+			        temp.v -
+			        f4_volt[2][pos[0]  ][pos[1]][numVectors-1].v +
+			        f4_volt[2][pos[0]+1][pos[1]][numVectors-1].v
+			    );
 
 			// z-pol
-			f4_curr[2][pos[0]][pos[1]][numVectors-1].v *= Op->f4_ii[2][pos[0]][pos[1]][numVectors-1].v;
-			f4_curr[2][pos[0]][pos[1]][numVectors-1].v += Op->f4_iv[2][pos[0]][pos[1]][numVectors-1].v * ( f4_volt[1][pos[0]][pos[1]][numVectors-1].v - f4_volt[1][pos[0]+1][pos[1]][numVectors-1].v - f4_volt[0][pos[0]][pos[1]][numVectors-1].v + f4_volt[0][pos[0]][pos[1]+1][numVectors-1].v);
+			f4_curr[2][pos[0]][pos[1]][numVectors-1].v *=
+			    Op->f4_ii[2][pos[0]][pos[1]][numVectors-1].v;
+			f4_curr[2][pos[0]][pos[1]][numVectors-1].v +=
+			    Op->f4_iv[2][pos[0]][pos[1]][numVectors-1].v * (
+			        f4_volt[1][pos[0]  ][pos[1]  ][numVectors-1].v -
+			        f4_volt[1][pos[0]+1][pos[1]  ][numVectors-1].v -
+			        f4_volt[0][pos[0]  ][pos[1]  ][numVectors-1].v +
+			        f4_volt[0][pos[0]  ][pos[1]+1][numVectors-1].v
+			    );
 		}
 		++pos[0];
 	}

FDTD/engine_sse_compressed.cpp

@@ -55,47 +55,93 @@ void Engine_SSE_Compressed::UpdateVoltages(unsigned int startX, unsigned int num
 			{
 				index = Op->m_Op_index[pos[0]][pos[1]][pos[2]];
 				// x-polarization
-				f4_volt[0][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv_Compressed[0][index].v;
-				f4_volt[0][pos[0]][pos[1]][pos[2]].v += Op->f4_vi_Compressed[0][index].v * ( f4_curr[2][pos[0]][pos[1]][pos[2]].v - f4_curr[2][pos[0]][pos[1]-shift[1]][pos[2]].v - f4_curr[1][pos[0]][pos[1]][pos[2]].v + f4_curr[1][pos[0]][pos[1]][pos[2]-1].v );
+				f4_volt[0][pos[0]][pos[1]][pos[2]].v *=
+				    Op->f4_vv_Compressed[0][index].v;
+				f4_volt[0][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_vi_Compressed[0][index].v * (
+				        f4_curr[2][pos[0]][pos[1]         ][pos[2]  ].v -
+				        f4_curr[2][pos[0]][pos[1]-shift[1]][pos[2]  ].v -
+				        f4_curr[1][pos[0]][pos[1]         ][pos[2]  ].v +
+				        f4_curr[1][pos[0]][pos[1]         ][pos[2]-1].v
+				    );
 
 				// y-polarization
-				f4_volt[1][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv_Compressed[1][index].v;
-				f4_volt[1][pos[0]][pos[1]][pos[2]].v += Op->f4_vi_Compressed[1][index].v * ( f4_curr[0][pos[0]][pos[1]][pos[2]].v - f4_curr[0][pos[0]][pos[1]][pos[2]-1].v - f4_curr[2][pos[0]][pos[1]][pos[2]].v + f4_curr[2][pos[0]-shift[0]][pos[1]][pos[2]].v);
+				f4_volt[1][pos[0]][pos[1]][pos[2]].v *=
+				    Op->f4_vv_Compressed[1][index].v;
+				f4_volt[1][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_vi_Compressed[1][index].v * (
+				        f4_curr[0][pos[0]         ][pos[1]][pos[2]  ].v -
+				        f4_curr[0][pos[0]         ][pos[1]][pos[2]-1].v -
+				        f4_curr[2][pos[0]         ][pos[1]][pos[2]  ].v +
+				        f4_curr[2][pos[0]-shift[0]][pos[1]][pos[2]  ].v
+				    );
 
 				// z-polarization
-				f4_volt[2][pos[0]][pos[1]][pos[2]].v *= Op->f4_vv_Compressed[2][index].v;
-				f4_volt[2][pos[0]][pos[1]][pos[2]].v += Op->f4_vi_Compressed[2][index].v * ( f4_curr[1][pos[0]][pos[1]][pos[2]].v - f4_curr[1][pos[0]-shift[0]][pos[1]][pos[2]].v - f4_curr[0][pos[0]][pos[1]][pos[2]].v + f4_curr[0][pos[0]][pos[1]-shift[1]][pos[2]].v);
+				f4_volt[2][pos[0]][pos[1]][pos[2]].v *=
+				    Op->f4_vv_Compressed[2][index].v;
+				f4_volt[2][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_vi_Compressed[2][index].v * (
+				        f4_curr[1][pos[0]         ][pos[1]]         [pos[2]].v -
+				        f4_curr[1][pos[0]-shift[0]][pos[1]]         [pos[2]].v -
+				        f4_curr[0][pos[0]         ][pos[1]]         [pos[2]].v +
+				        f4_curr[0][pos[0]         ][pos[1]-shift[1]][pos[2]].v
+				    );
 			}
 
 			// for pos[2] = 0
 			// x-polarization
 			index = Op->m_Op_index[pos[0]][pos[1]][0];
 #ifdef __SSE2__
-			temp.v = (__m128)_mm_slli_si128( (__m128i)f4_curr[1][pos[0]][pos[1]][numVectors-1].v, 4 );
+			temp.v = (__m128)_mm_slli_si128(
+			             (__m128i)f4_curr[1][pos[0]][pos[1]][numVectors-1].v, 4
+			         );
 #else
 			temp.f[0] = 0;
 			temp.f[1] = f4_curr[1][pos[0]][pos[1]][numVectors-1].f[0];
 			temp.f[2] = f4_curr[1][pos[0]][pos[1]][numVectors-1].f[1];
 			temp.f[3] = f4_curr[1][pos[0]][pos[1]][numVectors-1].f[2];
 #endif
-			f4_volt[0][pos[0]][pos[1]][0].v *= Op->f4_vv_Compressed[0][index].v;
-			f4_volt[0][pos[0]][pos[1]][0].v += Op->f4_vi_Compressed[0][index].v * ( f4_curr[2][pos[0]][pos[1]][0].v - f4_curr[2][pos[0]][pos[1]-shift[1]][0].v - f4_curr[1][pos[0]][pos[1]][0].v + temp.v );
+			f4_volt[0][pos[0]][pos[1]][0].v *=
+			    Op->f4_vv_Compressed[0][index].v;
+			f4_volt[0][pos[0]][pos[1]][0].v +=
+			    Op->f4_vi_Compressed[0][index].v * (
+			        f4_curr[2][pos[0]][pos[1]         ][0].v -
+			        f4_curr[2][pos[0]][pos[1]-shift[1]][0].v -
+			        f4_curr[1][pos[0]][pos[1]         ][0].v +
+			        temp.v
+			    );
 
 			// y-polarization
 #ifdef __SSE2__
-			temp.v = (__m128)_mm_slli_si128( (__m128i)f4_curr[0][pos[0]][pos[1]][numVectors-1].v, 4 );
+			temp.v = (__m128)_mm_slli_si128(
+			             (__m128i)f4_curr[0][pos[0]][pos[1]][numVectors-1].v, 4
+			         );
 #else
 			temp.f[0] = 0;
 			temp.f[1] = f4_curr[0][pos[0]][pos[1]][numVectors-1].f[0];
 			temp.f[2] = f4_curr[0][pos[0]][pos[1]][numVectors-1].f[1];
 			temp.f[3] = f4_curr[0][pos[0]][pos[1]][numVectors-1].f[2];
 #endif
-			f4_volt[1][pos[0]][pos[1]][0].v *= Op->f4_vv_Compressed[1][index].v;
-			f4_volt[1][pos[0]][pos[1]][0].v += Op->f4_vi_Compressed[1][index].v * ( f4_curr[0][pos[0]][pos[1]][0].v - temp.v - f4_curr[2][pos[0]][pos[1]][0].v + f4_curr[2][pos[0]-shift[0]][pos[1]][0].v);
+			f4_volt[1][pos[0]][pos[1]][0].v *=
+			    Op->f4_vv_Compressed[1][index].v;
+			f4_volt[1][pos[0]][pos[1]][0].v +=
+			    Op->f4_vi_Compressed[1][index].v * (
+			        f4_curr[0][pos[0]         ][pos[1]][0].v -
+			        temp.v -
+			        f4_curr[2][pos[0]         ][pos[1]][0].v +
+			        f4_curr[2][pos[0]-shift[0]][pos[1]][0].v
+			    );
 
 			// z-polarization
-			f4_volt[2][pos[0]][pos[1]][0].v *= Op->f4_vv_Compressed[2][index].v;
-			f4_volt[2][pos[0]][pos[1]][0].v += Op->f4_vi_Compressed[2][index].v * ( f4_curr[1][pos[0]][pos[1]][0].v - f4_curr[1][pos[0]-shift[0]][pos[1]][0].v - f4_curr[0][pos[0]][pos[1]][0].v + f4_curr[0][pos[0]][pos[1]-shift[1]][0].v);
+			f4_volt[2][pos[0]][pos[1]][0].v *=
+			    Op->f4_vv_Compressed[2][index].v;
+			f4_volt[2][pos[0]][pos[1]][0].v +=
+			    Op->f4_vi_Compressed[2][index].v * (
+			        f4_curr[1][pos[0]         ][pos[1]         ][0].v -
+			        f4_curr[1][pos[0]-shift[0]][pos[1]         ][0].v -
+			        f4_curr[0][pos[0]         ][pos[1]         ][0].v +
+			        f4_curr[0][pos[0]         ][pos[1]-shift[1]][0].v
+			    );
 		}
 		++pos[0];
 	}
@@ -116,47 +162,93 @@ void Engine_SSE_Compressed::UpdateCurrents(unsigned int startX, unsigned int num
 			{
 				index = Op->m_Op_index[pos[0]][pos[1]][pos[2]];
 				// x-pol
-				f4_curr[0][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii_Compressed[0][index].v;
-				f4_curr[0][pos[0]][pos[1]][pos[2]].v += Op->f4_iv_Compressed[0][index].v * ( f4_volt[2][pos[0]][pos[1]][pos[2]].v - f4_volt[2][pos[0]][pos[1]+1][pos[2]].v - f4_volt[1][pos[0]][pos[1]][pos[2]].v + f4_volt[1][pos[0]][pos[1]][pos[2]+1].v);
+				f4_curr[0][pos[0]][pos[1]][pos[2]].v *=
+				    Op->f4_ii_Compressed[0][index].v;
+				f4_curr[0][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_iv_Compressed[0][index].v * (
+				        f4_volt[2][pos[0]][pos[1]  ][pos[2]  ].v -
+				        f4_volt[2][pos[0]][pos[1]+1][pos[2]  ].v -
+				        f4_volt[1][pos[0]][pos[1]  ][pos[2]  ].v +
+				        f4_volt[1][pos[0]][pos[1]  ][pos[2]+1].v
+				    );
 
 				// y-pol
-				f4_curr[1][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii_Compressed[1][index].v;
-				f4_curr[1][pos[0]][pos[1]][pos[2]].v += Op->f4_iv_Compressed[1][index].v * ( f4_volt[0][pos[0]][pos[1]][pos[2]].v - f4_volt[0][pos[0]][pos[1]][pos[2]+1].v - f4_volt[2][pos[0]][pos[1]][pos[2]].v + f4_volt[2][pos[0]+1][pos[1]][pos[2]].v);
+				f4_curr[1][pos[0]][pos[1]][pos[2]].v *=
+				    Op->f4_ii_Compressed[1][index].v;
+				f4_curr[1][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_iv_Compressed[1][index].v * (
+				        f4_volt[0][pos[0]  ][pos[1]][pos[2]  ].v -
+				        f4_volt[0][pos[0]  ][pos[1]][pos[2]+1].v -
+				        f4_volt[2][pos[0]  ][pos[1]][pos[2]  ].v +
+				        f4_volt[2][pos[0]+1][pos[1]][pos[2]  ].v
+				    );
 
 				// z-pol
-				f4_curr[2][pos[0]][pos[1]][pos[2]].v *= Op->f4_ii_Compressed[2][index].v;
-				f4_curr[2][pos[0]][pos[1]][pos[2]].v += Op->f4_iv_Compressed[2][index].v * ( f4_volt[1][pos[0]][pos[1]][pos[2]].v - f4_volt[1][pos[0]+1][pos[1]][pos[2]].v - f4_volt[0][pos[0]][pos[1]][pos[2]].v + f4_volt[0][pos[0]][pos[1]+1][pos[2]].v);
+				f4_curr[2][pos[0]][pos[1]][pos[2]].v *=
+				    Op->f4_ii_Compressed[2][index].v;
+				f4_curr[2][pos[0]][pos[1]][pos[2]].v +=
+				    Op->f4_iv_Compressed[2][index].v * (
+				        f4_volt[1][pos[0]  ][pos[1]  ][pos[2]].v -
+				        f4_volt[1][pos[0]+1][pos[1]  ][pos[2]].v -
+				        f4_volt[0][pos[0]  ][pos[1]  ][pos[2]].v +
+				        f4_volt[0][pos[0]  ][pos[1]+1][pos[2]].v
+				    );
 			}
 
 			index = Op->m_Op_index[pos[0]][pos[1]][numVectors-1];
 			// for pos[2] = numVectors-1
 			// x-pol
 #ifdef __SSE2__
-			temp.v = (__m128)_mm_srli_si128( (__m128i)f4_volt[1][pos[0]][pos[1]][0].v, 4 );
+			temp.v = (__m128)_mm_srli_si128(
+			             (__m128i)f4_volt[1][pos[0]][pos[1]][0].v, 4
+			         );
 #else
 			temp.f[0] = f4_volt[1][pos[0]][pos[1]][0].f[1];
 			temp.f[1] = f4_volt[1][pos[0]][pos[1]][0].f[2];
 			temp.f[2] = f4_volt[1][pos[0]][pos[1]][0].f[3];
 			temp.f[3] = 0;
 #endif
-			f4_curr[0][pos[0]][pos[1]][numVectors-1].v *= Op->f4_ii_Compressed[0][index].v;
-			f4_curr[0][pos[0]][pos[1]][numVectors-1].v += Op->f4_iv_Compressed[0][index].v * ( f4_volt[2][pos[0]][pos[1]][numVectors-1].v - f4_volt[2][pos[0]][pos[1]+1][numVectors-1].v - f4_volt[1][pos[0]][pos[1]][numVectors-1].v + temp.v);
+			f4_curr[0][pos[0]][pos[1]][numVectors-1].v *=
+			    Op->f4_ii_Compressed[0][index].v;
+			f4_curr[0][pos[0]][pos[1]][numVectors-1].v +=
+			    Op->f4_iv_Compressed[0][index].v * (
+			        f4_volt[2][pos[0]][pos[1]  ][numVectors-1].v -
+			        f4_volt[2][pos[0]][pos[1]+1][numVectors-1].v -
+			        f4_volt[1][pos[0]][pos[1]  ][numVectors-1].v +
+			        temp.v
+			    );
 
 			// y-pol
 #ifdef __SSE2__
-			temp.v = (__m128)_mm_srli_si128( (__m128i)f4_volt[0][pos[0]][pos[1]][0].v, 4 );
+			temp.v = (__m128)_mm_srli_si128(
+			             (__m128i)f4_volt[0][pos[0]][pos[1]][0].v, 4
+			         );
 #else
 			temp.f[0] = f4_volt[0][pos[0]][pos[1]][0].f[1];
 			temp.f[1] = f4_volt[0][pos[0]][pos[1]][0].f[2];
 			temp.f[2] = f4_volt[0][pos[0]][pos[1]][0].f[3];
 			temp.f[3] = 0;
 #endif
-			f4_curr[1][pos[0]][pos[1]][numVectors-1].v *= Op->f4_ii_Compressed[1][index].v;
-			f4_curr[1][pos[0]][pos[1]][numVectors-1].v += Op->f4_iv_Compressed[1][index].v * ( f4_volt[0][pos[0]][pos[1]][numVectors-1].v - temp.v - f4_volt[2][pos[0]][pos[1]][numVectors-1].v + f4_volt[2][pos[0]+1][pos[1]][numVectors-1].v);
+			f4_curr[1][pos[0]][pos[1]][numVectors-1].v *=
+			    Op->f4_ii_Compressed[1][index].v;
+			f4_curr[1][pos[0]][pos[1]][numVectors-1].v +=
+			    Op->f4_iv_Compressed[1][index].v * (
+			        f4_volt[0][pos[0]  ][pos[1]][numVectors-1].v -
+			        temp.v -
+			        f4_volt[2][pos[0]  ][pos[1]][numVectors-1].v +
+			        f4_volt[2][pos[0]+1][pos[1]][numVectors-1].v
+			    );
 
 			// z-pol
-			f4_curr[2][pos[0]][pos[1]][numVectors-1].v *= Op->f4_ii_Compressed[2][index].v;
-			f4_curr[2][pos[0]][pos[1]][numVectors-1].v += Op->f4_iv_Compressed[2][index].v * ( f4_volt[1][pos[0]][pos[1]][numVectors-1].v - f4_volt[1][pos[0]+1][pos[1]][numVectors-1].v - f4_volt[0][pos[0]][pos[1]][numVectors-1].v + f4_volt[0][pos[0]][pos[1]+1][numVectors-1].v);
+			f4_curr[2][pos[0]][pos[1]][numVectors-1].v *=
+			    Op->f4_ii_Compressed[2][index].v;
+			f4_curr[2][pos[0]][pos[1]][numVectors-1].v +=
+			    Op->f4_iv_Compressed[2][index].v * (
+			        f4_volt[1][pos[0]  ][pos[1]  ][numVectors-1].v -
+			        f4_volt[1][pos[0]+1][pos[1]  ][numVectors-1].v -
+			        f4_volt[0][pos[0]  ][pos[1]  ][numVectors-1].v +
+			        f4_volt[0][pos[0]  ][pos[1]+1][numVectors-1].v
+			    );
 		}
 		++pos[0];
 	}

FDTD/excitation.cpp

@@ -53,7 +53,7 @@ void Excitation::Reset( double timestep )
 	m_foi = 0;
 }
 
-bool Excitation::SetupGaussianPulse(double f0, double fc)
+void Excitation::SetupGaussianPulse(double f0, double fc)
 {
 	m_Excit_Type = Excitation::GaissianPulse;
 	m_f0 = f0;
@@ -62,7 +62,7 @@ bool Excitation::SetupGaussianPulse(double f0, double fc)
 	m_SignalPeriod = 0;
 }
 
-bool Excitation::SetupSinusoidal(double f0)
+void Excitation::SetupSinusoidal(double f0)
 {
 	m_Excit_Type = Excitation::Sinusoidal;
 	m_f0 = f0;
@@ -70,21 +70,21 @@ bool Excitation::SetupSinusoidal(double f0)
 	m_SignalPeriod = 1/f0;
 }
 
-bool Excitation::SetupDiracPulse(double fmax)
+void Excitation::SetupDiracPulse(double fmax)
 {
 	m_Excit_Type = Excitation::DiracPulse;
 	m_SignalPeriod = 0;
 	m_f_max = fmax;
 }
 
-bool Excitation::SetupStepExcite(double fmax)
+void Excitation::SetupStepExcite(double fmax)
 {
 	m_Excit_Type = Excitation::Step;
 	m_SignalPeriod = 0;
 	m_f_max = fmax;
 }
 
-bool Excitation::SetupCustomExcite(string str, double f0, double fmax)
+void Excitation::SetupCustomExcite(string str, double f0, double fmax)
 {
 	m_Excit_Type = Excitation::CustomExcite;
 	m_CustomExc_Str = str;
@@ -137,7 +137,7 @@ unsigned int Excitation::GetMaxExcitationTimestep() const
 {
 	FDTD_FLOAT maxAmp=0;
 	unsigned int maxStep=0;
-	for (unsigned int n=1; n<Length+1; ++n)
+	for (unsigned int n=0; n<Length; ++n)
 	{
 		if (fabs(Signal_volt[n])>maxAmp)
 		{
@@ -152,7 +152,7 @@ void Excitation::CalcGaussianPulsExcitation(double f0, double fc, int nTS)
 {
 	if (dT==0) return;
 
-	Length = (unsigned int)(2.0 * 9.0/(2.0*PI*fc) / dT);
+	Length = (unsigned int)ceil(2.0 * 9.0/(2.0*PI*fc) / dT);
 	if (Length>(unsigned int)nTS)
 	{
 		cerr << "Operator::CalcGaussianPulsExcitation: Requested excitation pusle would be " << Length << " timesteps or " << Length * dT << " s long. Cutting to max number of timesteps!" << endl;
@@ -160,13 +160,11 @@ void Excitation::CalcGaussianPulsExcitation(double f0, double fc, int nTS)
 	}
 	delete[] Signal_volt;
 	delete[] Signal_curr;
-	Signal_volt = new FDTD_FLOAT[Length+1];
-	Signal_curr = new FDTD_FLOAT[Length+1];
-	Signal_volt[0]=0.0;
-	Signal_curr[0]=0.0;
-	for (unsigned int n=1; n<Length+1; ++n)
+	Signal_volt = new FDTD_FLOAT[Length];
+	Signal_curr = new FDTD_FLOAT[Length];
+	for (unsigned int n=0; n<Length; ++n)
 	{
-		double t = (n-1)*dT;
+		double t = n*dT;
 		Signal_volt[n] = cos(2.0*PI*f0*(t-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*t/3.0-3,2));
 		t += 0.5*dT;
 		Signal_curr[n] = cos(2.0*PI*f0*(t-9.0/(2.0*PI*fc)))*exp(-1*pow(2.0*PI*fc*t/3.0-3,2));
@@ -182,12 +180,12 @@ void Excitation::CalcDiracPulsExcitation()
 {
 	if (dT==0) return;
 
-	Length = 1;
+	Length = 2;
 //	cerr << "Operator::CalcDiracPulsExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
 	delete[] Signal_volt;
 	delete[] Signal_curr;
-	Signal_volt = new FDTD_FLOAT[Length+1];
-	Signal_curr = new FDTD_FLOAT[Length+1];
+	Signal_volt = new FDTD_FLOAT[Length];
+	Signal_curr = new FDTD_FLOAT[Length];
 	Signal_volt[0]=0.0;
 	Signal_volt[1]=1.0;
 	Signal_curr[0]=0.0;
@@ -203,11 +201,11 @@ void Excitation::CalcStepExcitation()
 {
 	if (dT==0) return;
 
-	Length = 1;
+	Length = 2;
 	delete[] Signal_volt;
 	delete[] Signal_curr;
-	Signal_volt = new FDTD_FLOAT[Length+1];
-	Signal_curr = new FDTD_FLOAT[Length+1];
+	Signal_volt = new FDTD_FLOAT[Length];
+	Signal_curr = new FDTD_FLOAT[Length];
 	Signal_volt[0]=1.0;
 	Signal_volt[1]=1.0;
 	Signal_curr[0]=1.0;
@@ -228,10 +226,9 @@ void Excitation::CalcCustomExcitation(double f0, int nTS, string signal)
 //	cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << ExciteLength << " timesteps" << endl;
 	delete[] Signal_volt;
 	delete[] Signal_curr;
-	Signal_volt = new FDTD_FLOAT[Length+1];
-	Signal_curr = new FDTD_FLOAT[Length+1];
-	Signal_volt[0]=0.0;
-	Signal_curr[0]=0.0;
+	Signal_volt = new FDTD_FLOAT[Length];
+	Signal_curr = new FDTD_FLOAT[Length];
+	setlocale(LC_NUMERIC, "en_US.UTF-8");
 	FunctionParser fParse;
 	fParse.AddConstant("pi", 3.14159265358979323846);
 	fParse.AddConstant("e", 2.71828182845904523536);
@@ -242,9 +239,9 @@ void Excitation::CalcCustomExcitation(double f0, int nTS, string signal)
 		exit(1);
 	}
 	double vars[1];
-	for (unsigned int n=1; n<Length+1; ++n)
+	for (unsigned int n=0; n<Length; ++n)
 	{
-		vars[0] = (n-1)*dT;
+		vars[0] = n*dT;
 		Signal_volt[n] = fParse.Eval(vars);
 		vars[0] += 0.5*dT;
 		Signal_curr[n] = fParse.Eval(vars);
@@ -260,17 +257,16 @@ void Excitation::CalcSinusExcitation(double f0, int nTS)
 	if (dT==0) return;
 	if (nTS<=0) return;
 
-	Length = (unsigned int)(2.0/f0/dT);
-	//cerr << "Operator::CalcSinusExcitation: Length of the excite signal: " << Length << " timesteps " << Length*dT << "s" << endl;
+	Length = (unsigned int)round(2.0/f0/dT);
 	delete[] Signal_volt;
 	delete[] Signal_curr;
-	Signal_volt = new FDTD_FLOAT[Length+1];
-	Signal_curr = new FDTD_FLOAT[Length+1];
+	Signal_volt = new FDTD_FLOAT[Length];
+	Signal_curr = new FDTD_FLOAT[Length];
 	Signal_volt[0]=0.0;
 	Signal_curr[0]=0.0;
-	for (unsigned int n=1; n<Length+1; ++n)
+	for (unsigned int n=1; n<Length; ++n)
 	{
-		double t = (n-1)*dT;
+		double t = n*dT;
 		Signal_volt[n] = sin(2.0*PI*f0*t);
 		t += 0.5*dT;
 		Signal_curr[n] = sin(2.0*PI*f0*t);
@@ -286,8 +282,8 @@ void Excitation::DumpVoltageExcite(string filename)
 	file.open( filename.c_str() );
 	if (file.fail())
 		return;
-	for (unsigned int n=1; n<Length+1; ++n)
-		file << (n-1)*dT << "\t" << Signal_volt[n] << "\n";
+	for (unsigned int n=0; n<Length; ++n)
+		file << n*dT << "\t" << Signal_volt[n] << "\n";
 	file.close();
 }
 
@@ -297,8 +293,8 @@ void Excitation::DumpCurrentExcite(string filename)
 	file.open( filename.c_str() );
 	if (file.fail())
 		return;
-	for (unsigned int n=1; n<Length+1; ++n)
-		file << (n-1)*dT + 0.5*dT << "\t" << Signal_curr[n] << "\n";
+	for (unsigned int n=0; n<Length; ++n)
+		file << n*dT + 0.5*dT << "\t" << Signal_curr[n] << "\n";
 	file.close();
 }
 

FDTD/excitation.h

@@ -31,11 +31,11 @@ public:
 
 	virtual void Reset( double timestep );
 
-	bool SetupGaussianPulse(double f0, double fc);
-	bool SetupSinusoidal(double f0);
-	bool SetupDiracPulse(double fmax);
-	bool SetupStepExcite(double fmax);
-	bool SetupCustomExcite(std::string str, double f0, double fmax);
+	void SetupGaussianPulse(double f0, double fc);
+	void SetupSinusoidal(double f0);
+	void SetupDiracPulse(double fmax);
+	void SetupStepExcite(double fmax);
+	void SetupCustomExcite(std::string str, double f0, double fmax);
 
 	double GetCenterFreq() {return m_f0;}
 	double GetCutOffFreq() {return m_fc;}
@@ -64,7 +64,7 @@ public:
 	//! Get the length of the excitation signal
 	unsigned int GetLength() const {return Length;}
 
-	//! Get the max frequeny excited by this signal
+	//! Get the max frequency excited by this signal
 	double GetMaxFrequency() const {return m_f_max;}
 
 	//! Get the frequency of interest

FDTD/extensions/CMakeLists.txt

@@ -24,6 +24,8 @@ set(SOURCES
   ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_tfsf.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_steadystate.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_steadystate.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/operator_ext_lumpedRLC.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/engine_ext_lumpedRLC.cpp
   PARENT_SCOPE
 )
 

FDTD/extensions/engine_ext_cylindermultigrid.cpp

@@ -55,7 +55,7 @@ void Engine_Ext_CylinderMultiGrid::DoPreVoltageUpdates()
 	if (!m_IsBase)
 	{
 		//cerr << "child: volt wait on base " << endl;
-		m_WaitOnBase->wait(); //wait on base to finisch current sync and/or to finisch voltage updates, than start child voltage updates
+		m_WaitOnBase->wait(); //wait on base to finish current sync and/or to finish voltage updates, than start child voltage updates
 	}
 }
 
@@ -69,14 +69,14 @@ void Engine_Ext_CylinderMultiGrid::Apply2Voltages()
 	if (m_IsBase)
 	{
 		m_WaitOnBase->wait();	//base voltage updates are done, tell child to start its voltage updates
-		m_WaitOnChild->wait();	//wait for child to finisch its updates
-		SyncVoltages();			//child is finisch, run sync and go to current updates next
+		m_WaitOnChild->wait();	//wait for child to finish its updates
+		SyncVoltages();			//child is finish, run sync and go to current updates next
 		m_WaitOnSync->wait();	//sync is done... move on and tell child to move on...
 	}
 	else
 	{
 		m_WaitOnChild->wait();	//child is finished voltage updates, will tell base to run sync
-		m_WaitOnSync->wait();	//wait for base to finisch sync before going to wait for current updates
+		m_WaitOnSync->wait();	//wait for base to finish sync before going to wait for current updates
 	}
 }
 
@@ -125,7 +125,7 @@ void Engine_Ext_CylinderMultiGrid::DoPreCurrentUpdates()
 	if (!m_IsBase)
 	{
 		//cerr << "child: curr wait on base " << endl;
-		m_WaitOnBase->wait(); //wait on base to finisch voltage sync and current updates, than start child current updates
+		m_WaitOnBase->wait(); //wait on base to finish voltage sync and current updates, than start child current updates
 	}
 }
 
@@ -140,14 +140,14 @@ void Engine_Ext_CylinderMultiGrid::Apply2Current()
 	{
 		//cerr << "Base: curr wait on base done, wait on sync" << endl;
 		m_WaitOnBase->wait();	//base current updates are done, tell child to start its current updates
-		m_WaitOnChild->wait();	//wait for child to finisch its updates
-		SyncCurrents();			//child is finisch, run sync and go to voltage updates next
+		m_WaitOnChild->wait();	//wait for child to finish its updates
+		SyncCurrents();			//child is finish, run sync and go to voltage updates next
 		m_WaitOnSync->wait();	//sync is done... move on and tell child to move on...
 	}
 	else
 	{
 		m_WaitOnChild->wait();	//child is finished current updates, will tell base to run sync...
-		m_WaitOnSync->wait();	//wait for base to finisch sync before going to wait for next voltage updates
+		m_WaitOnSync->wait();	//wait for base to finish sync before going to wait for next voltage updates
 		//cerr << "Child: curr done, wait on sync" << endl;
 	}
 }

FDTD/extensions/engine_ext_excitation.cpp

@@ -54,7 +54,7 @@ void Engine_Ext_Excitation::Apply2Voltages()
 				exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n];
 				exc_pos *= (exc_pos>0);
 				exc_pos %= p;
-				exc_pos *= (exc_pos<=(int)length);
+				exc_pos *= (exc_pos<(int)length);
 				ny = m_Op_Exc->Volt_dir[n];
 				pos[0]=m_Op_Exc->Volt_index[0][n];
 				pos[1]=m_Op_Exc->Volt_index[1][n];
@@ -71,7 +71,7 @@ void Engine_Ext_Excitation::Apply2Voltages()
 				exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n];
 				exc_pos *= (exc_pos>0);
 				exc_pos %= p;
-				exc_pos *= (exc_pos<=(int)length);
+				exc_pos *= (exc_pos<(int)length);
 				ny = m_Op_Exc->Volt_dir[n];
 				pos[0]=m_Op_Exc->Volt_index[0][n];
 				pos[1]=m_Op_Exc->Volt_index[1][n];
@@ -87,7 +87,7 @@ void Engine_Ext_Excitation::Apply2Voltages()
 				exc_pos = numTS - (int)m_Op_Exc->Volt_delay[n];
 				exc_pos *= (exc_pos>0);
 				exc_pos %= p;
-				exc_pos *= (exc_pos<=(int)length);
+				exc_pos *= (exc_pos<(int)length);
 				ny = m_Op_Exc->Volt_dir[n];
 				pos[0]=m_Op_Exc->Volt_index[0][n];
 				pos[1]=m_Op_Exc->Volt_index[1][n];
@@ -124,7 +124,7 @@ void Engine_Ext_Excitation::Apply2Current()
 				exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n];
 				exc_pos *= (exc_pos>0);
 				exc_pos %= p;
-				exc_pos *= (exc_pos<=(int)length);
+				exc_pos *= (exc_pos<(int)length);
 				ny = m_Op_Exc->Curr_dir[n];
 				pos[0]=m_Op_Exc->Curr_index[0][n];
 				pos[1]=m_Op_Exc->Curr_index[1][n];
@@ -141,7 +141,7 @@ void Engine_Ext_Excitation::Apply2Current()
 				exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n];
 				exc_pos *= (exc_pos>0);
 				exc_pos %= p;
-				exc_pos *= (exc_pos<=(int)length);
+				exc_pos *= (exc_pos<(int)length);
 				ny = m_Op_Exc->Curr_dir[n];
 				pos[0]=m_Op_Exc->Curr_index[0][n];
 				pos[1]=m_Op_Exc->Curr_index[1][n];
@@ -157,7 +157,7 @@ void Engine_Ext_Excitation::Apply2Current()
 				exc_pos = numTS - (int)m_Op_Exc->Curr_delay[n];
 				exc_pos *= (exc_pos>0);
 				exc_pos %= p;
-				exc_pos *= (exc_pos<=(int)length);
+				exc_pos *= (exc_pos<(int)length);
 				ny = m_Op_Exc->Curr_dir[n];
 				pos[0]=m_Op_Exc->Curr_index[0][n];
 				pos[1]=m_Op_Exc->Curr_index[1][n];

FDTD/extensions/engine_ext_lumpedRLC.cpp

@@ -0,0 +1,195 @@
+/*
+*	Additional
+*	Copyright (C) 2023 Gadi Lahav (gadi@rfwithcare.com)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "engine_ext_lumpedRLC.h"
+#include "operator_ext_lumpedRLC.h"
+
+#include "FDTD/engine_sse.h"
+
+Engine_Ext_LumpedRLC::Engine_Ext_LumpedRLC(Operator_Ext_LumpedRLC* op_ext_RLC) : Engine_Extension(op_ext_RLC)
+{
+	// Local pointer of the operator.
+	m_Op_Ext_RLC = op_ext_RLC;
+
+	v_Vdn		= new FDTD_FLOAT*[3];
+	v_Jn		= new FDTD_FLOAT*[3];
+
+	// No additional allocations are required if there are no actual lumped elements.
+	if (!(m_Op_Ext_RLC->RLC_count))
+		return;
+
+	// Initialize ADE containers for currents and voltages
+	v_Il 		= new FDTD_FLOAT[m_Op_Ext_RLC->RLC_count];
+
+	for (unsigned int posIdx = 0 ; posIdx < m_Op_Ext_RLC->RLC_count ; ++posIdx)
+		v_Il[posIdx] 	= 0.0;
+
+	for (unsigned int k = 0 ; k < 3 ; k++)
+	{
+		v_Vdn[k] = new FDTD_FLOAT[m_Op_Ext_RLC->RLC_count];
+		v_Jn[k] = new FDTD_FLOAT[m_Op_Ext_RLC->RLC_count];
+
+		for (unsigned int posIdx = 0 ; posIdx < m_Op_Ext_RLC->RLC_count ; ++posIdx)
+		{
+			v_Jn[k][posIdx] = 0.0;
+			v_Vdn[k][posIdx] = 0.0;;
+		}
+	}
+
+}
+
+Engine_Ext_LumpedRLC::~Engine_Ext_LumpedRLC()
+{
+	// Only delete if values were allocated in the first place
+	if (m_Op_Ext_RLC->RLC_count)
+	{
+		delete[] v_Il;
+
+		for (unsigned int k = 0 ; k < 3 ; k++)
+		{
+			delete[] v_Vdn[k];
+			delete[] v_Jn[k];
+		}
+	}
+
+	delete[] v_Vdn;
+	delete[] v_Jn;
+
+	v_Il	= NULL;
+
+	v_Vdn	= NULL;
+	v_Jn	= NULL;
+
+	m_Op_Ext_RLC = NULL;
+
+
+}
+
+void Engine_Ext_LumpedRLC::DoPreVoltageUpdates()
+{
+	unsigned int **pos = m_Op_Ext_RLC->v_RLC_pos;
+	int *dir = m_Op_Ext_RLC->v_RLC_dir;
+
+	// Iterate Vd containers
+	FDTD_FLOAT	*v_temp;
+	v_temp = v_Vdn[2];
+	v_Vdn[2] = v_Vdn[1];
+	v_Vdn[1] = v_Vdn[0];
+	v_Vdn[0] = v_temp;
+
+	// In pre-process, only update the parallel inductor current:
+	for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+		v_Il[pIdx] += (m_Op_Ext_RLC->v_RLC_i2v[pIdx])*(m_Op_Ext_RLC->v_RLC_ilv[pIdx])*v_Vdn[1][pIdx];
+
+	return;
+}
+
+void Engine_Ext_LumpedRLC::Apply2Voltages()
+{
+	unsigned int **pos = m_Op_Ext_RLC->v_RLC_pos;
+	int *dir = m_Op_Ext_RLC->v_RLC_dir;
+
+	// Iterate J containers
+	FDTD_FLOAT	*v_temp;
+	v_temp = v_Jn[2];
+	v_Jn[2] = v_Jn[1];
+	v_Jn[1] = v_Jn[0];
+	v_Jn[0] = v_temp;
+
+
+	// Read engine calculated node voltage
+	switch (m_Eng->GetType())
+	{
+		case Engine::BASIC:
+		{
+			for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+				v_Vdn[0][pIdx] = m_Eng->Engine::GetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx]);
+
+			break;
+		}
+		case Engine::SSE:
+		{
+			Engine_sse* eng_sse = (Engine_sse*)m_Eng;
+			for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+				v_Vdn[0][pIdx] = eng_sse->Engine_sse::GetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx]);
+
+			break;
+		}
+		default:
+		{
+			for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+				v_Vdn[0][pIdx] = m_Eng->GetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx]);;
+
+			break;
+		}
+	}
+
+	// Post process: Calculate node voltage with respect to the lumped RLC auxilliary quantity, J
+	for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+	{
+		// Calculate updated node voltage, with series and parallel additions
+		v_Vdn[0][pIdx] = 	(m_Op_Ext_RLC->v_RLC_vvd[pIdx])*(
+							v_Vdn[0][pIdx] - v_Il[pIdx]						// Addition for Parallel inductor
+							+
+							(m_Op_Ext_RLC->v_RLC_vv2[pIdx])*v_Vdn[2][pIdx]	// Vd[n-2] addition
+							+
+							(m_Op_Ext_RLC->v_RLC_vj1[pIdx])*v_Jn[1][pIdx]	// J[n-1] addition
+							+
+							(m_Op_Ext_RLC->v_RLC_vj2[pIdx])*v_Jn[2][pIdx]);	// J[n-2] addition
+
+		// Update J[0]
+		v_Jn[0][pIdx] =	(m_Op_Ext_RLC->v_RLC_ib0[pIdx])*(v_Vdn[0][pIdx] - v_Vdn[2][pIdx])
+						-
+						((m_Op_Ext_RLC->v_RLC_b1[pIdx])*(m_Op_Ext_RLC->v_RLC_ib0[pIdx]))*v_Jn[1][pIdx]
+						-
+						((m_Op_Ext_RLC->v_RLC_b2[pIdx])*(m_Op_Ext_RLC->v_RLC_ib0[pIdx]))*v_Jn[2][pIdx];
+	}
+
+
+	// Update node voltage
+	switch (m_Eng->GetType())
+	{
+		case Engine::BASIC:
+		{
+			for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+				m_Eng->Engine::SetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx],v_Vdn[0][pIdx]);
+
+			break;
+		}
+		case Engine::SSE:
+		{
+			Engine_sse* eng_sse = (Engine_sse*)m_Eng;
+			for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+				eng_sse->Engine_sse::SetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx],v_Vdn[0][pIdx]);
+
+			break;
+		}
+		default:
+		{
+			for (unsigned int pIdx = 0 ; pIdx < m_Op_Ext_RLC->RLC_count ; pIdx++)
+				m_Eng->SetVolt(dir[pIdx],pos[0][pIdx],pos[1][pIdx],pos[2][pIdx],v_Vdn[0][pIdx]);
+
+			break;
+		}
+	}
+
+
+	return;
+}
+
+

FDTD/extensions/engine_ext_lumpedRLC.h

@@ -0,0 +1,54 @@
+/*
+*	Copyright (C) 2023 Gadi Lahav (gadi@rfwithcare.com)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef ENGINE_EXT_LUMPEDRLC_H
+#define ENGINE_EXT_LUMPEDRLC_H
+
+#include "engine_extension.h"
+#include "FDTD/engine.h"
+#include "FDTD/operator.h"
+
+class Operator_Ext_LumpedRLC;
+
+class Engine_Ext_LumpedRLC : public Engine_Extension
+{
+	friend class Operator_Ext_LumpedRLC;
+	friend class Operator;
+	friend class ContinuousStructure;
+
+public:
+
+	Engine_Ext_LumpedRLC(Operator_Ext_LumpedRLC *op_ext_RLC);
+	virtual ~Engine_Ext_LumpedRLC();
+
+	virtual void DoPreVoltageUpdates();
+	virtual void Apply2Voltages();
+
+protected:
+	Operator_Ext_LumpedRLC* m_Op_Ext_RLC;
+
+	// Auxilliary containers
+
+	// Array setup: volt_C_ADE[mesh_pos]
+	FDTD_FLOAT *v_Il;		// Container for current on inductor- Parallel RLC
+
+	FDTD_FLOAT **v_Vdn;		// Container for nodal vd at [n],[n-1],[n-2]
+	FDTD_FLOAT **v_Jn;		// Container for nodal J at [n],[n-1],[n-2]
+
+};
+
+#endif // ENGINE_EXT_LORENTZMATERIAL_H

FDTD/extensions/engine_ext_tfsf.cpp

@@ -44,7 +44,7 @@ void Engine_Ext_TFSF::DoPostVoltageUpdates()
 	{
 		if ( numTS < n )
 			m_DelayLookup[n]=0;
-		else if ((numTS-n > length) && (p==0))
+		else if ((numTS-n >= length) && (p==0))
 			m_DelayLookup[n]=0;
 		else
 			m_DelayLookup[n] = numTS - n;
@@ -134,7 +134,7 @@ void Engine_Ext_TFSF::DoPostCurrentUpdates()
 	{
 		if ( numTS < n )
 			m_DelayLookup[n]=0;
-		else if ((numTS-n > length) && (p==0))
+		else if ((numTS-n >= length) && (p==0))
 			m_DelayLookup[n]=0;
 		else
 			m_DelayLookup[n] = numTS - n;

FDTD/extensions/engine_extension.h

@@ -41,27 +41,27 @@ public:
 
 	virtual void SetNumberOfThreads(int nrThread);
 
-	//! This methode will be called __before__ the main engine does the usual voltage updates. This methode may __not__ change the engine voltages!!!
+	//! This method will be called __before__ the main engine does the usual voltage updates. This method may __not__ change the engine voltages!!!
 	virtual void DoPreVoltageUpdates() {}
 	virtual void DoPreVoltageUpdates(int threadID);
-	//! This methode will be called __after__ the main engine does the usual voltage updates. This methode may __not__ change the engine voltages!!!
+	//! This method will be called __after__ the main engine does the usual voltage updates. This method may __not__ change the engine voltages!!!
 	virtual void DoPostVoltageUpdates() {}
 	virtual void DoPostVoltageUpdates(int threadID);
-	//! This methode will be called __after__ all updates to the voltages and extensions and may add/set its results to the engine voltages, but may __not__ rely on the current value of the engine voltages!!!
+	//! This method will be called __after__ all updates to the voltages and extensions and may add/set its results to the engine voltages, but may __not__ rely on the current value of the engine voltages!!!
 	virtual void Apply2Voltages() {}
 	virtual void Apply2Voltages(int threadID);
 
-	//! This methode will be called __before__ the main engine does the usual current updates. This methode may __not__ change the engine current!!!
+	//! This method will be called __before__ the main engine does the usual current updates. This method may __not__ change the engine current!!!
 	virtual void DoPreCurrentUpdates() {}
 	virtual void DoPreCurrentUpdates(int threadID);
-	//! This methode will be called __after__ the main engine does the usual current updates. This methode may __not__ change the engine current!!!
+	//! This method will be called __after__ the main engine does the usual current updates. This method may __not__ change the engine current!!!
 	virtual void DoPostCurrentUpdates() {}
 	virtual void DoPostCurrentUpdates(int threadID);
-	//! This methode will be called __after__ all updates to the current and extensions and may add/set its results to the engine current, but may __not__ rely on the current value of the engine current!!!
+	//! This method will be called __after__ all updates to the current and extensions and may add/set its results to the engine current, but may __not__ rely on the current value of the engine current!!!
 	virtual void Apply2Current() {}
 	virtual void Apply2Current(int threadID);
 
-	//! Set the Engine to this extention. This will usually done automatically by Engine::AddExtension
+	//! Set the Engine to this extension. This will usually done automatically by Engine::AddExtension
 	virtual void SetEngine(Engine* eng) {m_Eng=eng;}
 
 	//! Get the priority for this extension

FDTD/extensions/operator_ext_excitation.cpp

@@ -135,14 +135,13 @@ bool Operator_Ext_Excitation::BuildExtension()
 
 	CSPropExcitation* elec=NULL;
 	CSProperties* prop=NULL;
-	int priority=0;
 
 	unsigned int numLines[] = {m_Op->GetNumberOfLines(0,true),m_Op->GetNumberOfLines(1,true),m_Op->GetNumberOfLines(2,true)};
-
 	for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
 	{
 		for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
 		{
+			vector<CSPrimitives*> vPrims = m_Op->GetPrimitivesBoundBox(-1, pos[1], pos[2], CSProperties::EXCITATION);
 			for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
 			{
 				//electric field excite
@@ -156,14 +155,12 @@ bool Operator_Ext_Excitation::BuildExtension()
 					if (m_CC_R0_included && (n==1) && (pos[0]==0))
 						continue;
 
-					for (size_t p=0; p<vec_prop.size(); ++p)
+					CSProperties* prop = CSX->GetPropertyByCoordPriority(volt_coord, vPrims, true);
+					if (prop)
 					{
-						prop = vec_prop.at(p);
 						elec = prop->ToExcitation();
 						if (elec==NULL)
 							continue;
-						if (prop->CheckCoordInPrimitive(volt_coord,priority,true))
-						{
 						if ((elec->GetActiveDir(n)) && ( (elec->GetExcitType()==0) || (elec->GetExcitType()==1) ))//&& (pos[n]<numLines[n]-1))
 						{
 							amp = elec->GetWeightedExcitation(n,volt_coord)*m_Op->GetEdgeLength(n,pos);// delta[n]*gridDelta;
@@ -184,7 +181,6 @@ bool Operator_Ext_Excitation::BuildExtension()
 						}
 					}
 				}
-				}
 
 				//magnetic field excite
 				for (int n=0; n<3; ++n)
@@ -193,14 +189,12 @@ bool Operator_Ext_Excitation::BuildExtension()
 						continue;  //skip the last H-Line which is outside the FDTD-domain
 					if (m_Op->GetYeeCoords(n,pos,curr_coord,true)==false)
 						continue;
-					for (size_t p=0; p<vec_prop.size(); ++p)
+					CSProperties* prop = CSX->GetPropertyByCoordPriority(curr_coord, vPrims, true);
+					if (prop)
 					{
-						prop = vec_prop.at(p);
 						elec = prop->ToExcitation();
 						if (elec==NULL)
 							continue;
-						if (prop->CheckCoordInPrimitive(curr_coord,priority,true))
-						{
 						if ((elec->GetActiveDir(n)) && ( (elec->GetExcitType()==2) || (elec->GetExcitType()==3) ))
 						{
 							amp = elec->GetWeightedExcitation(n,curr_coord)*m_Op->GetEdgeLength(n,pos,true);// delta[n]*gridDelta;
@@ -221,7 +215,6 @@ bool Operator_Ext_Excitation::BuildExtension()
 						}
 					}
 				}
-				}
 
 			}
 		}

FDTD/extensions/operator_ext_lorentzmaterial.cpp

@@ -225,6 +225,9 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
 
 //						CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), true);
 						CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, vPrims, true);
+
+						if (prop==NULL) continue;
+
 						if ((mat = prop->ToLorentzMaterial()))
 						{
 							w_plasma = mat->GetEpsPlasmaFreqWeighted(order,n,coord) * 2 * PI;
@@ -277,6 +280,9 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
 
 //						CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::METAL | CSProperties::MATERIAL), true);
 						CSProperties* prop = m_Op->GetGeometryCSX()->GetPropertyByCoordPriority(coord, vPrims, true);
+
+						if (prop==NULL) continue;
+
 						if ((mat = prop->ToLorentzMaterial()))
 						{
 							w_plasma = mat->GetMuePlasmaFreqWeighted(order,n,coord) * 2 * PI;
@@ -308,7 +314,7 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
 							if (L_D[n]>0)
 							{
 								v_int[n].push_back((2.0*L_D[n]-dT*R_D[n])/(2.0*L_D[n]+dT*R_D[n]));
-								// check for r==0 in clyindrical coords and get special VI cooefficient
+								// check for r==0 in clyindrical coords and get special VI coefficient
 								if (m_CC_R0_included && n==2 && pos[0]==0)
 									v_ext[n].push_back(dT/(L_D[n]+dT*R_D[n]/2.0)*m_Op_Cyl->m_Cyl_Ext->vi_R0[pos[2]]);
 								else

FDTD/extensions/operator_ext_lumpedRLC.cpp

@@ -0,0 +1,514 @@
+/*
+*	Copyright (C) 2023 Gadi Lahav (gadi@rfwithcare.com)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "../operator.h"
+#include "tools/array_ops.h"
+#include "tools/constants.h"
+//#include "cond_sheet_parameter.h"
+#include "tools/AdrOp.h"
+
+#include "operator_ext_lumpedRLC.h"
+#include "engine_ext_lumpedRLC.h"
+
+#include "CSPrimBox.h"
+#include "CSProperties.h"
+#include "CSPropLumpedElement.h"
+
+#define COPY_V2A(V,A) std::copy(V.begin(),V.end(),A)
+
+Operator_Ext_LumpedRLC::Operator_Ext_LumpedRLC(Operator* op) : Operator_Extension(op)
+{
+	// Parallel circuit coefficients
+	v_RLC_ilv = NULL;
+	v_RLC_i2v = NULL;
+
+	// Series circuit coefficients
+	v_RLC_vv2	= NULL;	// Coefficient for [n-2] time of Vd update in Vd equation
+	v_RLC_vj1	= NULL;	// Coefficient for [n-1] time of J update in Vd equation
+	v_RLC_vj2	= NULL;	// Coefficient for [n-2] time of J update in Vd equation
+	v_RLC_vvd	= NULL;	// Coefficient to multiply all Vd in the Vd update equation
+	v_RLC_ib0	= NULL;	// Inverse of beta_0
+	v_RLC_b1	= NULL;	// beta_1
+	v_RLC_b2	= NULL;	// beta_2
+
+	// Additional containers
+	v_RLC_dir = NULL;
+	v_RLC_pos = NULL;
+
+	RLC_count = 0;
+}
+
+Operator_Ext_LumpedRLC::Operator_Ext_LumpedRLC(Operator* op, Operator_Ext_LumpedRLC* op_ext) : Operator_Extension(op,op_ext)
+{
+	// Parallel circuit coefficients
+	v_RLC_ilv = NULL;
+	v_RLC_i2v = NULL;
+
+	// Series circuit coefficients
+	v_RLC_vv2	= NULL;	// Coefficient for [n-2] time of Vd update in Vd equation
+	v_RLC_vj1	= NULL;	// Coefficient for [n-1] time of J update in Vd equation
+	v_RLC_vj2	= NULL;	// Coefficient for [n-2] time of J update in Vd equation
+	v_RLC_vvd	= NULL;	// Coefficient to multiply all Vd in the Vd update equation
+	v_RLC_ib0	= NULL;	// Inverse of beta_0
+	v_RLC_b1	= NULL;	// beta_1
+	v_RLC_b2	= NULL;	// beta_2
+
+	// Additional containers
+	v_RLC_dir = NULL;
+	v_RLC_pos = NULL;
+
+	RLC_count = 0;
+}
+
+Operator_Ext_LumpedRLC::~Operator_Ext_LumpedRLC()
+{
+	if (this->RLC_count)
+	{
+		// Parallel circuit coefficients
+		delete[] v_RLC_ilv;
+		delete[] v_RLC_i2v;
+
+		// Series circuit coefficients
+		delete[] v_RLC_vv2;
+		delete[] v_RLC_vj1;
+		delete[] v_RLC_vj2;
+		delete[] v_RLC_vvd;
+		delete[] v_RLC_ib0;
+		delete[] v_RLC_b1;
+		delete[] v_RLC_b2;
+
+		// Additional containers
+		delete[] v_RLC_dir;
+
+		for (unsigned int dIdx = 0 ; dIdx < 3 ; dIdx++)
+			delete[] v_RLC_pos[dIdx];
+
+		delete[] v_RLC_pos;
+
+	}
+}
+
+Operator_Extension* Operator_Ext_LumpedRLC::Clone(Operator* op)
+{
+	if (dynamic_cast<Operator_Ext_LumpedRLC*>(this)==NULL)
+		return NULL;
+	return new Operator_Ext_LumpedRLC(op, this);
+}
+
+bool Operator_Ext_LumpedRLC::BuildExtension()
+{
+	double 	dT = m_Op->GetTimestep();
+
+	double	fMax = 	m_Op->GetExcitationSignal()->GetCenterFreq()
+					+
+					m_Op->GetExcitationSignal()->GetCutOffFreq();
+
+	unsigned int 	pos[] = {0,0,0};
+
+	vector<CSProperties*> cs_props;
+
+	int 			dir;
+	CSPropLumpedElement::LEtype lumpedType;
+
+	vector<unsigned int> 	v_pos[3];
+
+	vector<int>		v_dir;
+
+	vector<double>  v_ilv;
+	vector<double>	v_i2v;
+
+	vector<double>	v_vv2;
+	vector<double>	v_vj1;
+	vector<double>	v_vj2;
+	vector<double>	v_vvd;
+	vector<double>	v_ib0;
+	vector<double>	v_b1;
+	vector<double>	v_b2;
+
+	// Lumped RLC parameters
+	double R, L, C;
+
+	// clear all vectors to initialize them
+	for (unsigned int dIdx = 0 ; dIdx < 3 ; dIdx++)
+		v_pos[dIdx].clear();
+
+	v_dir.clear();
+
+	v_ilv.clear();
+	v_i2v.clear();
+
+	v_vv2.clear();
+	v_vj1.clear();
+	v_vj2.clear();
+	v_vvd.clear();
+	v_ib0.clear();
+	v_b1.clear();
+	v_b2.clear();
+
+	// Obtain from CSX (continuous structure) all the lumped RLC properties
+	// Properties are material properties, not the objects themselves
+	cs_props = m_Op->CSX->GetPropertyByType(CSProperties::LUMPED_ELEMENT);
+
+	// Iterate through various properties. In theory, there should be a property set per-
+	// primitive, as each lumped element should have it's own unique properties.
+	for(size_t n = 0 ; n < cs_props.size() ; ++n)
+	{
+		// Cast current property to lumped RLC property continuous structure properties
+		CSPropLumpedElement* cs_RLC_props = dynamic_cast<CSPropLumpedElement*>(cs_props.at(n));
+		if (cs_RLC_props==NULL)
+			return false; //sanity check: this should never happen!
+
+		// Store direction and type
+		dir = cs_RLC_props->GetDirection();
+		lumpedType = cs_RLC_props->GetLEtype();
+
+		//		if (lumpedType == LEtype::INVALID
+		if (lumpedType == CSPropLumpedElement::INVALID)
+		{
+			cerr << "Operator_Ext_LumpedRLC::BuildExtension(): Warning: RLCtype is invalid! considering as parallel. "
+					<< " ID: " << cs_RLC_props->GetID() << " @ Property: " << cs_RLC_props->GetName() << endl;
+			lumpedType = CSPropLumpedElement::PARALLEL;
+		}
+
+		// Extract R, L and C from property class
+		C = cs_RLC_props->GetCapacity();
+		if (C < 0)
+			C = NAN;
+		R = cs_RLC_props->GetResistance();
+		if (R < 0)
+			R = NAN;
+		L = cs_RLC_props->GetInductance();
+		if (L < 0)
+			L = NAN;
+
+		// Check that this is a lumped RLC
+		if (!(this->IsLElumpedRLC(cs_RLC_props)))
+			continue;
+
+		if ((dir < 0) || (dir > 2))
+		{
+			cerr << "Operator_Ext_LumpedRLC::Calc_LumpedElements(): Warning: Lumped Element direction is invalid! skipping. "
+					<< " ID: " << cs_RLC_props->GetID() << " @ Property: " << cs_RLC_props->GetName() << endl;
+			continue;
+		}
+
+		// Initialize other two direction containers
+		int dir_p1 = (dir + 1) % 3;
+		int dir_p2 = (dir + 2) % 3;
+
+		// Now iterate through primitive(s). I still think there should be only one per-
+		// material definition, but maybe I'm wrong...
+		vector<CSPrimitives*> cs_RLC_prims = cs_RLC_props->GetAllPrimitives();
+
+		for (size_t boxIdx = 0 ; boxIdx < cs_RLC_prims.size() ; ++boxIdx)
+		{
+			CSPrimBox* cBox = dynamic_cast<CSPrimBox*>(cs_RLC_prims.at(boxIdx));
+
+			if (cBox)
+			{
+
+				// Get box start and stop positions
+				unsigned int 	uiStart[3],
+								uiStop[3];
+
+
+				// snap to the native coordinate system
+				int Snap_Dimension =
+						m_Op->SnapBox2Mesh(
+								cBox->GetStartCoord()->GetCoords(m_Op->m_MeshType),	// Start Coord
+								cBox->GetStopCoord()->GetCoords(m_Op->m_MeshType),	// Stop Coord
+								uiStart,	// Start Index
+								uiStop,		// Stop Index
+								false,		// Dual (doublet) Grid?
+								true);		// Full mesh?
+
+				// Verify that snapped dimension is correct
+				if (Snap_Dimension<=0)
+				{
+					if (Snap_Dimension>=-1)
+						cerr << "Operator_Ext_LumpedRLC::BuildExtension(): Warning: Lumped RLC snapping failed! Dimension is: " << Snap_Dimension << " skipping. "
+								<< " ID: " << cs_RLC_prims.at(boxIdx)->GetID() << " @ Property: " << cs_RLC_props->GetName() << endl;
+					// Snap_Dimension == -2 means outside the simulation domain --> no special warning, but box probably marked as unused!
+					continue;
+				}
+
+				// Verify that in the direction of the current propagation, the size isn't zero.
+				if (uiStart[dir]==uiStop[dir])
+				{
+					cerr << "Operator_Ext_LumpedRLC::BuildExtension(): Warning: Lumped RLC with zero (snapped) length is invalid! skipping. "
+							<< " ID: " << cs_RLC_prims.at(boxIdx)->GetID() << " @ Property: " << cs_RLC_props->GetName() << endl;
+					continue;
+				}
+
+				// Calculate number of cells per-direction
+				unsigned int 	Ncells_0 = uiStop[dir] - uiStart[dir],
+								Ncells_1 = uiStop[dir_p1] - uiStart[dir_p1] + 1,
+								Ncells_2 = uiStop[dir_p2] - uiStart[dir_p2] + 1;
+
+				// All cells in directions 1 and 2 are considered parallel connection
+				unsigned int Npar = Ncells_1*Ncells_2;
+
+				// Separate elements such that individual elements can be calculated.
+				double	dL = L*Npar/Ncells_0,
+						dR = R*Npar/Ncells_0,
+						dG = (1/R)*Ncells_0/Npar,
+						dC = C*Ncells_0/Npar;
+
+				double	ib0	= 2.0*dT*dC/(4.0*dL*dC + 2.0*dT*dR*dC + dT*dT),
+						b1	= (dT*dT - 4.0*dL*dC)/(dT*dC),
+						b2	= (4.0*dL*dC - 2.0*dT*dR*dC + dT*dT)/(2.0*dT*dC);
+
+				// Special case: If this is a parallel resonant circuit, and there is no
+				// parallel resistor, use zero conductivity. May be risky when low-loss
+				// simulations are involved
+				if (lumpedType == CSPropLumpedElement::PARALLEL)
+					if (R == 0.0)
+						dG = 0.0;
+
+				int iPos = 0;
+
+				double Zmin,Zcd_min;
+
+				// In the various positions, update the capacitors and "inverse" resistors
+				for (pos[dir] = uiStart[dir] ; pos[dir] < uiStop[dir] ; ++pos[dir])
+				{
+					for (pos[dir_p1] = uiStart[dir_p1] ; pos[dir_p1] <= uiStop[dir_p1] ; ++pos[dir_p1])
+					{
+						for (pos[dir_p2] = uiStart[dir_p2] ; pos[dir_p2] <= uiStop[dir_p2] ; ++pos[dir_p2])
+						{
+							iPos = m_Op->MainOp->SetPos(pos[0],pos[1],pos[2]);
+
+
+							// Separate to two different cases. Parallel and series
+							switch (lumpedType)
+							{
+								case CSPropLumpedElement::PARALLEL:
+									// Update capacitor either way.
+									if (dC > 0)
+										m_Op->EC_C[dir][iPos] = dC;
+									else
+										// This case takes the "natural" capacitor into account.
+										dC = m_Op->EC_C[dir][iPos];
+
+									v_i2v.push_back((dT/dC)/(1.0 + dT*dG/(2.0*dC)));
+
+									// Update conductivity
+									if (R >= 0)
+										m_Op->EC_G[dir][iPos] = dG;
+
+									// Update coefficients with respect to the parallel inductance
+									if (L > 0)
+										v_ilv.push_back(dT/dL);
+									else
+										v_ilv.push_back(0.0);
+
+									// Take into account the case that the "natural" capacitor is too small
+									// with respect to the inductor or the resistor, and add a warning.
+									if (dC == 0)
+									{
+										double Cd = m_Op->EC_C[dir][iPos];
+										Zmin = max(dR,2*PI*fMax*dL);
+										Zcd_min = 1.0/(2.0*PI*fMax*Cd);
+
+										// Check if the "parasitic" capcitance is not small enough
+										if (Zcd_min < LUMPED_RLC_Z_FACT*Zmin)
+										{
+											Cd = 1.0/(2*PI*fMax*Zmin*LUMPED_RLC_Z_FACT);
+											m_Op->EC_C[dir][iPos] = Cd;
+										}
+									}
+
+									v_vv2.push_back(0.0);
+									v_vj1.push_back(0.0);
+									v_vj2.push_back(0.0);
+									v_vvd.push_back(1.0);
+									v_ib0.push_back(0.0);
+									v_b1.push_back(0.0);
+									v_b2.push_back(0.0);
+
+									// Update with discrete component values of
+									m_Op->Calc_ECOperatorPos(dir,pos);
+
+									v_dir.push_back(dir);
+
+									break;
+
+								case CSPropLumpedElement::SERIES:
+									m_Op->EC_G[dir][iPos] = 0.0;
+
+									// is a series inductor, modeled separately.
+									FDTD_FLOAT Cd = m_Op->EC_C[dir][iPos];
+
+									// Calculate minimum impedance, at maximum frequency
+									Zmin = sqrt(pow(dR,2) + pow(2*PI*fMax*dL - 1.0/(dC*2*PI*fMax),2));
+									Zcd_min = 1.0/(2.0*PI*fMax*Cd);
+
+									// Check if the "parasitic" capcitance is not small enough
+									if (Zcd_min < LUMPED_RLC_Z_FACT*Zmin)
+									{
+										Cd = 1.0/(2*PI*fMax*Zmin*LUMPED_RLC_Z_FACT);
+										m_Op->EC_C[dir][iPos] = Cd;
+									}
+
+									// No contribution from parallel inductor
+									v_ilv.push_back(0.0);
+									v_i2v.push_back(0.0);
+
+									// Contributions from series resistor and inductor
+									v_vv2.push_back(0.5*dT*ib0/Cd);
+									v_vj1.push_back(0.5*dT*(b1*ib0 - 1.0)/Cd);
+									v_vj2.push_back(0.5*dT*b2*ib0/Cd);
+									v_vvd.push_back(1.0/(1.0 + 0.5*dT*ib0/Cd));
+									v_ib0.push_back(ib0);
+									v_b1.push_back(b1);
+									v_b2.push_back(b2);
+
+									m_Op->Calc_ECOperatorPos(dir,pos);
+
+									v_dir.push_back(dir);
+
+									break;
+							}
+
+							// Store position and direction
+							for (unsigned int dIdx = 0 ; dIdx < 3 ; ++dIdx)
+								v_pos[dIdx].push_back(pos[dIdx]);
+
+						}
+					}
+				}
+
+				// Build metallic caps
+				if (cs_RLC_props->GetCaps())
+					for (pos[dir_p1] = uiStart[dir_p1] ; pos[dir_p1] <= uiStop[dir_p1] ; ++pos[dir_p1])
+					{
+						for (pos[dir_p2] = uiStart[dir_p2] ; pos[dir_p2] <= uiStop[dir_p2] ; ++pos[dir_p2])
+						{
+							pos[dir]=uiStart[dir];
+							if (pos[dir_p1]<uiStop[dir_p1])
+							{
+								m_Op->SetVV(dir_p1,pos[0],pos[1],pos[2], 0 );
+								m_Op->SetVI(dir_p1,pos[0],pos[1],pos[2], 0 );
+								++(m_Op->m_Nr_PEC[dir_p1]);
+							}
+
+							if (pos[dir_p2]<uiStop[dir_p2])
+							{
+								m_Op->SetVV(dir_p2,pos[0],pos[1],pos[2], 0 );
+								m_Op->SetVI(dir_p2,pos[0],pos[1],pos[2], 0 );
+								++(m_Op->m_Nr_PEC[dir_p2]);
+							}
+
+							pos[dir]=uiStop[dir];
+							if (pos[dir_p1]<uiStop[dir_p1])
+							{
+								m_Op->SetVV(dir_p1,pos[0],pos[1],pos[2], 0 );
+								m_Op->SetVI(dir_p1,pos[0],pos[1],pos[2], 0 );
+								++(m_Op->m_Nr_PEC[dir_p1]);
+							}
+
+							if (pos[dir_p2]<uiStop[dir_p2])
+							{
+								m_Op->SetVV(dir_p2,pos[0],pos[1],pos[2], 0 );
+								m_Op->SetVI(dir_p2,pos[0],pos[1],pos[2], 0 );
+								++(m_Op->m_Nr_PEC[dir_p2]);
+							}
+						}
+					}
+
+
+				// Mark as used
+				cBox->SetPrimitiveUsed(true);
+			}
+		}
+
+	}
+
+	// Start data storage
+	RLC_count = v_dir.size();
+
+	// values
+	if (RLC_count)
+	{
+		// Allocate space to all variables
+		v_RLC_dir 	= new int[RLC_count];
+
+		// Parallel circuit coefficients
+		v_RLC_ilv 	= new FDTD_FLOAT[RLC_count];
+		v_RLC_i2v 	= new FDTD_FLOAT[RLC_count];
+
+		// Series circuit coefficients
+		v_RLC_vv2 = new FDTD_FLOAT[RLC_count];
+		v_RLC_vj1 = new FDTD_FLOAT[RLC_count];
+		v_RLC_vj2 = new FDTD_FLOAT[RLC_count];
+		v_RLC_vvd = new FDTD_FLOAT[RLC_count];
+		v_RLC_ib0 = new FDTD_FLOAT[RLC_count];
+		v_RLC_b1 = new FDTD_FLOAT[RLC_count];
+		v_RLC_b2 = new FDTD_FLOAT[RLC_count];
+
+		v_RLC_pos = new unsigned int*[3];
+		for (unsigned int dIdx = 0 ; dIdx < 3 ; ++dIdx)
+			v_RLC_pos[dIdx] = new unsigned int[RLC_count];
+
+		// Copy all vectors to arrays
+		COPY_V2A(v_dir, v_RLC_dir);
+
+		COPY_V2A(v_ilv, v_RLC_ilv);
+		COPY_V2A(v_i2v, v_RLC_i2v);
+
+		COPY_V2A(v_vv2,v_RLC_vv2);
+		COPY_V2A(v_vj1,v_RLC_vj1);
+		COPY_V2A(v_vj2,v_RLC_vj2);
+		COPY_V2A(v_vvd,v_RLC_vvd);
+		COPY_V2A(v_ib0,v_RLC_ib0);
+		COPY_V2A(v_b1,v_RLC_b1);
+		COPY_V2A(v_b2,v_RLC_b2);
+
+		for (unsigned int dIdx = 0 ; dIdx < 3 ; ++dIdx)
+			COPY_V2A(v_pos[dIdx],v_RLC_pos[dIdx]);
+	}
+
+	return true;
+}
+
+Engine_Extension* Operator_Ext_LumpedRLC::CreateEngineExtention()
+{
+	Engine_Ext_LumpedRLC* eng_ext_RLC = new Engine_Ext_LumpedRLC(this);
+	return eng_ext_RLC;
+}
+
+void Operator_Ext_LumpedRLC::ShowStat(ostream &ostr)  const
+{
+	Operator_Extension::ShowStat(ostr);
+	string On_Off[2] = {"Off", "On"};
+
+	ostr << "Active cells\t\t: " << RLC_count << endl;
+}
+
+bool Operator_Ext_LumpedRLC::IsLElumpedRLC(const CSPropLumpedElement* const p_prop)
+{
+	CSPropLumpedElement::LEtype lumpedType = p_prop->GetLEtype();
+
+	double L = p_prop->GetInductance();
+
+	bool isParallelRLC = (lumpedType == CSPropLumpedElement::PARALLEL) && (L > 0.0);
+	bool isSeriesRLC = lumpedType == CSPropLumpedElement::SERIES;
+
+	// This needs to be something that isn't a parallel RC circuit to add data to this extension.
+	return isParallelRLC || isSeriesRLC;
+}
+

FDTD/extensions/operator_ext_lumpedRLC.h

@@ -0,0 +1,88 @@
+/*
+*	Copyright (C) 2023 Gadi Lahav (gadi@rfwithcare.com)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPERATOR_EXT_LUMPEDRLC_H_
+#define OPERATOR_EXT_LUMPEDRLC_H_
+
+#include "vector"
+
+#include "FDTD/operator.h"
+#include "operator_extension.h"
+#include "operator_ext_cylinder.h"
+
+#include "engine_ext_lumpedRLC.h"
+
+#define LUMPED_RLC_Z_FACT 20.0
+
+class Operator_Ext_LumpedRLC : public Operator_Extension
+{
+	friend class Engine_Ext_LumpedRLC;
+public:
+	Operator_Ext_LumpedRLC(Operator* op);
+
+	virtual ~Operator_Ext_LumpedRLC();
+
+	virtual Operator_Extension* Clone(Operator* op);
+
+	virtual bool BuildExtension();
+
+	virtual Engine_Extension* CreateEngineExtention();
+
+	virtual bool IsCylinderCoordsSave(bool closedAlpha, bool R0_included) const {UNUSED(closedAlpha); UNUSED(R0_included); return true;}
+	virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child); return true;}
+	virtual bool IsMPISave() const {return true;}
+
+	virtual string GetExtensionName() const {return string("Series\\Parallel Lumped RLC load");}
+
+	virtual void ShowStat(ostream &ostr) const;
+
+	virtual bool IsLElumpedRLC(const CSPropLumpedElement* const p_prop);
+
+protected:
+	//! Copy constructor
+	Operator_Ext_LumpedRLC(Operator* op, Operator_Ext_LumpedRLC* op_ext);
+
+	// ADE update coefficients, array setup: coeff[mesh_pos_index]
+
+	// Parallel circuit coefficients
+	FDTD_FLOAT		*v_RLC_ilv;
+	FDTD_FLOAT		*v_RLC_i2v;
+
+	// Series circuit coefficients
+	FDTD_FLOAT		*v_RLC_vv2;		// Coefficient for [n-2] time of Vd update in Vd equation
+	FDTD_FLOAT		*v_RLC_vj1;		// Coefficient for [n-1] time of J update in Vd equation
+	FDTD_FLOAT		*v_RLC_vj2;		// Coefficient for [n-2] time of J update in Vd equation
+	FDTD_FLOAT		*v_RLC_vvd;		// Coefficient to multiply all Vd in the Vd update equation
+	FDTD_FLOAT		*v_RLC_ib0;		// Inverse of beta_0
+	FDTD_FLOAT		*v_RLC_b1;		// beta_1
+	FDTD_FLOAT		*v_RLC_b2;		// beta_2
+
+	// Additional containers
+	int				*v_RLC_dir;
+	unsigned int			**v_RLC_pos;
+
+	// Vector length indicator
+	unsigned int			RLC_count;
+
+
+
+
+};
+
+
+
+#endif /* OPERATOR_EXT_LUMPEDRLC_H_ */

FDTD/extensions/operator_ext_tfsf.cpp

@@ -159,7 +159,7 @@ bool Operator_Ext_TFSF::BuildExtension()
 		else
 			m_PhVel=m_Op->CalcNumericPhaseVelocity(m_Start,m_Stop,m_PropDir,m_Frequency);
 
-		if ((m_PhVel<0) || (m_PhVel>__C0__/ref_index) || isnan(m_PhVel))
+		if ((m_PhVel<0) || (m_PhVel>__C0__/ref_index) || std::isnan(m_PhVel))
 		{
 			cerr << "Operator_Ext_TFSF::BuildExtension: Warning, invalid phase velocity found, resetting to c0! " << endl;
 			m_PhVel = __C0__/ref_index;

FDTD/extensions/operator_ext_upml.cpp

@@ -25,6 +25,7 @@ using namespace std;
 
 Operator_Ext_UPML::Operator_Ext_UPML(Operator* op) : Operator_Extension(op)
 {
+	setlocale(LC_NUMERIC, "en_US.UTF-8");
 	m_GradingFunction = new FunctionParser();
 	//default grading function
 	SetGradingFunction(" -log(1e-6)*log(2.5)/(2*dl*Z*(pow(2.5,W/dl)-1)) * pow(2.5, D/dl) ");
@@ -280,7 +281,7 @@ bool Operator_Ext_UPML::SetGradingFunction(string func)
 	int res = m_GradingFunction->Parse(m_GradFunc.c_str(), "D,dl,W,Z,N");
 	if (res < 0) return true;
 
-	cerr << "Operator_Ext_UPML::SetGradingFunction: Warning, an error occured parsing the pml grading function (see below) ..." << endl;
+	cerr << "Operator_Ext_UPML::SetGradingFunction: Warning, an error occurred parsing the pml grading function (see below) ..." << endl;
 	cerr << func << "\n" << string(res, ' ') << "^\n" << m_GradingFunction->ErrorMsg() << "\n";
 	return false;
 }
@@ -403,7 +404,9 @@ bool Operator_Ext_UPML::BuildExtension()
 					CalcGradingKappa(n, pos,__Z0__ ,kappa_v ,kappa_i);
 					nP = (n+1)%3;
 					nPP = (n+2)%3;
-					if ((kappa_v[0]+kappa_v[1]+kappa_v[2])!=0)
+					// if eff_Mat[1] > 1e3 assume a metal and disable PML to continue a signal layer
+					// sigma > 1000 S/m is "very lossy" already...
+					if (((kappa_v[0]+kappa_v[1]+kappa_v[2])!=0) && (eff_Mat[1]<1e3))
 					{
 						//check if pos is on PEC
 						if ( (m_Op->GetVV(n,pos[0],pos[1],pos[2]) + m_Op->GetVI(n,pos[0],pos[1],pos[2])) != 0 )

FDTD/extensions/operator_ext_upml.h

@@ -23,11 +23,11 @@
 
 class FunctionParser;
 
-//! Operator extension implemention an uniaxial perfectly matched layer (upml)
+//! Operator extension implementation an uniaxial perfectly matched layer (upml)
 /*
   The priority for this extension should be the highest of all extensions since this operator will use the main engine to perform vital parts in the upml implementation.
   Therefore the voltages and currents as well as the operator are replaced during these update process.
-  This extension is propably incompatible with the most other extensions operating in the same regions.
+  This extension is probably incompatible with the most other extensions operating in the same regions.
   */
 class Operator_Ext_UPML : public Operator_Extension
 {

FDTD/extensions/operator_extension.h

@@ -40,7 +40,7 @@ public:
 	//! Create a clone of this extension, will return NULL if this is impossible
 	/*!
 		Create a clone of this extension, will return NULL if this is impossible (e.g. derived extension has no clone method and copy-constructor)...
-		BuildExtension has to be called separatly!
+		BuildExtension has to be called separately!
 	*/
 	virtual Operator_Extension* Clone(Operator* op) {UNUSED(op); return NULL;}
 

FDTD/openems_fdtd_mpi.cpp

@@ -210,7 +210,7 @@ bool openEMS_FDTD_MPI::SetupMPI()
 	if (numProcs!=m_NumProc)
 	{
 		if (m_MyID==0)
-			cerr << "openEMS_FDTD_MPI::SetupMPI: Error: Requested splits require " << numProcs << " processes, but only " << m_NumProc << " were found! Exit! " << endl;
+			cerr << "openEMS_FDTD_MPI::SetupMPI: Error: Requested splits require " << numProcs << " processes, but " << m_NumProc << " were found! Exit! " << endl;
 		exit(10);
 	}
 
@@ -256,8 +256,8 @@ bool openEMS_FDTD_MPI::SetupMPI()
 						grid->AddDiscLine(2, m_Original_Grid->GetLine(2,n) );
 
 					m_MPI_Op->SetSplitPos(0,m_SplitNumber[0].at(i));
-					m_MPI_Op->SetSplitPos(1,m_SplitNumber[1].at(i));
-					m_MPI_Op->SetSplitPos(2,m_SplitNumber[2].at(i));
+					m_MPI_Op->SetSplitPos(1,m_SplitNumber[1].at(j));
+					m_MPI_Op->SetSplitPos(2,m_SplitNumber[2].at(k));
 
 					if (i>0)
 						m_MPI_Op->SetNeighborDown(0,procTable[i-1][j][k]);
@@ -334,7 +334,7 @@ unsigned int openEMS_FDTD_MPI::GetNextStep()
 
 	int local_step=step;
 
-	//find the smallest next step requestes by all processings
+	//find the smallest next step requests by all processings
 	MPI_Reduce(&local_step, &step, 1, MPI_INT, MPI_MIN, 0, MPI_COMM_WORLD);
 	//send the smallest next step to all
 	MPI_Bcast(&step, 1, MPI_INT, 0, MPI_COMM_WORLD);
@@ -397,19 +397,20 @@ bool openEMS_FDTD_MPI::SetupProcessing()
 			if (dynamic_cast<ProcessIntegral*>(proc)!=NULL)
 			{
 				//type is integral processing --> disable! Needs to be fixed!
-				cerr << "openEMS_FDTD_MPI::SetupProcessing(): Warning: Processing: " << proc->GetName() << " occures multiple times and is being deactivated..." << endl;
+				cerr << "openEMS_FDTD_MPI::SetupProcessing(): Warning: Processing: " << proc->GetName() << " occurs multiple times and is being deactivated..." << endl;
+				cerr << "openEMS_FDTD_MPI::SetupProcessing(): Note: Processing: Make sure that there are no splits inside probes or sources." << endl;
 				deactivate = true;
 				rename = false;
 			}
 			if (dynamic_cast<ProcessFields*>(proc)!=NULL)
 			{
 				//type is field processing --> renameing! Needs to be fixed!
-				cerr << "openEMS_FDTD_MPI::SetupProcessing(): Warning: Processing: " << proc->GetName() << " occures multiple times and is being renamed..." << endl;
+				cerr << "openEMS_FDTD_MPI::SetupProcessing(): Warning: Processing: " << proc->GetName() << " occurs multiple times and is being renamed..." << endl;
 				deactivate = false;
 				rename = true;
 			}
 		}
-		//broadcast informations to all
+		//broadcast information to all
 		MPI_Bcast(&deactivate, 1, MPI::BOOL, 0, MPI_COMM_WORLD);
 		MPI_Bcast(&rename, 1, MPI::BOOL, 0, MPI_COMM_WORLD);
 		if (deactivate)

FDTD/operator.cpp

@@ -373,13 +373,12 @@ int Operator::SnapLine2Mesh(const double* start, const double* stop, unsigned in
 		return ret;
 
 	//fixme, do we need to do something about start or stop being outside the field domain?
-	//maybe caclulate the intersection point and snap to that?
+	//maybe calculate the intersection point and snap to that?
 	//it seems to work like this as well...
 
 	return ret;
 }
 
-
 Grid_Path Operator::FindPath(double start[], double stop[])
 {
 	Grid_Path path;
@@ -450,7 +449,7 @@ Grid_Path Operator::FindPath(double start[], double stop[])
 			currPos[minDir]+=-1;
 			minPos[minDir]-=1;
 		}
-		//check validity of current postion
+		//check validity of current position
 		for (int n=0;n<3;++n)
 			if (currPos[n]>=numLines[n])
 			{
@@ -790,7 +789,7 @@ void Operator::DumpMaterial2File(string filename)
 	delete vtk_Writer;
 }
 
- bool Operator::SetupCSXGrid(CSRectGrid* grid)
+bool Operator::SetupCSXGrid(CSRectGrid* grid)
  {
 	 for (int n=0; n<3; ++n)
 	 {
@@ -1395,9 +1394,9 @@ bool Operator::AverageMatCellCenter(int ny, const unsigned int* pos, double* Eff
 	if (EffMat[3]) EffMat[3]=length / EffMat[3];
 
 	for (int n=0; n<4; ++n)
-		if (isnan(EffMat[n]) || isinf(EffMat[n]))
+		if (std::isnan(EffMat[n]) || std::isinf(EffMat[n]))
 		{
-			cerr << "Operator::" << __func__ << ": Error, an effective material parameter is not a valid result, this should NOT have happend... exit..." << endl;
+			cerr << "Operator::" << __func__ << ": Error, an effective material parameter is not a valid result, this should NOT have happened... exit..." << endl;
 			cerr << ny << "@" << n << " : " << pos[0] << "," << pos[1] << ","  << pos[2] << endl;
 			exit(0);
 		}
@@ -1508,9 +1507,9 @@ bool Operator::AverageMatQuarterCell(int ny, const unsigned int* pos, double* Ef
 	if (EffMat[3]) EffMat[3]=length / EffMat[3];
 
 	for (int n=0; n<4; ++n)
-		if (isnan(EffMat[n]) || isinf(EffMat[n]))
+		if (std::isnan(EffMat[n]) || std::isinf(EffMat[n]))
 		{
-			cerr << "Operator::" << __func__ << ": Error, An effective material parameter is not a valid result, this should NOT have happend... exit..." << endl;
+			cerr << "Operator::" << __func__ << ": Error, An effective material parameter is not a valid result, this should NOT have happened... exit..." << endl;
 			cerr << ny << "@" << n << " : " << pos[0] << "," << pos[1] << ","  << pos[2] << endl;
 			exit(0);
 		}
@@ -1536,11 +1535,15 @@ bool Operator::Calc_EffMatPos(int ny, const unsigned int* pos, double* EffMat, v
 bool Operator::Calc_LumpedElements()
 {
 	vector<CSProperties*> props = CSX->GetPropertyByType(CSProperties::LUMPED_ELEMENT);
+
 	for (size_t i=0;i<props.size();++i)
 	{
+
 		CSPropLumpedElement* PLE = dynamic_cast<CSPropLumpedElement*>(props.at(i));
+
 		if (PLE==NULL)
 			return false; //sanity check: this should never happen!
+
 		vector<CSPrimitives*> prims = PLE->GetAllPrimitives();
 		for (size_t bn=0;bn<prims.size();++bn)
 		{
@@ -1555,7 +1558,12 @@ bool Operator::Calc_LumpedElements()
 				if (R<0)
 					R = NAN;
 
-				if ((isnan(R)) && (isnan(C)))
+				// If this is not a parallel RC, skip this.
+				if (!(this->IsLEparRC(PLE)))
+					continue;
+
+
+				if ((std::isnan(R)) && (std::isnan(C)))
 				{
 					cerr << "Operator::Calc_LumpedElements(): Warning: Lumped Element R or C not specified! skipping. "
 							<< " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
@@ -1696,13 +1704,25 @@ bool Operator::Calc_LumpedElements()
 
 			}
 			else
-				cerr << "Operator::Calc_LumpedElements(): Warning: Primitves other than boxes are not supported for lumped elements! skipping "
+				cerr << "Operator::Calc_LumpedElements(): Warning: Primitives other than boxes are not supported for lumped elements! skipping "
 						<< prims.at(bn)->GetTypeName() << " ID: " << prims.at(bn)->GetID() << " @ Property: " << PLE->GetName() << endl;
 		}
 	}
 	return true;
 }
 
+bool Operator::IsLEparRC(const CSPropLumpedElement* const p_prop)
+{
+	CSPropLumpedElement::LEtype lumpedType = p_prop->GetLEtype();
+
+	double L = p_prop->GetInductance();
+
+	bool IsParallelRC = (lumpedType == CSPropLumpedElement::PARALLEL) && !(L > 0.0);
+
+	// This needs to be something that isn't a parallel RC circuit to add data to this extension.
+	return IsParallelRC;
+}
+
 void Operator::Init_EC()
 {
 	for (int n=0; n<3; ++n)

FDTD/operator.h

@@ -33,12 +33,16 @@ class Operator : public Operator_Base
 {
 	friend class Engine;
 	friend class Engine_Interface_FDTD;
-	friend class Operator_Ext_LorentzMaterial; //we need to find a way around this... friend class Operator_Extension only would be nice
-	friend class Operator_Ext_ConductingSheet; //we need to find a way around this... friend class Operator_Extension only would be nice
+	friend class Operator_Ext_LorentzMaterial; 	// We need to find a way around this... friend class Operator_Extension only would be nice
+	friend class Operator_Ext_ConductingSheet; 	// We need to find a way around this... friend class Operator_Extension only would be nice
 	friend class Operator_Ext_PML_SF_Plane;
 	friend class Operator_Ext_Excitation;
 	friend class Operator_Ext_UPML;
 	friend class Operator_Ext_Cylinder;
+	friend class Operator_Ext_LumpedRLC;		// Gadi: I now know why the two previous remarks are here.
+
+	// So apparaently I have to use functionality from operator
+	// in my "lumpedRLC" class. This is ugly...
 public:
 	enum DebugFlags {None=0,debugMaterial=1,debugOperator=2,debugPEC=4};
 
@@ -134,7 +138,7 @@ public:
 	//! Get the volume of an FDTD cell
 	virtual double GetCellVolume(const unsigned int pos[3], bool dualMesh = false) const;
 
-	//! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos
+	//! Get the area around an edge for a given direction \a n and a given mesh position \a pos
 	/*!
 		This will return the area around an edge with a given direction, measured at the middle of the edge.
 		In a cartesian mesh this is equal to the NodeArea, may be different in other coordinate systems.
@@ -143,7 +147,7 @@ public:
 
 	virtual unsigned int SnapToMeshLine(int ny, double coord, bool &inside, bool dualMesh=false, bool fullMesh=false) const;
 
-	//! Snap the given coodinates to mesh indices
+	//! Snap the given coordinates to mesh indices
 	virtual bool SnapToMesh(const double* coord, unsigned int* uicoord, bool dualMesh=false, bool fullMesh=false, bool* inside=NULL) const;
 
 	//! Snap a given box to the FDTD mesh
@@ -223,7 +227,7 @@ protected:
 	/*!
 	  Get the raw disc delta for a given position and direction.
 	  The result will be positive if a disc delta inside the simulation domain is requested.
-	  The result will be the negative value of the first or last disc delta respectivly if the position is outside the field domain.
+	  The result will be the negative value of the first or last disc delta respectively if the position is outside the field domain.
 	  */
 	virtual double GetRawDiscDelta(int ny, const int pos) const;
 
@@ -244,6 +248,9 @@ protected:
 	//! Calculate and setup lumped elements
 	virtual bool Calc_LumpedElements();
 
+	//! Condition to determine if this is a lumped RC, to invoke Calc_LumpedElements
+	virtual bool IsLEparRC(const CSPropLumpedElement* const p_prop);
+
 	//! Store the size of the applied boundary conditions
 	int m_BC_Size[6];
 

FDTD/operator_cylinder.cpp

@@ -46,7 +46,7 @@ Operator_Cylinder::~Operator_Cylinder()
 Engine* Operator_Cylinder::CreateEngine()
 {
 	//! create a special cylindrical-engine
-	m_Engine = Engine_Cylinder::New(this, m_numThreads);
+	m_Engine = Engine_Cylinder::New(this, m_orig_numThreads);
 	return m_Engine;
 }
 

FDTD/operator_cylinder.h

@@ -48,14 +48,14 @@ public:
 	//! Get the coordinates for a given node index and component, according to the cylindrical yee-algorithm. Returns true if inside the FDTD domain.
 	virtual bool GetYeeCoords(int ny, unsigned int pos[3], double* coords, bool dualMesh) const;
 
-	//! Get the node width for a given direction \a n and a given mesh posisition \a pos
+	//! Get the node width for a given direction \a n and a given mesh position \a pos
 	virtual double GetNodeWidth(int ny, const unsigned int pos[3], bool dualMesh = false) const;
-	//! Get the node width for a given direction \a n and a given mesh posisition \a pos
+	//! Get the node width for a given direction \a n and a given mesh position \a pos
 	virtual double GetNodeWidth(int ny, const int pos[3], bool dualMesh = false) const;
 
-	//! Get the node area for a given direction \a n and a given mesh posisition \a pos
+	//! Get the node area for a given direction \a n and a given mesh position \a pos
 	virtual double GetNodeArea(int n, const unsigned int* pos, bool dualMesh=false) const;
-	//! Get the node area for a given direction \a n and a given mesh posisition \a pos
+	//! Get the node area for a given direction \a n and a given mesh position \a pos
 	virtual double GetNodeArea(int ny, const int pos[3], bool dualMesh = false) const;
 
 	//! Get the length of an FDTD edge, including radius corrected alpha-mesh width.
@@ -64,7 +64,7 @@ public:
 	//! Get the volume of an FDTD cell
 	virtual double GetCellVolume(const unsigned int pos[3], bool dualMesh = false) const;
 
-	//! Get the area around an edge for a given direction \a n and a given mesh posisition \a pos
+	//! Get the area around an edge for a given direction \a n and a given mesh position \a pos
 	/*!
 		This will return the area around an edge with a given direction, measured at the middle of the edge.
 		In a cartesian mesh this is equal to the NodeArea, may be different in other coordinate systems.

FDTD/operator_cylindermultigrid.cpp

@@ -51,7 +51,7 @@ Operator_CylinderMultiGrid* Operator_CylinderMultiGrid::New(vector<double> Split
 
 Engine* Operator_CylinderMultiGrid::CreateEngine()
 {
-	m_Engine = Engine_CylinderMultiGrid::New(this,m_numThreads);
+	m_Engine = Engine_CylinderMultiGrid::New(this, m_orig_numThreads);
 	return m_Engine;
 }
 

FDTD/operator_mpi.cpp

@@ -53,7 +53,7 @@ double Operator_MPI::CalcTimestep()
 		return ret;
 
 	double local_dT = dT;
-	//find the smallest time-step requestes by all processings
+	//find the smallest time-step requests by all processings
 	MPI_Reduce(&local_dT, &dT, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD);
 	//send the smallest time-step to all
 	MPI_Bcast(&dT, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

FDTD/operator_multithread.cpp

@@ -36,11 +36,12 @@ Operator_Multithread::~Operator_Multithread()
 void Operator_Multithread::setNumThreads( unsigned int numThreads )
 {
 	m_numThreads = numThreads;
+	m_orig_numThreads = numThreads;
 }
 
 Engine* Operator_Multithread::CreateEngine()
 {
-	m_Engine = Engine_Multithread::New(this,m_numThreads);
+	m_Engine = Engine_Multithread::New(this, m_orig_numThreads);
 	return m_Engine;
 }
 
@@ -106,7 +107,7 @@ void Operator_Multithread::CalcStartStopLines(unsigned int &numThreads, vector<u
 
 int Operator_Multithread::CalcECOperator( DebugFlags debugFlags )
 {
-	if (m_numThreads == 0)
+	if ((m_numThreads == 0) || (m_numThreads > boost::thread::hardware_concurrency()))
 		m_numThreads = boost::thread::hardware_concurrency();
 
 	vector<unsigned int> m_Start_Lines;

FDTD/operator_multithread.h

@@ -64,6 +64,7 @@ protected:
 
 	boost::thread_group m_thread_group;
 	unsigned int m_numThreads; // number of worker threads
+	unsigned int m_orig_numThreads;
 
 	//! Calculate the start/stop lines for the multithreading operator and engine.
 	/*!

FDTD/operator_sse.cpp

@@ -41,11 +41,11 @@ Operator_sse::~Operator_sse()
 	Delete();
 }
 
-Engine* Operator_sse::CreateEngine() const
+Engine* Operator_sse::CreateEngine()
 {
 	//! create a special sse-engine
-	Engine_sse* eng = Engine_sse::New(this);
-	return eng;
+	m_Engine = Engine_sse::New(this);
+	return m_Engine;
 }
 
 void Operator_sse::Init()

FDTD/operator_sse.h

@@ -29,7 +29,7 @@ public:
 	static Operator_sse* New();
 	virtual ~Operator_sse();
 
-	virtual Engine* CreateEngine() const;
+	virtual Engine* CreateEngine();
 
 	inline virtual FDTD_FLOAT GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_vv[n][x][y][z%numVectors].f[z/numVectors]; }
 	inline virtual FDTD_FLOAT GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_vi[n][x][y][z%numVectors].f[z/numVectors]; }

FDTD/operator_sse_compressed.cpp

@@ -42,16 +42,14 @@ Operator_SSE_Compressed::~Operator_SSE_Compressed()
 	Delete();
 }
 
-Engine* Operator_SSE_Compressed::CreateEngine() const
+Engine* Operator_SSE_Compressed::CreateEngine()
 {
 	if (!m_Use_Compression)
-	{
 		//! create a default sse-engine
-		Engine_sse* eng = Engine_sse::New(this);
-		return eng;
-	}
-	Engine_SSE_Compressed* eng = Engine_SSE_Compressed::New(this);
-	return eng;
+		m_Engine = Engine_sse::New(this);
+  else
+	  m_Engine = Engine_SSE_Compressed::New(this);
+	return m_Engine;
 }
 
 int Operator_SSE_Compressed::CalcECOperator( DebugFlags debugFlags )

FDTD/operator_sse_compressed.h

@@ -43,7 +43,7 @@ public:
 	static Operator_SSE_Compressed* New();
 	virtual ~Operator_SSE_Compressed();
 
-	virtual Engine* CreateEngine() const;
+	virtual Engine* CreateEngine();
 
 	inline virtual FDTD_FLOAT GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_vv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetVV(n,x,y,z);}
 	inline virtual FDTD_FLOAT GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_vi_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetVI(n,x,y,z);}

INSTALL

@@ -22,6 +22,6 @@ Install instructions for openEMS:
     make install (may require root)
 
     Note:
-      - all path informations may be stored in a localConfig.cmake
+      - all path information may be stored in a localConfig.cmake
       - the default "prefix" is /usr/local
 

TODO

@@ -3,10 +3,10 @@
 
 todo / mandatory for v0.1.0:
 - more examples and lots of testing...
-- improvments and testing for MPI engine
+- improvements and testing for MPI engine
 
 wishes:
-- location dependend excitation given by a file (e.g. mode-profile simulated with comsol)
+- location dependent excitation given by a file (e.g. mode-profile simulated with comsol)
 - more import filter (e.g. gerber import)
 
 outlook:

linux/debian.changelog

@@ -17,7 +17,7 @@ openems (0.0.25-4) stable; urgency=low
   * Upstream fixes
  -- Sebastian Held <sebastian.held@gmx.de>  Sun, 25 Dec 2011 21:23:38 +0100
 openems (0.0.25-3) stable; urgency=low
-  * Changed dependancy on boost to enable build on oneiric
+  * Changed dependency on boost to enable build on oneiric
  -- Sebastian Held <sebastian.held@gmx.de>  Sun, 18 Dec 2011 21:39:38 +0100
 openems (0.0.25-2) stable; urgency=low
   * New upstream release

main.cpp

@@ -68,8 +68,12 @@ int main(int argc, char *argv[])
 	}
 
 	int EC = FDTD.ParseFDTDSetup(argv[1]);
+	if(!EC) {
+	  cerr << "openEMS - ParseFDTDSetup failed." << endl;
+	  exit(1);
+	}
 	EC = FDTD.SetupFDTD();
-	if (EC) return EC;
+	if (EC) exit(EC);
 	FDTD.RunFDTD();
 
 #ifdef MPI_SUPPORT
@@ -77,5 +81,5 @@ int main(int argc, char *argv[])
 	MPI::Finalize();
 #endif
 
-	return 0;
+	exit(0);
 }

matlab/AR_estimate.m

@@ -20,7 +20,7 @@ function [val_ar t_ar f_val_ar EC] = AR_estimate( t, val, freq, nu, mu, expand_f
 %           0 --> no error
 %           1 --> input error: t and val mismatch
 %           2 --> input error: mu has to be larger than 2*nu
-%           3 --> inout error: expand_factor has to be larger than 1
+%           3 --> input error: expand_factor has to be larger than 1
 %           10 --> AR error: signal is to short for AR estimate --> decrease AR order
 %           11 --> AR error: estimated signal appears to be unstable --> use a different mu
 %

matlab/Add2Queue.m

@@ -1,7 +1,7 @@
 function [queue] = Add2Queue(queue,func_name, func_args, varargin)
 % function [queue] = Add2Queue(queue,func_name, func_args, varargin)
 %
-% Use this function to add a funtion to the queue.
+% Use this function to add a function to the queue.
 %
 % For more details see: InitQueue
 %

matlab/AddCPWPort.m

@@ -17,7 +17,7 @@ function [CSX,port] = AddCPWPort( CSX, prio, portnr, materialname, start, stop,
 %               is false)
 % 'FeedShift'   shift to port from start by a given distance in drawing
 %               units. Default is 0. Only active if 'ExcitePort' is set!
-% 'Feed_R'      Specifiy a lumped port resistance. Default is no lumped
+% 'Feed_R'      Specify a lumped port resistance. Default is no lumped
 %               port resistance --> port has to end in an ABC.
 % 'MeasPlaneShift'  Shift the measurement plane from start t a given distance
 %               in drawing units. Default is the middle of start/stop.
@@ -68,7 +68,7 @@ evec0 = evec ./ sum(evec); % evec0 is a unit vector
 
 %set defaults
 feed_shift = 0;
-feed_R = inf; %(default is open, no resitance)
+feed_R = inf; %(default is open, no resistance)
 excite = false;
 measplanepos = nan;
 PortNamePrefix = '';
@@ -84,12 +84,12 @@ for n=1:2:numel(varargin)
         end
     elseif (strcmp(varargin{n},'Feed_R')==1);
         feed_R = varargin{n+1};
-        if (numel(feed_shift)>1)
+        if (numel(feed_R)>1)
             error 'Feed_R must be a scalar value'
         end
     elseif (strcmp(varargin{n},'MeasPlaneShift')==1);
         measplanepos = varargin{n+1};
-        if (numel(feed_shift)>1)
+        if (numel(measplanepos)>1)
             error 'MeasPlaneShift must be a scalar value'
         end
     elseif (strcmp(varargin{n},'ExcitePort')==1);
@@ -152,28 +152,32 @@ else
 end
 
 % calculate position of the voltage probes
-mesh{1} = sort(CSX.RectilinearGrid.XLines);
-mesh{2} = sort(CSX.RectilinearGrid.YLines);
-mesh{3} = sort(CSX.RectilinearGrid.ZLines);
-meshlines = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), measplanepos, 'nearest' );
-meshlines = mesh{idx_prop}(meshlines-1:meshlines+1); % get three lines (approx. at center)
-if direction == -1
+try
+	mesh{1} = sort(CSX.RectilinearGrid.XLines);
+	mesh{2} = sort(CSX.RectilinearGrid.YLines);
+	mesh{3} = sort(CSX.RectilinearGrid.ZLines);
+	meshlines = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), measplanepos, 'nearest' );
+	meshlines = mesh{idx_prop}(meshlines-1:meshlines+1); % get three lines (approx. at center)
+	if direction == -1
 		meshlines = fliplr(meshlines);
+	end
+	SL_w2 = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), (nstart(idx_width)+nstop(idx_width))/2, 'nearest' );
+	SL_w2 = mesh{idx_width}(SL_w2); % get e-line at center of CPW (SL_width/2)
+	v1_start(idx_prop)   = meshlines(1);
+	v1_start(idx_width)  = (nstart(idx_width)+nstop(idx_width))/2;
+	v1_start(idx_height) = start(idx_height);
+	v1_stop  = v1_start;
+	v2_start = v1_start;
+	v2_stop  = v1_stop;
+	v2_start(idx_prop)   = meshlines(2);
+	v2_stop(idx_prop)    = meshlines(2);
+	v3_start = v2_start;
+	v3_stop  = v2_stop;
+	v3_start(idx_prop)   = meshlines(3);
+	v3_stop(idx_prop)    = meshlines(3);
+catch
+	error('Unable to place voltage probe on mesh; check the location of the CPW and the probe (MeasPlaneShift), and make sure that the mesh is large enough');
 end
-SL_w2 = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), (nstart(idx_width)+nstop(idx_width))/2, 'nearest' );
-SL_w2 = mesh{idx_width}(SL_w2); % get e-line at center of CPW (SL_width/2)
-v1_start(idx_prop)   = meshlines(1);
-v1_start(idx_width)  = (nstart(idx_width)+nstop(idx_width))/2;
-v1_start(idx_height) = start(idx_height);
-v1_stop  = v1_start;
-v2_start = v1_start;
-v2_stop  = v1_stop;
-v2_start(idx_prop)   = meshlines(2);
-v2_stop(idx_prop)    = meshlines(2);
-v3_start = v2_start;
-v3_stop  = v2_stop;
-v3_start(idx_prop)   = meshlines(3);
-v3_stop(idx_prop)    = meshlines(3);
 
 width_add_start = [0 0 0];
 width_add_stop  = [0 0 0];
@@ -183,7 +187,7 @@ width_add_stop(idx_width)  = (nstop(idx_width)-nstart(idx_width))/2+gap_width;
 weight = 0.5;
 % create the voltage-probes
 port.U_filename{1,1} = [PortNamePrefix 'port_ut' num2str(portnr) 'A1'];
-CSX = AddProbe( CSX, port.U_filename{1,1}, 0, 'weight', -1*weight );
+CSX = AddProbe( CSX, port.U_filename{1,1}, 0, 'weight', weight );
 CSX = AddBox( CSX, port.U_filename{1,1}, prio, v1_start-width_add_start, v1_stop-width_add_stop);
 
 port.U_filename{1,2} = [PortNamePrefix 'port_ut' num2str(portnr) 'A2'];
@@ -192,7 +196,7 @@ CSX = AddBox( CSX, port.U_filename{1,2}, prio, v1_start+width_add_start, v1_stop
 
 
 port.U_filename{2,1} = [PortNamePrefix 'port_ut' num2str(portnr) 'B1'];
-CSX = AddProbe( CSX, port.U_filename{2,1}, 0, 'weight', -1*weight );
+CSX = AddProbe( CSX, port.U_filename{2,1}, 0, 'weight', weight );
 CSX = AddBox( CSX, port.U_filename{2,1}, prio, v2_start-width_add_start, v2_stop-width_add_stop );
 
 port.U_filename{2,2} = [PortNamePrefix 'port_ut' num2str(portnr) 'B2'];
@@ -201,7 +205,7 @@ CSX = AddBox( CSX, port.U_filename{2,2}, prio, v2_start+width_add_start, v2_stop
 
 
 port.U_filename{3,1} = [PortNamePrefix 'port_ut' num2str(portnr) 'C1'];
-CSX = AddProbe( CSX, port.U_filename{3,1}, 0, 'weight', -1*weight );
+CSX = AddProbe( CSX, port.U_filename{3,1}, 0, 'weight', weight );
 CSX = AddBox( CSX, port.U_filename{3,1}, prio, v3_start-width_add_start, v3_stop-width_add_stop );
 
 port.U_filename{3,2} = [PortNamePrefix 'port_ut' num2str(portnr) 'C2'];
@@ -209,19 +213,23 @@ CSX = AddProbe( CSX, port.U_filename{3,2}, 0, 'weight', weight );
 CSX = AddBox( CSX, port.U_filename{3,2}, prio, v3_start+width_add_start, v3_stop+width_add_stop );
 
 % calculate position of the current probes
-idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstart(idx_width), 'nearest' );
-i1_start(idx_width)  = mesh{idx_width}(idx) - diff(mesh{idx_width}(idx-1:idx))/2;
-idx = interp1( mesh{idx_height}, 1:numel(mesh{idx_height}), start(idx_height), 'nearest' );
-i1_start(idx_height) = mesh{idx_height}(idx-1) - diff(mesh{idx_height}(idx-2:idx-1))/2;
-i1_stop(idx_height)  = mesh{idx_height}(idx+1) + diff(mesh{idx_height}(idx+1:idx+2))/2;
-i1_start(idx_prop)   = sum(meshlines(1:2))/2;
-i1_stop(idx_prop)    = i1_start(idx_prop);
-idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstop(idx_width), 'nearest' );
-i1_stop(idx_width)   = mesh{idx_width}(idx) + diff(mesh{idx_width}(idx:idx+1))/2;
-i2_start = i1_start;
-i2_stop  = i1_stop;
-i2_start(idx_prop)   = sum(meshlines(2:3))/2;
-i2_stop(idx_prop)    = i2_start(idx_prop);
+try
+	idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstart(idx_width), 'nearest' );
+	i1_start(idx_width)  = mesh{idx_width}(idx) - diff(mesh{idx_width}(idx-1:idx))/2;
+	idx = interp1( mesh{idx_height}, 1:numel(mesh{idx_height}), start(idx_height), 'nearest' );
+	i1_start(idx_height) = mesh{idx_height}(idx-1) - diff(mesh{idx_height}(idx-2:idx-1))/2;
+	i1_stop(idx_height)  = mesh{idx_height}(idx+1) + diff(mesh{idx_height}(idx+1:idx+2))/2;
+	i1_start(idx_prop)   = sum(meshlines(1:2))/2;
+	i1_stop(idx_prop)    = i1_start(idx_prop);
+	idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstop(idx_width), 'nearest' );
+	i1_stop(idx_width)   = mesh{idx_width}(idx) + diff(mesh{idx_width}(idx:idx+1))/2;
+	i2_start = i1_start;
+	i2_stop  = i1_stop;
+	i2_start(idx_prop)   = sum(meshlines(2:3))/2;
+	i2_stop(idx_prop)    = i2_start(idx_prop);
+catch
+	error('Unable to place current probe on mesh; check the location of the CPW and the probe (MeasPlaneShift), and make sure that the mesh is large enough');
+end
 
 % create the curr-probes
 weight = direction;
@@ -248,13 +256,17 @@ port.measplanepos = abs(v2_start(idx_prop) - start(idx_prop))*port.LengthScale;
 % port
 
 % create excitation (if enabled) and port resistance
-meshline = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), start(idx_prop) + feed_shift*direction, 'nearest' );
-ex_start(idx_prop)   = mesh{idx_prop}(meshline) ;
-ex_start(idx_width)  = (nstart(idx_width)+nstop(idx_width))/2;
-ex_start(idx_height) = nstart(idx_height);
-ex_stop(idx_prop)    = ex_start(idx_prop);
-ex_stop(idx_width)   = (nstart(idx_width)+nstop(idx_width))/2;
-ex_stop(idx_height)  = nstop(idx_height);
+try
+	meshline = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), start(idx_prop) + feed_shift*direction, 'nearest' );
+	ex_start(idx_prop)   = mesh{idx_prop}(meshline) ;
+	ex_start(idx_width)  = (nstart(idx_width)+nstop(idx_width))/2;
+	ex_start(idx_height) = nstart(idx_height);
+	ex_stop(idx_prop)    = ex_start(idx_prop);
+	ex_stop(idx_width)   = (nstart(idx_width)+nstop(idx_width))/2;
+	ex_stop(idx_height)  = nstop(idx_height);
+catch
+	error('Unable to place excitation on mesh; check the location of the CPW and the excitation (FeedShift), and make sure that the mesh is large enough');
+end
 
 port.excite = 0;
 if excite
@@ -265,7 +277,7 @@ if excite
     CSX = AddBox( CSX, [PortNamePrefix 'port_excite_2_' num2str(portnr)], prio, ex_start+width_add_start, ex_stop+width_add_stop );
 end
 
-%% CPW resitance at start of CPW line
+%% CPW resistance at start of CPW line
 ex_start(idx_prop) = start(idx_prop);
 ex_stop(idx_prop) = ex_start(idx_prop);
 
@@ -280,6 +292,6 @@ elseif feed_R == 0
     CSX = AddBox( CSX, materialname, prio, ex_start-width_add_start, ex_stop-width_add_stop );
     CSX = AddBox( CSX, materialname, prio, ex_start+width_add_start, ex_stop+width_add_stop );
 else
-    error('openEMS:AddCPWPort','CPW port with resitance <= 0 it not possible');
+    error('openEMS:AddCPWPort','CPW port with resistance <= 0 it not possible');
 end
 end

matlab/AddCoaxialPort.m

@@ -20,7 +20,7 @@ function [CSX,port] = AddCoaxialPort( CSX, prio, portnr, pec_name, materialname,
 %               set to 0 (default) for a passive port
 % 'FeedShift'   shift to port from start by a given distance in drawing
 %               units. Default is 0. Only active if 'ExciteAmp' is set!
-% 'Feed_R'      Specifiy a lumped port resistance. Default is no lumped
+% 'Feed_R'      Specify a lumped port resistance. Default is no lumped
 %               port resistance --> port has to end in an ABC.
 % 'MeasPlaneShift'  Shift the measurement plane from start t a given distance
 %               in drawing units. Default is the middle of start/stop.
@@ -50,7 +50,7 @@ dir = DirChar2Int(dir);
 
 %set defaults
 feed_shift = 0;
-feed_R = inf; %(default is open, no resitance)
+feed_R = inf; %(default is open, no resistance)
 excite_amp = 0;
 measplanepos = nan;
 PortNamePrefix = '';
@@ -66,12 +66,12 @@ for n=1:2:numel(varargin)
         end
     elseif (strcmp(varargin{n},'Feed_R')==1);
         feed_R = varargin{n+1};
-        if (numel(feed_shift)>1)
+        if (numel(feed_R)>1)
             error 'Feed_R must be a scalar value'
         end
     elseif (strcmp(varargin{n},'MeasPlaneShift')==1);
         measplanepos = varargin{n+1};
-        if (numel(feed_shift)>1)
+        if (numel(measplanepos)>1)
             error 'MeasPlaneShift must be a scalar value'
         end
     elseif (strcmp(varargin{n},'ExciteAmp')==1);
@@ -114,27 +114,31 @@ else
 end
 
 % calculate position of the voltage probes
-mesh{1} = sort(unique(CSX.RectilinearGrid.XLines));
-mesh{2} = sort(unique(CSX.RectilinearGrid.YLines));
-mesh{3} = sort(unique(CSX.RectilinearGrid.ZLines));
-meshlines = interp1( mesh{idx_prop_n}, 1:numel(mesh{idx_prop_n}), measplanepos, 'nearest' );
-meshlines = mesh{idx_prop_n}(meshlines-1:meshlines+1); % get three lines (approx. at center)
-if direction == -1
+try
+	mesh{1} = sort(unique(CSX.RectilinearGrid.XLines));
+	mesh{2} = sort(unique(CSX.RectilinearGrid.YLines));
+	mesh{3} = sort(unique(CSX.RectilinearGrid.ZLines));
+	meshlines = interp1( mesh{idx_prop_n}, 1:numel(mesh{idx_prop_n}), measplanepos, 'nearest' );
+	meshlines = mesh{idx_prop_n}(meshlines-1:meshlines+1); % get three lines (approx. at center)
+	if direction == -1
 		meshlines = fliplr(meshlines);
+	end
+	v1_start(idx_prop_n)   = meshlines(1);
+	v1_start(idx_prop_nP)  = start(idx_prop_nP)+r_i;
+	v1_start(idx_prop_nPP) = start(idx_prop_nPP);
+	v1_stop  = v1_start;
+	v1_stop(idx_prop_nP)  = start(idx_prop_nP)+r_o;
+	v2_start = v1_start;
+	v2_stop  = v1_stop;
+	v2_start(idx_prop_n)   = meshlines(2);
+	v2_stop(idx_prop_n)    = meshlines(2);
+	v3_start = v2_start;
+	v3_stop  = v2_stop;
+	v3_start(idx_prop_n)   = meshlines(3);
+	v3_stop(idx_prop_n)    = meshlines(3);
+catch
+	error('Unable to place voltage probe on mesh; check the location of the port and the probe (MeasPlaneShift), and make sure that the mesh is large enough');
 end
-v1_start(idx_prop_n)   = meshlines(1);
-v1_start(idx_prop_nP)  = start(idx_prop_nP)+r_i;
-v1_start(idx_prop_nPP) = start(idx_prop_nPP);
-v1_stop  = v1_start;
-v1_stop(idx_prop_nP)  = start(idx_prop_nP)+r_o;
-v2_start = v1_start;
-v2_stop  = v1_stop;
-v2_start(idx_prop_n)   = meshlines(2);
-v2_stop(idx_prop_n)    = meshlines(2);
-v3_start = v2_start;
-v3_stop  = v2_stop;
-v3_start(idx_prop_n)   = meshlines(3);
-v3_stop(idx_prop_n)    = meshlines(3);
 
 % calculate position of the current probes
 i1_start(idx_prop_n)   = 0.5*(meshlines(1)+meshlines(2));
@@ -184,12 +188,16 @@ port.r_i = r_i;
 port.r_o = r_o;
 
 % create excitation (if enabled) and port resistance
-meshline = interp1( mesh{idx_prop_n}, 1:numel(mesh{idx_prop_n}), start(idx_prop_n) + feed_shift*direction, 'nearest' );
-min_cell_prop = min(diff(mesh{idx_prop_n}));
-ex_start = start;
-ex_start(idx_prop_n)   = mesh{idx_prop_n}(meshline) - 0.01*min_cell_prop;
-ex_stop  = ex_start;
-ex_stop(idx_prop_n)   = mesh{idx_prop_n}(meshline) + 0.01*min_cell_prop;
+try
+	meshline = interp1( mesh{idx_prop_n}, 1:numel(mesh{idx_prop_n}), start(idx_prop_n) + feed_shift*direction, 'nearest' );
+	min_cell_prop = min(diff(mesh{idx_prop_n}));
+	ex_start = start;
+	ex_start(idx_prop_n)   = mesh{idx_prop_n}(meshline) - 0.01*min_cell_prop;
+	ex_stop  = ex_start;
+	ex_stop(idx_prop_n)   = mesh{idx_prop_n}(meshline) + 0.01*min_cell_prop;
+catch
+	error('Unable to place excitation on mesh; check the location of the port and the excitation (FeedShift), and make sure that the mesh is large enough');
+end
 
 port.excite = 0;
 if (excite_amp~=0)
@@ -213,7 +221,7 @@ if (excite_amp~=0)
     CSX = AddCylindricalShell(CSX,[PortNamePrefix 'port_excite_' num2str(portnr)],0 ,ex_start,ex_stop,0.5*(r_i+r_o),(r_o-r_i));
 end
 
-%% resitance at start of coaxial line
+%% resistance at start of coaxial line
 ex_start = start;
 ex_stop = stop;
 ex_stop(idx_prop_n) = ex_start(idx_prop_n);
@@ -227,6 +235,6 @@ elseif feed_R == 0
     CSX = AddBox( CSX, pec_name, prio, ex_start, ex_stop );
     CSX = AddCylindricalShell(CSX, pec_name, prio ,ex_start, ex_stop, 0.5*(r_i+r_o),(r_o-r_i));
 else
-    error('openEMS:AddMSLPort','MSL port with resitance <= 0 it not possible');
+    error('openEMS:AddCoaxialPort','Coaxial port with resistance <= 0 it not possible');
 end
 end

matlab/AddCurvePort.m

@@ -78,9 +78,10 @@ for n=1:3
 end
 
 % calculate position
-port_start_idx = start_idx;
-port_stop_idx  = stop_idx;
-if abs(start_idx(dir) - stop_idx(dir)) ~= 1
+try
+	port_start_idx = start_idx;
+	port_stop_idx  = stop_idx;
+	if abs(start_idx(dir) - stop_idx(dir)) ~= 1
 		% calc port position
 		idx  = interp1( mesh{dir},  1:numel(mesh{dir}),  (nstart(dir)+nstop(dir))/2,   'nearest' );
 		idx1 = interp1( mesh{dir1}, 1:numel(mesh{dir1}), (nstart(dir1)+nstop(dir1))/2, 'nearest' );
@@ -95,6 +96,9 @@ if abs(start_idx(dir) - stop_idx(dir)) ~= 1
 		CSX = AddMetal( CSX, metalname );
 		CSX = AddCurve( CSX, metalname, prio, [nstart.' [mesh{1}(port_start_idx(1));mesh{2}(port_start_idx(2));mesh{3}(port_start_idx(3))]] );
 		CSX = AddCurve( CSX, metalname, prio, [nstop.' [mesh{1}(port_stop_idx(1));mesh{2}(port_stop_idx(2));mesh{3}(port_stop_idx(3))]] );
+	end
+catch
+	error('Unable to place port on mesh; check the location of the port, and make sure that the mesh is large enough');
 end
 
 % calculate position of resistive material

matlab/AddLumpedPort.m

@@ -15,8 +15,8 @@ function [CSX, port] = AddLumpedPort( CSX, prio, portnr, R, start, stop, dir, ex
 % dir:      direction/amplitude of port (e.g.: [1 0 0], [0 1 0] or [0 0 1])
 % excite (optional): if true, the port will be switched on (see AddExcitation())
 %                       Note: for legacy support a string will be accepted
-% V_Probe_Weight:  additional weigth for the voltage probes
-% I_Probe_Weight:  additional weigth for the current probes
+% V_Probe_Weight:  additional weight for the voltage probes
+% I_Probe_Weight:  additional weight for the current probes
 % optional (key/values):
 %   'PortNamePrefix': an prefix to the port name
 % varargin (optional): additional excitations options, see also AddExcitation
@@ -25,7 +25,7 @@ function [CSX, port] = AddLumpedPort( CSX, prio, portnr, R, start, stop, dir, ex
 %   start = [0 -width/2 0];
 %   stop  = [0  width/2 height];
 %   [CSX] = AddLumpedPort(CSX, 5 ,1 , 50, start, stop, [0 0 1], true);
-%   %this defines an active lumped port in z-direction with a 50 Ohm port impedence
+%   %this defines an active lumped port in z-direction with a 50 Ohm port impedance
 %
 % openEMS matlab interface
 % -----------------------
@@ -94,7 +94,7 @@ end
 
 % legacy support, will be removed at some point
 if ischar(excite)
-    warning('CSXCAD:AddLumpedPort','depreceated: a string as excite option is no longer supported and will be removed in the future, please use true or false');
+    warning('CSXCAD:AddLumpedPort','deprecated: a string as excite option is no longer supported and will be removed in the future, please use true or false');
     if ~isempty(excite)
         excite = true;
     else
@@ -115,7 +115,7 @@ u_start(n_dir) = start(n_dir);
 u_stop(n_dir)  = stop(n_dir);
 
 port.U_filename = [PortNamePrefix 'port_ut' int2str(portnr)];
-CSX = AddProbe(CSX, port.U_filename, 0, 'weight', -direction*V_Probe_Weight);
+CSX = AddProbe(CSX, port.U_filename, 0, 'weight', -1*V_Probe_Weight);
 CSX = AddBox(CSX, port.U_filename, prio, u_start, u_stop);
 
 i_start = start;

matlab/AddMRStub.m

@@ -70,4 +70,4 @@ p(1,:) = p(1,:) + position(dim1);
 p(2,:) = p(2,:) + position(dim2);
 
 elevation = position(idx_elevation);
-CSX = AddPolygon( CSX, materialname, prio, normVector, elevation, p );
+CSX = AddPolygon( CSX, materialname, prio, idx_elevation-1, elevation, p );

matlab/AddMSLPort.m

@@ -16,7 +16,7 @@ function [CSX,port] = AddMSLPort( CSX, prio, portnr, materialname, start, stop,
 %               is false)
 % 'FeedShift'   shift to port from start by a given distance in drawing
 %               units. Default is 0. Only active if 'ExcitePort' is set!
-% 'Feed_R'      Specifiy a lumped port resistance. Default is no lumped
+% 'Feed_R'      Specify a lumped port resistance. Default is no lumped
 %               port resistance --> port has to end in an ABC.
 % 'MeasPlaneShift'  Shift the measurement plane from start t a given distance 
 %               in drawing units. Default is the middle of start/stop.
@@ -37,7 +37,7 @@ function [CSX,port] = AddMSLPort( CSX, prio, portnr, materialname, start, stop,
 %   - the excitation is active and placed at x=start(1) ('ExcitePort', true)
 %   - a 50 Ohm lumped port resistance is placed at x=start(1) ('Feed_R', 50)
 %   - the width-direction is determined by the cross product of the
-%       direction of propagtion (dir='x') and the excitation vector
+%       direction of propagation (dir='x') and the excitation vector
 %       (evec=[0 0 -1]), in this case it is the y-direction
 %   - the MSL-metal is created in a xy-plane at a height at z=start(3)
 %     --> It is important to define the MSL height in the start coordinate!
@@ -69,7 +69,7 @@ evec0 = evec ./ sum(evec); % evec0 is a unit vector
 
 %set defaults
 feed_shift = 0;
-feed_R = inf; %(default is open, no resitance)
+feed_R = inf; %(default is open, no resistance)
 excite = false;
 measplanepos = nan;
 PortNamePrefix = '';
@@ -85,17 +85,17 @@ for n=1:2:numel(varargin)
         end
     elseif (strcmp(varargin{n},'Feed_R')==1);
         feed_R = varargin{n+1};
-        if (numel(feed_shift)>1)
+        if (numel(feed_R)>1)
             error 'Feed_R must be a scalar value'
         end
     elseif (strcmp(varargin{n},'MeasPlaneShift')==1);
         measplanepos = varargin{n+1};
-        if (numel(feed_shift)>1)
+        if (numel(measplanepos)>1)
             error 'MeasPlaneShift must be a scalar value'
         end
     elseif (strcmp(varargin{n},'ExcitePort')==1);
         if ischar(varargin{n+1})
-            warning('CSXCAD:AddMSLPort','depreceated: a string as excite option is no longer supported and will be removed in the future, please use true or false');
+            warning('CSXCAD:AddMSLPort','deprecated: a string as excite option is no longer supported and will be removed in the future, please use true or false');
             if ~isempty(excite)
                 excite = true;
             else
@@ -136,13 +136,6 @@ else
     direction = -1;
 end
 
-% direction of propagation
-if stop(idx_height)-start(idx_height) > 0
-    upsidedown = +1;
-else
-    upsidedown = -1;
-end
-
 % create the metal/material for the MSL
 MSL_start = start;
 MSL_stop = stop;
@@ -156,49 +149,57 @@ else
 end
 
 % calculate position of the voltage probes
-mesh{1} = sort(CSX.RectilinearGrid.XLines);
-mesh{2} = sort(CSX.RectilinearGrid.YLines);
-mesh{3} = sort(CSX.RectilinearGrid.ZLines);
-meshlines = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), measplanepos, 'nearest' );
-meshlines = mesh{idx_prop}(meshlines-1:meshlines+1); % get three lines (approx. at center)
-if direction == -1
+try
+	mesh{1} = sort(CSX.RectilinearGrid.XLines);
+	mesh{2} = sort(CSX.RectilinearGrid.YLines);
+	mesh{3} = sort(CSX.RectilinearGrid.ZLines);
+	meshlines = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), measplanepos, 'nearest' );
+	meshlines = mesh{idx_prop}(meshlines-1:meshlines+1); % get three lines (approx. at center)
+	if direction == -1
 		meshlines = fliplr(meshlines);
+	end
+	MSL_w2 = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), (nstart(idx_width)+nstop(idx_width))/2, 'nearest' );
+	MSL_w2 = mesh{idx_width}(MSL_w2); % get e-line at center of MSL (MSL_width/2)
+	v1_start(idx_prop)   = meshlines(1);
+	v1_start(idx_width)  = MSL_w2;
+	v1_start(idx_height) = start(idx_height);
+	v1_stop  = v1_start;
+	v1_stop(idx_height)  = stop(idx_height);
+	v2_start = v1_start;
+	v2_stop  = v1_stop;
+	v2_start(idx_prop)   = meshlines(2);
+	v2_stop(idx_prop)    = meshlines(2);
+	v3_start = v2_start;
+	v3_stop  = v2_stop;
+	v3_start(idx_prop)   = meshlines(3);
+	v3_stop(idx_prop)    = meshlines(3);
+catch
+	error('Unable to place voltage probe on mesh; check the location of the MSL and the probe (MeasPlaneShift), and make sure that the mesh is large enough');
 end
-MSL_w2 = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), (nstart(idx_width)+nstop(idx_width))/2, 'nearest' );
-MSL_w2 = mesh{idx_width}(MSL_w2); % get e-line at center of MSL (MSL_width/2)
-v1_start(idx_prop)   = meshlines(1);
-v1_start(idx_width)  = MSL_w2;
-v1_start(idx_height) = start(idx_height);
-v1_stop  = v1_start;
-v1_stop(idx_height)  = stop(idx_height);
-v2_start = v1_start;
-v2_stop  = v1_stop;
-v2_start(idx_prop)   = meshlines(2);
-v2_stop(idx_prop)    = meshlines(2);
-v3_start = v2_start;
-v3_stop  = v2_stop;
-v3_start(idx_prop)   = meshlines(3);
-v3_stop(idx_prop)    = meshlines(3);
 
 % calculate position of the current probes
-idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstart(idx_width), 'nearest' );
-i1_start(idx_width)  = mesh{idx_width}(idx) - diff(mesh{idx_width}(idx-1:idx))/2;
-idx = interp1( mesh{idx_height}, 1:numel(mesh{idx_height}), start(idx_height), 'nearest' );
-i1_start(idx_height) = mesh{idx_height}(idx-1) - diff(mesh{idx_height}(idx-2:idx-1))/2;
-i1_stop(idx_height)  = mesh{idx_height}(idx+1) + diff(mesh{idx_height}(idx+1:idx+2))/2;
-i1_start(idx_prop)   = sum(meshlines(1:2))/2;
-i1_stop(idx_prop)    = i1_start(idx_prop);
-idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstop(idx_width), 'nearest' );
-i1_stop(idx_width)   = mesh{idx_width}(idx) + diff(mesh{idx_width}(idx:idx+1))/2;
-i2_start = i1_start;
-i2_stop  = i1_stop;
-i2_start(idx_prop)   = sum(meshlines(2:3))/2;
-i2_stop(idx_prop)    = i2_start(idx_prop);
+try
+	idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstart(idx_width), 'nearest' );
+	i1_start(idx_width)  = mesh{idx_width}(idx) - diff(mesh{idx_width}(idx-1:idx))/2;
+	idx = interp1( mesh{idx_height}, 1:numel(mesh{idx_height}), start(idx_height), 'nearest' );
+	i1_start(idx_height) = mesh{idx_height}(idx-1) - diff(mesh{idx_height}(idx-2:idx-1))/2;
+	i1_stop(idx_height)  = mesh{idx_height}(idx+1) + diff(mesh{idx_height}(idx+1:idx+2))/2;
+	i1_start(idx_prop)   = sum(meshlines(1:2))/2;
+	i1_stop(idx_prop)    = i1_start(idx_prop);
+	idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstop(idx_width), 'nearest' );
+	i1_stop(idx_width)   = mesh{idx_width}(idx) + diff(mesh{idx_width}(idx:idx+1))/2;
+	i2_start = i1_start;
+	i2_stop  = i1_stop;
+	i2_start(idx_prop)   = sum(meshlines(2:3))/2;
+	i2_stop(idx_prop)    = i2_start(idx_prop);
+catch
+  error('Unable to place current probe on mesh; check the location of the MSL, and make sure that the mesh is large enough');
+end
 
 % create the probes
 port.U_filename{1} = [PortNamePrefix 'port_ut' num2str(portnr) 'A'];
 % weight = sign(stop(idx_height)-start(idx_height))
-weight = upsidedown;
+weight = 1;
 CSX = AddProbe( CSX, port.U_filename{1}, 0, 'weight', weight );
 CSX = AddBox( CSX, port.U_filename{1}, prio, v1_start, v1_stop );
 port.U_filename{2} = [PortNamePrefix 'port_ut' num2str(portnr) 'B'];
@@ -232,13 +233,17 @@ port.measplanepos = abs(v2_start(idx_prop) - start(idx_prop))*port.LengthScale;
 % port
 
 % create excitation (if enabled) and port resistance
-meshline = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), start(idx_prop) + feed_shift*direction, 'nearest' );
-ex_start(idx_prop)   = mesh{idx_prop}(meshline) ;
-ex_start(idx_width)  = nstart(idx_width);
-ex_start(idx_height) = nstart(idx_height);
-ex_stop(idx_prop)    = ex_start(idx_prop);
-ex_stop(idx_width)   = nstop(idx_width);
-ex_stop(idx_height)  = nstop(idx_height);
+try
+	meshline = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), start(idx_prop) + feed_shift*direction, 'nearest' );
+	ex_start(idx_prop)   = mesh{idx_prop}(meshline) ;
+	ex_start(idx_width)  = nstart(idx_width);
+	ex_start(idx_height) = nstart(idx_height);
+	ex_stop(idx_prop)    = ex_start(idx_prop);
+	ex_stop(idx_width)   = nstop(idx_width);
+	ex_stop(idx_height)  = nstop(idx_height);
+catch
+  error('Unable to place excitation on mesh; check the location of the MSL and the excitation (FeedShift), and make sure that the mesh is large enough');
+end
 
 port.excite = 0;
 if excite
@@ -247,7 +252,7 @@ if excite
     CSX = AddBox( CSX, [PortNamePrefix 'port_excite_' num2str(portnr)], prio, ex_start, ex_stop );
 end
 
-%% MSL resitance at start of MSL line
+%% MSL resistance at start of MSL line
 ex_start(idx_prop) = start(idx_prop);
 ex_stop(idx_prop) = ex_start(idx_prop);
 
@@ -260,6 +265,6 @@ elseif feed_R == 0
     %port "resistance" as metal
     CSX = AddBox( CSX, materialname, prio, ex_start, ex_stop );
 else
-    error('openEMS:AddMSLPort','MSL port with resitance <= 0 it not possible');
+    error('openEMS:AddMSLPort','MSL port with resistance <= 0 it not possible');
 end
 end

matlab/AddStripLinePort.m

@@ -4,7 +4,7 @@ function [CSX,port] = AddStripLinePort( CSX, prio, portnr, materialname, start,
 % CSX:          CSX-object created by InitCSX()
 % prio:         priority for excitation and probe boxes
 % portnr:       (integer) number of the port
-% materialname: property for the MSL (created by AddMetal())
+% materialname: property for the stripline (created by AddMetal())
 % start:        3D start rowvector for port definition
 % stop:         3D end rowvector for port definition
 % height:       height of the stripline (top and bottom)
@@ -17,7 +17,7 @@ function [CSX,port] = AddStripLinePort( CSX, prio, portnr, materialname, start,
 %               is false)
 % 'FeedShift'   shift to port from start by a given distance in drawing
 %               units. Default is 0. Only active if 'ExcitePort' is set!
-% 'Feed_R'      Specifiy a lumped port resistance. Default is no lumped
+% 'Feed_R'      Specify a lumped port resistance. Default is no lumped
 %               port resistance --> port has to end in an ABC.
 % 'MeasPlaneShift'  Shift the measurement plane from start t a given distance
 %               in drawing units. Default is the middle of start/stop.
@@ -68,7 +68,7 @@ evec0 = evec ./ sum(evec); % evec0 is a unit vector
 
 %set defaults
 feed_shift = 0;
-feed_R = inf; %(default is open, no resitance)
+feed_R = inf; %(default is open, no resistance)
 excite = false;
 measplanepos = nan;
 PortNamePrefix = '';
@@ -84,17 +84,17 @@ for n=1:2:numel(varargin)
         end
     elseif (strcmp(varargin{n},'Feed_R')==1);
         feed_R = varargin{n+1};
-        if (numel(feed_shift)>1)
+        if (numel(feed_R)>1)
             error 'Feed_R must be a scalar value'
         end
     elseif (strcmp(varargin{n},'MeasPlaneShift')==1);
         measplanepos = varargin{n+1};
-        if (numel(feed_shift)>1)
+        if (numel(measplanepos)>1)
             error 'MeasPlaneShift must be a scalar value'
         end
     elseif (strcmp(varargin{n},'ExcitePort')==1);
         if ischar(varargin{n+1})
-            warning('CSXCAD:AddMSLPort','depreceated: a string as excite option is no longer supported and will be removed in the future, please use true or false');
+            warning('CSXCAD:AddStripLinePort','depreceated: a string as excite option is no longer supported and will be removed in the future, please use true or false');
             if ~isempty(excite)
                 excite = true;
             else
@@ -117,10 +117,10 @@ end
 nstart = min( [start;stop] );
 nstop  = max( [start;stop] );
 
-% determine index (1, 2 or 3) of propagation (length of MSL)
+% determine index (1, 2 or 3) of propagation (length of stripline)
 idx_prop = dir + 1;
 
-% determine index (1, 2 or 3) of width of MSL
+% determine index (1, 2 or 3) of width of stripline
 dir = [0 0 0];
 dir(idx_prop) = 1;
 idx_width = abs(cross(dir,evec0)) * [1;2;3];
@@ -139,7 +139,7 @@ else
     direction = -1;
 end
 
-% create the metal/material for the MSL
+% create the metal/material for the stripline
 SL_start = start;
 SL_stop = stop;
 CSX = AddBox( CSX, materialname, prio, SL_start, SL_stop );
@@ -151,29 +151,33 @@ else
 end
 
 % calculate position of the voltage probes
-mesh{1} = sort(CSX.RectilinearGrid.XLines);
-mesh{2} = sort(CSX.RectilinearGrid.YLines);
-mesh{3} = sort(CSX.RectilinearGrid.ZLines);
-meshlines = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), measplanepos, 'nearest' );
-meshlines = mesh{idx_prop}(meshlines-1:meshlines+1); % get three lines (approx. at center)
-if direction == -1
+try
+	mesh{1} = sort(CSX.RectilinearGrid.XLines);
+	mesh{2} = sort(CSX.RectilinearGrid.YLines);
+	mesh{3} = sort(CSX.RectilinearGrid.ZLines);
+	meshlines = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), measplanepos, 'nearest' );
+	meshlines = mesh{idx_prop}(meshlines-1:meshlines+1); % get three lines (approx. at center)
+	if direction == -1
 		meshlines = fliplr(meshlines);
+	end
+	SL_w2 = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), (nstart(idx_width)+nstop(idx_width))/2, 'nearest' );
+	SL_w2 = mesh{idx_width}(SL_w2); % get e-line at center of stripline (SL_width/2)
+	v1_start(idx_prop)   = meshlines(1);
+	v1_start(idx_width)  = SL_w2;
+	v1_start(idx_height) = start(idx_height);
+	v1_stop  = v1_start;
+	v1_stop(idx_height)  = v1_start(idx_height);
+	v2_start = v1_start;
+	v2_stop  = v1_stop;
+	v2_start(idx_prop)   = meshlines(2);
+	v2_stop(idx_prop)    = meshlines(2);
+	v3_start = v2_start;
+	v3_stop  = v2_stop;
+	v3_start(idx_prop)   = meshlines(3);
+	v3_stop(idx_prop)    = meshlines(3);
+catch
+	error('Unable to place voltage probe on mesh; check the location of the stripline and the probe (MeasPlaneShift), and make sure that the mesh is large enough');
 end
-SL_w2 = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), (nstart(idx_width)+nstop(idx_width))/2, 'nearest' );
-SL_w2 = mesh{idx_width}(SL_w2); % get e-line at center of MSL (SL_width/2)
-v1_start(idx_prop)   = meshlines(1);
-v1_start(idx_width)  = SL_w2;
-v1_start(idx_height) = start(idx_height);
-v1_stop  = v1_start;
-v1_stop(idx_height)  = v1_start(idx_height);
-v2_start = v1_start;
-v2_stop  = v1_stop;
-v2_start(idx_prop)   = meshlines(2);
-v2_stop(idx_prop)    = meshlines(2);
-v3_start = v2_start;
-v3_stop  = v2_stop;
-v3_start(idx_prop)   = meshlines(3);
-v3_stop(idx_prop)    = meshlines(3);
 
 height_vector = [0 0 0];
 height_vector(idx_height) = height;
@@ -185,7 +189,7 @@ CSX = AddProbe( CSX, port.U_filename{1,1}, 0, 'weight', weight );
 CSX = AddBox( CSX, port.U_filename{1,1}, prio, v1_start, v1_stop+height_vector);
 
 port.U_filename{1,2} = [PortNamePrefix 'port_ut' num2str(portnr) 'A2'];
-CSX = AddProbe( CSX, port.U_filename{1,2}, 0, 'weight', -1*weight );
+CSX = AddProbe( CSX, port.U_filename{1,2}, 0, 'weight', weight );
 CSX = AddBox( CSX, port.U_filename{1,2}, prio, v1_start, v1_stop-height_vector);
 
 
@@ -194,7 +198,7 @@ CSX = AddProbe( CSX, port.U_filename{2,1}, 0, 'weight', weight );
 CSX = AddBox( CSX, port.U_filename{2,1}, prio, v2_start, v2_stop+height_vector );
 
 port.U_filename{2,2} = [PortNamePrefix 'port_ut' num2str(portnr) 'B2'];
-CSX = AddProbe( CSX, port.U_filename{2,2}, 0, 'weight', -1*weight );
+CSX = AddProbe( CSX, port.U_filename{2,2}, 0, 'weight', weight );
 CSX = AddBox( CSX, port.U_filename{2,2}, prio, v2_start, v2_stop-height_vector );
 
 
@@ -203,23 +207,27 @@ CSX = AddProbe( CSX, port.U_filename{3,1}, 0, 'weight', weight );
 CSX = AddBox( CSX, port.U_filename{3,1}, prio, v3_start, v3_stop+height_vector );
 
 port.U_filename{3,2} = [PortNamePrefix 'port_ut' num2str(portnr) 'C2'];
-CSX = AddProbe( CSX, port.U_filename{3,2}, 0, 'weight', -1*weight );
+CSX = AddProbe( CSX, port.U_filename{3,2}, 0, 'weight', weight );
 CSX = AddBox( CSX, port.U_filename{3,2}, prio, v3_start, v3_stop-height_vector );
 
 % calculate position of the current probes
-idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstart(idx_width), 'nearest' );
-i1_start(idx_width)  = mesh{idx_width}(idx) - diff(mesh{idx_width}(idx-1:idx))/2;
-idx = interp1( mesh{idx_height}, 1:numel(mesh{idx_height}), start(idx_height), 'nearest' );
-i1_start(idx_height) = mesh{idx_height}(idx-1) - diff(mesh{idx_height}(idx-2:idx-1))/2;
-i1_stop(idx_height)  = mesh{idx_height}(idx+1) + diff(mesh{idx_height}(idx+1:idx+2))/2;
-i1_start(idx_prop)   = sum(meshlines(1:2))/2;
-i1_stop(idx_prop)    = i1_start(idx_prop);
-idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstop(idx_width), 'nearest' );
-i1_stop(idx_width)   = mesh{idx_width}(idx) + diff(mesh{idx_width}(idx:idx+1))/2;
-i2_start = i1_start;
-i2_stop  = i1_stop;
-i2_start(idx_prop)   = sum(meshlines(2:3))/2;
-i2_stop(idx_prop)    = i2_start(idx_prop);
+try
+	idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstart(idx_width), 'nearest' );
+	i1_start(idx_width)  = mesh{idx_width}(idx) - diff(mesh{idx_width}(idx-1:idx))/2;
+	idx = interp1( mesh{idx_height}, 1:numel(mesh{idx_height}), start(idx_height), 'nearest' );
+	i1_start(idx_height) = mesh{idx_height}(idx-1) - diff(mesh{idx_height}(idx-2:idx-1))/2;
+	i1_stop(idx_height)  = mesh{idx_height}(idx+1) + diff(mesh{idx_height}(idx+1:idx+2))/2;
+	i1_start(idx_prop)   = sum(meshlines(1:2))/2;
+	i1_stop(idx_prop)    = i1_start(idx_prop);
+	idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstop(idx_width), 'nearest' );
+	i1_stop(idx_width)   = mesh{idx_width}(idx) + diff(mesh{idx_width}(idx:idx+1))/2;
+	i2_start = i1_start;
+	i2_stop  = i1_stop;
+	i2_start(idx_prop)   = sum(meshlines(2:3))/2;
+	i2_stop(idx_prop)    = i2_start(idx_prop);
+catch
+	error('Unable to place current probe on mesh; check the location of the stripline and the probe (MeasPlaneShift), and make sure that the mesh is large enough');
+end
 
 % create the curr-probes
 weight = direction;
@@ -246,13 +254,17 @@ port.measplanepos = abs(v2_start(idx_prop) - start(idx_prop))*port.LengthScale;
 % port
 
 % create excitation (if enabled) and port resistance
-meshline = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), start(idx_prop) + feed_shift*direction, 'nearest' );
-ex_start(idx_prop)   = mesh{idx_prop}(meshline) ;
-ex_start(idx_width)  = nstart(idx_width);
-ex_start(idx_height) = nstart(idx_height);
-ex_stop(idx_prop)    = ex_start(idx_prop);
-ex_stop(idx_width)   = nstop(idx_width);
-ex_stop(idx_height)  = nstop(idx_height);
+try
+	meshline = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), start(idx_prop) + feed_shift*direction, 'nearest' );
+	ex_start(idx_prop)   = mesh{idx_prop}(meshline) ;
+	ex_start(idx_width)  = nstart(idx_width);
+	ex_start(idx_height) = nstart(idx_height);
+	ex_stop(idx_prop)    = ex_start(idx_prop);
+	ex_stop(idx_width)   = nstop(idx_width);
+	ex_stop(idx_height)  = nstop(idx_height);
+catch
+	error('Unable to place excitation on mesh; check the location of the stripline and the probe (MeasPlaneShift), and make sure that the mesh is large enough');
+end
 
 port.excite = 0;
 if excite
@@ -263,7 +275,7 @@ if excite
     CSX = AddBox( CSX, [PortNamePrefix 'port_excite_2_' num2str(portnr)], prio, ex_start, ex_stop-height_vector );
 end
 
-%% MSL resitance at start of MSL line
+%% stripline resistance at start of stripline line
 ex_start(idx_prop) = start(idx_prop);
 ex_stop(idx_prop) = ex_start(idx_prop);
 
@@ -278,6 +290,6 @@ elseif feed_R == 0
     CSX = AddBox( CSX, materialname, prio, ex_start, ex_stop+height_vector );
     CSX = AddBox( CSX, materialname, prio, ex_start, ex_stop-height_vector );
 else
-    error('openEMS:AddMSLPort','MSL port with resitance <= 0 it not possible');
+    error('openEMS:AddStripLinePort','stripline port with resistance <= 0 it not possible');
 end
 end

matlab/AddWaveGuidePort.m

@@ -74,7 +74,7 @@ for n=1:2:numel(varargin_tmp)
     end
 end
 
-% matlab adressing
+% matlab addressing
 dir = dir + 1;
 dir_sign = sign(stop(dir) - start(dir));
 if (dir_sign==0)

matlab/CalcNF2FF.m

@@ -14,7 +14,7 @@ function nf2ff = CalcNF2FF(nf2ff, Sim_Path, freq, theta, phi, varargin)
 % freq:     array of frequencies to analyse
 % theta,phi: spherical coordinates to evaluate the far-field on (in radians)
 %
-% optional paramater:
+% optional parameter:
 % 'Center': nf2ff phase center, default is [0 0 0]
 %           !! Make sure the center is never outside of your nf2ff box!!
 %           Definition is the correct coordinate system necessary

matlab/ConvertHDF5_VTK.m

@@ -15,7 +15,7 @@ function ConvertHDF5_VTK(hdf_file, vtk_prefix, varargin)
 %   'FieldName': field name written to vtk, e.g. 'E-Field'
 %   'weight':   field weighting
 %
-%   for more optional aguments have a look at ReadHDF5Dump
+%   for more optional augments have a look at ReadHDF5Dump
 %
 % example:
 %   % read time-domian data from hdf5, perform dft and dump as vtk
@@ -56,7 +56,7 @@ end
 
 if (do_FD_dump)
     if (~isfield(field,'FD'))
-        warning('openEMS:ConvertHDF5_VTK','no FD data found skipping frequency domian vtk dump...');
+        warning('openEMS:ConvertHDF5_VTK','no FD data found skipping frequency domain vtk dump...');
     else
         %set weighting
         if (numel(weight)~=numel(field.FD.frequency))
@@ -88,7 +88,7 @@ end
 
 if (do_TD_dump)
     if (~isfield(field,'TD'))
-        warning('openEMS:ConvertHDF5_VTK','no TD data found skipping time domian vtk dump...');
+        warning('openEMS:ConvertHDF5_VTK','no TD data found skipping time domain vtk dump...');
     else
         disp('dumping time domain data...')
         acc = ['%0' int2str(ceil(log10(numel(field.TD.time)+1))) 'd'];

matlab/DelayFidelity.m

@@ -0,0 +1,94 @@
+function [delay, fidelity, nf2ff_out] = DelayFidelity(nf2ff, port, path, weight_theta, weight_phi, theta, phi, f_0, f_c, varargin)
+% [delay, fidelity] = DelayFidelity(nf2ff, port, path, weight_theta, weight_phi, theta, phi, f_lo, f_hi, varargin)
+% 
+% 
+% This function calculates the time delay from the source port to the phase center of the antenna and the fidelity. 
+% The fidelity is the similarity between the excitation pulse and the radiated pulse (normalized scalar product).
+% The resolution of the delay will be equal to or better than ((f_0 + f_c)*Oversampling)^-1 when using Gaussian excitation.
+% Oversampling is an input parameter to InitFDTD. The rows of delay and fidelity correspond to theta and the columns to phi.
+% 
+% input:
+%   nf2ff: return value of CreateNF2FFBox.
+%   port: return value of AddLumpedPort
+%   path: path of the simulation results.
+%   weight_theta: weight of the E_theta component 
+%   weight_phi: weight of the E_phi component
+%     -> with both (possibly complex) parameters any polarization can be examined
+%   theta: theta values to be simulated
+%   phi: phi values to be simulated
+%   f_0: center frequency of SetGaussExcite
+%   f_c: cutoff frequency of SetGaussExcite
+%
+% variable input:
+%   'Center': phase center of the antenna for CalcNF2FF
+%   'Radius': radius for CalcNF2FF
+%   'Mode': mode CalcNF2FF
+%   
+% example:
+%   theta = [-180:10:180] * pi / 180;
+%   phi = [0, 90] * pi / 180;
+%   % use circular right handed polarization
+%   [delay, fidelity] = DelayFidelity2(nf2ff, port, Sim_Path, -1i, 1, theta, phi, f_0, f_c, 'Mode', 1);
+%   figure
+%   polar(theta.', delay(:,1) * 3e11); % delay in mm
+%   figure
+%   polar(theta', (fidelity(:,1)-0.95)/0.05); % last 5 percent of fidelity
+% 
+% Author: Georg Michel 
+
+C0 = 299792458;
+center = [0, 0, 0];
+radius = 1;
+nf2ff_mode = 0;
+
+for n=1:2:numel(varargin)
+    if (strcmp(varargin{n},'Center')==1);
+        center = varargin{n+1};
+    elseif (strcmp(varargin{n},'Radius')==1);
+        radius = varargin{n+1};
+    elseif (strcmp(varargin{n},'Mode')==1);
+        nf2ff_mode = varargin{n+1};
+    end
+end
+
+
+port_ut = load(fullfile(path, port.U_filename));
+port_it = load(fullfile(path, port.I_filename));
+dt = port_ut(2,1) - port_ut(1,1);
+fftsize = 2^(nextpow2(size(port_ut)(1)) + 1);
+df = 1 / (dt * fftsize);
+uport = fft(port_ut(:, 2), fftsize)(1:fftsize/2+1);
+iport = fft(port_it(:, 2), fftsize)(1:fftsize/2+1);
+fport = df * (0:fftsize/2);
+f_ind = find(fport > (f_0 - f_c ) & fport < (f_0 + f_c));
+disp(["frequencies: ", num2str(numel(f_ind))]);
+exc_f = uport.' + iport.' * port.Feed_R; %excitation in freq domain
+exc_f(!f_ind) = 0;
+exc_f /= sqrt(exc_f * exc_f'); % normalization (transposing also conjugates)
+
+nf2ff = CalcNF2FF(nf2ff, path, fport(f_ind), theta, phi, ...
+        'Center', center, 'Radius', radius, 'Mode', nf2ff_mode);
+radfield = weight_theta * cell2mat(nf2ff.E_theta) + weight_phi * cell2mat(nf2ff.E_phi); % rows: theta(f1), columns: phi(f1), phi(f2), ...phi(fn)
+radfield = reshape(radfield, [length(nf2ff.theta), length(nf2ff.phi), length(nf2ff.freq)]);
+correction = reshape(exp(-2i*pi*nf2ff.r/C0*nf2ff.freq), 1,1,numel(nf2ff.freq)); %dimensions: theta, phi, frequencies
+radfield = radfield./correction; % correct for radius delay
+% normalize radfield
+radnorm = sqrt(dot(radfield, radfield, 3));
+radfield ./= radnorm;
+
+%initialize radiated field in fully populated frequency domain
+rad_f = zeros([numel(nf2ff.theta), numel(nf2ff.phi), numel(fport)]);
+rad_f(:, :, f_ind) = radfield; % assign selected frequencies
+exc_f = reshape(exc_f, [1,1,numel(exc_f)]); %make exc_f conformant with rad_f
+
+cr_f = rad_f .* conj(exc_f); % calculate cross correlation
+% calculate the cross correlation in time domain (analytic signal)
+cr = ifft(cr_f(:, :, 1:end-1), [], 3) * (numel(fport) -1); % twice the FFT normalization (sqrt^2) because product of two normalized functions 
+%search for the maximum of the envelope
+[fidelity, delay_ind] = max(abs(cr), [], 3);
+delay = (delay_ind - 1) * dt * 2; % double time step because of single-sided FFT
+nf2ff_out = nf2ff; %possibly needed for plotting the far field and other things
+disp(["DelayFidelity: delay resolution = ", num2str(dt*2e9), "ns"]);
+return;
+
+

matlab/Dump2VTK.m

@@ -1,7 +1,7 @@
 function Dump2VTK(filename, fields, mesh, fieldname, varargin)
 % Dump2VTK(filename, fields, mesh, fieldname, varargin)
 %
-%   Dump fields extraced from an hdf5 file to a vtk file format
+%   Dump fields extracted from an hdf5 file to a vtk file format
 %
 %   possible arguments:
 %       'NativeDump': 0 (default) / 1, dump in native coordinate system

matlab/DumpFF2VTK.m

@@ -5,16 +5,16 @@ function DumpFF2VTK(filename, farfield, thetaRange, phiRange, varargin)
 %
 % input:
 %   filename:      filename of VTK file, existing file will be overwritten
-%   farfield:      farfield in V/m
+%   farfield:      far field in V/m
 %   thetaRange:    theta range in deg
 %   phiRange:      phi range in deg
 %
 % variable input:
 %   'scale':       - linear scale of plot, doesn't affect gain values
-%   'logscale':    - if set, show farfield with logarithmic scale
+%   'logscale':    - if set, show far field with logarithmic scale
 %                  - set the dB value for point of origin
 %                  - values below will be clamped
-%   'maxgain':     - add max gain in dB to normalized farfield
+%   'maxgain':     - add max gain in dB to normalized far field
 %                  - only valid if logscale is set
 %                  - default is 0dB
 %

matlab/FindFreeSSH.m

@@ -44,7 +44,7 @@ end
 if ischar(host_list)
     fid=fopen(host_list);
     if (fid==-1)
-        error('FindFreeSSH: cant open host file');
+        error('FindFreeSSH: cannot open host file');
     end
     clear host_list;
     host_list = {};

matlab/GetField_Interpolation.m

@@ -3,9 +3,9 @@ function [field_i mesh_i] = GetField_Interpolation(field, mesh, lines, varargin)
 %
 %   Get an interpolated field, e.g. read by ReadHDF5Dump 
 % 
-%   homogen interpolation given by a 3x1 vector: e.g. [21,1,101]
+%   homogeneous interpolation given by a 3x1 vector: e.g. [21,1,101]
 % 
-%   abitrary interpolation on a given mesh:
+%   arbitrary interpolation on a given mesh:
 %               e.g.:   mesh_interp{1} = linspace(0,  1,101) * 1e-3;
 %                       mesh_interp{2} = linspace(0,0.5, 51) * 1e-3;
 %                       mesh_interp{3} = linspace(0,0.2, 21) * 1e-3;

matlab/GetField_TD2FD.m

@@ -2,7 +2,7 @@ function field = GetField_TD2FD(field, freq)
 % function field = GetField_TD2FD(field, freq)
 %
 % Transforms time-domain field data into the frequency domain
-% Autocorrects the half-timestep offset of the H-field
+% Auto-corrects the half-timestep offset of the H-field
 % 
 % example:
 %   freq = linspace(0,1e9,100); %target frequency vector (Hz)

matlab/InitFDTD.m

@@ -1,7 +1,7 @@
 function FDTD = InitFDTD(varargin)
 % function FDTD = InitFDTD(varargin)
 %
-% Inititalize the FDTD data-structure.
+% Initialize the FDTD data-structure.
 %
 % optional field arguments for usage with openEMS:
 %   NrTS:           max. number of timesteps to simulate (e.g. default=1e9)

matlab/PlotHDF5FieldData.m

@@ -25,7 +25,7 @@ mesh = ReadHDF5Mesh(file);
 fields = ReadHDF5FieldData(file);
 
 if (mesh.type==0)
-    % cartesian mesh
+    % Cartesian mesh
     [X Y Z] = meshgrid(mesh.lines{1},mesh.lines{2},mesh.lines{3});
     for n=1:numel(fields.TD.values)
         % since Matlab 7.1SP3 the field needs to be reordered

matlab/ReadHDF5Mesh.m

@@ -4,7 +4,7 @@ function hdf_mesh = ReadHDF5Mesh(file)
 %   Get the raw mesh data stored in the hdf5 dump file created by openEMS
 %
 % returns:
-% hdf_mesh.type     (0-> cartesian, 1-> cylindrical mesh type)
+% hdf_mesh.type     (0-> Cartesian, 1-> cylindrical mesh type)
 % hdf_mesh.names    (e.g. 'Mesh/y')
 % hdf_mesh.lines    (e.g. [0,1,2,3,4])
 %
@@ -56,7 +56,7 @@ hdf_mesh.type=0;
 function hdf_mesh = ReadHDF5Mesh_octave(file)
 hdf = load( '-hdf5', file );
 hdf_mesh.names = fieldnames(hdf.Mesh);
-hdf_mesh.type = 0; % cartesian mesh
+hdf_mesh.type = 0; % Cartesian mesh
 for n=1:numel(hdf_mesh.names)
     hdf_mesh.lines{n} = hdf.Mesh.(hdf_mesh.names{n});
     hdf_mesh.names{n} = ['/Mesh/' hdf_mesh.names{n}];

matlab/RunOpenEMS.m

@@ -4,7 +4,7 @@ function RunOpenEMS(Sim_Path, Sim_File, opts, Settings)
 % Run an openEMS simulation.
 %
 % arguments:
-% Sim_Path: specifiy the simulation folder (folder must exist!)
+% Sim_Path: specify the simulation folder (folder must exist!)
 % Sim_File: xml-filename to simulate, created by WriteOpenEMS
 %
 % optional arguments
@@ -157,7 +157,7 @@ if (enable_ssh)
     disp( 'Remote simulation done... copying back results and cleaning up...' );
 
     %copy back all results
-    [stat, res] = system([scp_command ' -r ' scp_options ' ' Settings.SSH.host ':' ssh_work_path '/* ' pwd '/']);
+    [stat, res] = system([scp_command ' -r ' scp_options ' ' Settings.SSH.host ':' ssh_work_path '/* ./']);
     if (stat~=0);
         disp(res);
         error('openEMS:RunOpenEMS','scp failed!');
@@ -170,11 +170,11 @@ if (enable_ssh)
         warning('openEMS:RunOpenEMS','remote cleanup failed!');
     end
 else
-    args = [Sim_File ' ' opts];
+    args = ['"' Sim_File '" ' opts]
     if isfield(Settings,'LogFile') && isfield(Settings,'Silent')
-        invoke_openEMS(args,Settings.LogFile,Settings.Silent);
+        invoke_openEMS(args,['"' Settings.LogFile '"'],Settings.Silent);
     elseif isfield(Settings,'LogFile')
-        invoke_openEMS(args,Settings.LogFile);
+        invoke_openEMS(args,['"' Settings.LogFile '"']);
     elseif isfield(Settings,'Silent')
         invoke_openEMS(args,[],Settings.Silent);
     else

matlab/RunOpenEMS_MPI.m

@@ -87,7 +87,7 @@ end
 
 if isfield(Settings.MPI,'Hosts')
     disp(['Running remote openEMS_MPI in working dir: ' work_path]);
-    [status]  = system(['mpiexec -host ' HostList ' -n ' int2str(NrProc) ' -wdir ' work_path ' ' Settings.MPI.Binary ' ' Sim_File ' ' opts ' ' append_unix]);
+    [status]  = system(['mpiexec ' Settings.MPI.GlobalArgs ' -host ' HostList ' -n ' int2str(NrProc) ' -wdir ' work_path ' ' Settings.MPI.Binary ' ' Sim_File ' ' opts ' ' append_unix]);
 else
     disp('Running local openEMS_MPI');
     [status]  = system(['mpiexec ' Settings.MPI.GlobalArgs ' -n ' int2str(NrProc) ' ' Settings.MPI.Binary ' ' Sim_File ' ' opts ' ' append_unix]);
@@ -100,7 +100,7 @@ end
 if isfield(Settings.MPI,'Hosts')
     disp( 'Remote simulation done... copying back results and cleaning up...' );
     
-    if (strncmp(work_path,'/tmp/',5)~=1) % savety precaution...
+    if (strncmp(work_path,'/tmp/',5)~=1) % safety precaution...
         error('openEMS:RunOpenEMS','working path invalid for deletion');
     end
     

matlab/RunOpenEMS_Parallel.m

@@ -7,9 +7,9 @@ function [stdout, stderr] = RunOpenEMS_Parallel(Sim_Paths, Sim_Files, opts, Sett
 % This function relies on InitQueue etc.
 % 
 % input:
-%       Sim_Paths:  cell array of pathes to simulate by RunOpenEMS
+%       Sim_Paths:  cell array of paths to simulate by RunOpenEMS
 %       Sim_Files:  filename or cell array of filenames to simulate
-%       opts:       openEMS options. sa RunOpenEMS
+%       opts:       openEMS options. see also RunOpenEMS
 %       Settings:   use the settings to define multiple host for simulation
 %                   e.g.: Settings.SSH.bin ='<path_to_openEMS>/openEMS.sh';
 %                         Settings.SSH.host_list = {'list','of','hosts'};

matlab/SetBoundaryCond.m

@@ -11,7 +11,7 @@ function FDTD = SetBoundaryCond(FDTD, BC, varargin)
 % 
 % example:
 % BC = [ 1     1     0     0     2     3     ]  %using numbers or
-% BC = {'PMC' 'PMC' 'PEC' 'PEC' 'MUR' 'PML_8'}  %usign equivalent strings
+% BC = {'PMC' 'PMC' 'PEC' 'PEC' 'MUR' 'PML_8'}  %using equivalent strings
 %
 % mur-abc definitions
 % define a phase-velocity to be used by the mur-abc

matlab/SetCustomExcite.m

@@ -2,7 +2,7 @@ function FDTD = SetCustomExcite(FDTD,f0,funcStr)
 % function FDTD = SetCustomExcite(FDTD,f0,funcStr)
 %
 % f0 : nyquist rate
-% funcStr : string desribing the excitation function e(t)
+% funcStr : string describing the excitation function e(t)
 %
 % see also SetSinusExcite SetGaussExcite
 %

matlab/Tutorials/Bent_Patch_Antenna.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / bent patch antenna
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_Bent_Patch_Antenna
 %
 % Tested with
@@ -76,7 +76,7 @@ start = [patch.radius+substrate.thickness -substr_ang_width/2 -substrate.length/
 stop  = [patch.radius+substrate.thickness +substr_ang_width/2  substrate.length/2];
 CSX = AddBox( CSX, 'Jt_patch', 0, start, stop );
 
-%% create ground (not really necessary, only for esthetic reasons)
+%% create ground (not really necessary, only for aesthetic reasons)
 CSX = AddMetal( CSX, 'gnd' ); % create a perfect electric conductor (PEC)
 start = [patch.radius -substr_ang_width/2 -substrate.length/2];
 stop  = [patch.radius +substr_ang_width/2 +substrate.length/2];
@@ -158,7 +158,7 @@ ylabel( 'reflection coefficient |S_{11}|' );
 
 drawnow
 
-%find resonance frequncy from s11
+%find resonance frequency from s11
 f_res_ind = find(s11==min(s11));
 f_res = freq(f_res_ind);
 

matlab/Tutorials/CRLH_Extraction.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / CRLH_Extraction
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_CRLH_Extraction
 %
 % Tested with

matlab/Tutorials/CRLH_LeakyWaveAnt.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / CRLH_LeakyWaveAnt
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_CRLH_Leaky_Wave_Antenna
 %
 % Tested with

matlab/Tutorials/Circ_Waveguide.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / Circ_Waveguide
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_Circular_Waveguide
 %
 % Tested with

matlab/Tutorials/Conical_Horn_Antenna.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / conical horn antenna
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_Conical_Horn_Antenna
 %
 % Tested with

matlab/Tutorials/CylindricalWave_CC.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / CylindricalWave_CC
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_2D_Cylindrical_Wave
 %
 % Tested with
@@ -28,7 +28,7 @@ exite_offset = 1300;
 excite_angle = 45;
 
 %% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-FDTD = InitFDTD(100000,1e-4,'CoordSystem',1,'MultiGrid',split);
+FDTD = InitFDTD('NrTS', 100000, 'EndCriteria', 1e-4, 'CoordSystem', 1, 'MultiGrid', split);
 FDTD = SetGaussExcite(FDTD,f0,f0/2);
 BC = [0 3 0 0 0 0];             % pml in positive r-direction
 FDTD = SetBoundaryCond(FDTD,BC);

matlab/Tutorials/Dipole_SAR.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / Dipole SAR + Power budget
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_Dipole_SAR
 %
 % Tested with

matlab/Tutorials/Helical_Antenna.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / helical antenna
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_Helical_Antenna
 %
 % Tested with

matlab/Tutorials/Horn_Antenna.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / horn antenna
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_Horn_Antenna
 %
 % Tested with

matlab/Tutorials/MRI_LP_Birdcage.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / 3T MRI Low Pass Birdcage coil
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_MRI_LP_Birdcage
 %
 % Estimated time to run:    ~7h @ ~65MC/s
@@ -69,7 +69,7 @@ body_model_transform = {'Rotate_X',pi,'Rotate_Z',pi, ...
                         'Translate',[0,5,-720]};
 
 %% some internal parameter
-physical_constants % load important physical constans
+physical_constants % load important physical constants
 end_crit = 1e-5;    %abort simulation at -50dB energy drop
 unit = 1e-3;        %drawing unit used
 
@@ -164,7 +164,7 @@ for n=1:BC.N_rungs
     stop  = [BC.rad a0+da_Segs/2+da_Strip/2 -0.5*BC.portlength];
     CSX = AddBox(CSX,'metal',1, start, stop);
     
-    % some additonal mesh lines
+    % some additional mesh lines
     mesh.a = [mesh.a a0+da_Segs/2];
     
     a0 = a0 + da_Segs;

matlab/Tutorials/MRI_Loop_Coil.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / 7T MRI Loop Coil
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_MRI_Loop_Coil
 %
 % Tested with
@@ -145,7 +145,7 @@ CSX = AddBox(CSX, 'body_model', 0, body_box.start, body_box.stop);
 % create loop mesh
 mesh = DetectEdges(CSX);
 
-% add a dense homegeneous mesh inside the human body model
+% add a dense homogeneous mesh inside the human body model
 mesh.x = [mesh.x mesh_box.start(1) mesh_box.stop(1)];
 mesh.y = [mesh.y mesh_box.start(2) mesh_box.stop(2)];
 mesh.z = [mesh.z mesh_box.start(3) mesh_box.stop(3)];
@@ -178,7 +178,7 @@ CSX = AddBox(CSX,'SAR_xz',0, body_box.start.*[1 0 1], body_box.stop.*[1 0 1]);
 %% boundary conditions
 mesh = AddPML(mesh, 10);
 
-%% finaly define the FDTD mesh grid
+%% finally define the FDTD mesh grid
 disp(['number of cells: ' num2str(1e-6*numel(mesh.x)*numel(mesh.y)*numel(mesh.z)) ' Mcells'])
 CSX = DefineRectGrid( CSX, unit, mesh );
 

matlab/Tutorials/MSL_NotchFilter.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / MSL_NotchFilter
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_Microstrip_Notch_Filter
 %
 % Tested with

matlab/Tutorials/Parallel_Plate_Waveguide.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / Parallel_Plate_Waveguide
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_Parallel_Plate_Waveguide
 %
 % Tested with

matlab/Tutorials/Patch_Antenna_Array.m

@@ -6,10 +6,10 @@ function [port nf2ff] = Patch_Antenna_Array(Sim_Path, postproc_only, show_struct
 %
 % Sim_Path: Simulation path
 % postproc_only: set to post process only 0/1
-% show_structure: show the strucuture in AppCSXCAD 0/1
+% show_structure: show the structure in AppCSXCAD 0/1
 % xpos: the x-position for each antenna is defined
 % caps: the port capacity (will override active port)
-% resist: port resitance
+% resist: port resistance
 % active: switch port active
 %
 % References:

matlab/Tutorials/Patch_Antenna_Phased_Array.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / Patch Antenna Phased Array
 %
-% Describtion at:
+% Description at:
 %
 % Tested with
 %  - Matlab 2011a
@@ -122,7 +122,7 @@ disp(['I2/I1: Matlab: ' num2str(I_out(2)/I_out(1))])
 disp(['I3/I1: Matlab: ' num2str(I_out(3)/I_out(1))])
 
 
-%% do a referenc simulation for the given C2/C3 values
+%% do a reference simulation for the given C2/C3 values
 if (do_reference_simulation)
     active = [1 0 0];
     caps = [0 C2 C3];
@@ -132,7 +132,7 @@ if (do_reference_simulation)
     port_ref = calcPort( port_ref, Sim_Path, f, 'RefImpedance', 50);
     nf2ff_ref = CalcNF2FF(nf2ff_ref, Sim_Path, f0, [-180:2:180]*pi/180, 0);
 
-    % extract currents from referenc simulation
+    % extract currents from reference simulation
     for p=1:3
         I_ref(p,1) = interp1(f, port_ref{p}.if.tot,f0);
     end
@@ -141,7 +141,7 @@ if (do_reference_simulation)
     disp(['I3/I1: openEMS: ' num2str(I_ref(3)/I_ref(1))])
 end
 
-%% calculate and apply weighting cooefficients [3]
+%% calculate and apply weighting coefficients [3]
 % calculate
 coeff = I\I_out;
 

matlab/Tutorials/RCS_Sphere.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / radar cross section of a metal sphere
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_RCS_Sphere
 %
 % Tested with

matlab/Tutorials/RadarUWBTutorial.m

@@ -0,0 +1,236 @@
+% Tutorial on time delay and signal integrity for radar 
+% and UWB applications
+% 
+% Tested with
+%  - Octave 4.0
+%  - openEMS v0.0.35
+% 
+% Author: Georg Michel, 2016
+
+ clear;
+ close all;
+
+physical_constants;
+
+% --- start of configuration section ---
+
+% In radar and ultrawideband applications it is important to know the
+% delay and fidelity of RF pulses. The delay is the retardation of the
+% signal from the source to the phase center of the antenna. It is
+% composed out of linear delay, dispersion and minimum-phase
+% delay. Dispersion due to waveguides or frequency-dependent
+% permittivity and minimum-phase delay due to resonances will degrade
+% the fidelity which is the normalized similarity between excitation and
+% radiated signal. In this tutorial you can examine the performance of a
+% simple ultrawideband (UWB) monopole. The delay and fidelity of this
+% antenna are calculated and plotted. You can compare these properties
+% in different channels.
+% 
+% The Gaussian excitation is set to the same 3dB bandwidth as the
+% channels of the IEEE 802.15.4 UWB PHY. One exception is channel4twice
+% which has the double bandwidth of channel 4. It can be seen that the
+% delay is larger and the fidelity is smaller in the vicinity of the
+% (undesired) resonances of the antenna. Note that for a real UWB system
+% the total delay and fidelity result from both the transmitting and
+% receiving antenna or twice the delay and the square of the fidelity
+% for monostatic radar. 
+% 
+% The resolution of the delay will depend on the 'Oversampling'
+% parameter to InitFDTD.  See the description of DelayFidelity
+% 
+% In the configuration section below you can uncomment the respective
+% parameter settings. As an exercise, you can examine how the permittivity 
+% of the substrate influences gain, delay and fidelity.  
+
+
+%suffix = "channel1";
+%f_0 = 3.5e9; % center frequency of the channel
+%f_c = 0.25e9 / 0.3925; % 3dB bandwidth is 0.3925 times 20dB bandwidth for Gaussian excitation
+
+%suffix = "channel2";
+%f_0 = 4.0e9; % center frequency of the channel
+%f_c = 0.25e9 / 0.3925; 
+
+%suffix = "channel3";
+%f_0 = 4.5e9; % center frequency of the channel
+%f_c = 0.25e9 / 0.3925; 
+
+suffix = "channel4";
+f_0 = 4.0e9; % center frequency of the channel
+f_c = 0.5e9 / 0.3925; 
+
+%suffix = "channel5";
+%f_0 = 6.5e9; % center frequency of the channel
+%f_c = 0.25e9 / 0.3925; 
+
+%suffix = "channel7";
+%f_0 = 6.5e9; % center frequency of the channel
+%f_c = 0.5e9 / 0.3925; 
+
+%suffix = "channel4twice"; % this is just to demonstrate the degradation of the fidelity with increasing bandwidth
+%f_0 = 4.0e9; % center frequency of the channel
+%f_c = 1e9 / 0.3925; 
+
+tilt = 45 * pi / 180; % polarization tilt angle against co-polarization (90DEG is cross polarized) 
+
+% --- end of configuration section ---
+
+% path and filename setup
+Sim_Path = 'tmp';
+Sim_CSX = 'uwb.xml';
+
+% properties of the substrate
+substrate.epsR = 4; % FR4
+substrate.height = 0.707;
+substrate.cells = 3; % thickness in cells
+
+% size of the monopole and the gap to the ground plane
+gap = 0.62; % 0.5
+patchsize = 14;
+
+% we will use millimeters
+unit = 1e-3;
+
+% set the resolution for the finer structures, e.g. the antenna gap
+fineResolution = C0 / (f_0 + f_c) / sqrt(substrate.epsR) / unit / 40;
+% set the resolution for the coarser structures, e.g. the surrounding air
+coarseResolution = C0/(f_0 + f_c) / unit / 20;
+
+
+% initialize the CSX structure
+CSX = InitCSX();
+
+% add the properties which are used to model the antenna
+CSX = AddMetal(CSX, 'Ground' );
+CSX = AddMetal(CSX, 'Patch');
+CSX = AddMetal(CSX, 'Line');
+CSX = AddMaterial(CSX, 'Substrate' );
+CSX = SetMaterialProperty(CSX, 'Substrate', 'Epsilon', substrate.epsR);
+
+% define the supstrate and sheet-like primitives for the properties
+CSX = AddBox(CSX, 'Substrate', 1, [-16, -16, -substrate.height], [16, 18, 0]);
+CSX = AddBox(CSX, 'Ground', 2, [-16, -16, -substrate.height], [16, 0, -substrate.height]);
+CSX = AddBox(CSX, 'Line', 2, [-1.15, -16, 0], [1.15, gap, 0]);
+CSX = AddBox(CSX, 'Patch', 2, [-patchsize/2, gap, 0], [patchsize/2, gap + patchsize, 0]);
+
+% setup a mesh
+mesh.x = [];
+mesh.y = [];
+
+% two mesh lines for the metal coatings of the substrate
+mesh.z = linspace(-substrate.height, 0, substrate.cells +1);
+
+% find optimal mesh lines for the patch and ground, not yes the microstrip line
+mesh = DetectEdges(CSX, mesh, 'SetProperty',{'Patch', 'Ground'}, '2D_Metal_Edge_Res', fineResolution/2);
+
+%replace gap mesh lines which are too close by a single mesh line
+tooclose = find (diff(mesh.y) < fineResolution/4);
+if ~isempty(tooclose)
+  mesh.y(tooclose) = (mesh.y(tooclose) + mesh.y(tooclose+1))/2;
+  mesh.y(tooclose + 1) = [];
+endif
+
+% store the microstrip  edges in a temporary variable
+meshline = DetectEdges(CSX, [], 'SetProperty', 'Line', '2D_Metal_Edge_Res', fineResolution/2);
+% as well as the edges of the substrate (without 1/3 - 2/3 rule)
+meshsubstrate = DetectEdges(CSX, [], 'SetProperty', 'Substrate');
+% add only the x mesh lines of the microstrip 
+mesh.x = [mesh.x meshline.x];
+% and only the top of the substrate, the other edges are covered by the ground plane
+mesh.y = [mesh.y, meshsubstrate.y(end)]; % top of substrate
+
+% for now we have only the edges, now calculate mesh lines in between
+mesh = SmoothMesh(mesh, fineResolution);
+
+% add the outer boundary
+mesh.x = [mesh.x -60, 60];
+mesh.y = [mesh.y, -60, 65];
+mesh.z = [mesh.z, -46, 45];
+
+% add coarse mesh lines for the free space
+mesh = SmoothMesh(mesh, coarseResolution);
+
+% define the grid
+CSX = DefineRectGrid( CSX, unit, mesh);
+% and the feeding port
+[CSX, port] = AddLumpedPort( CSX, 999, 1, 50, [-1.15, meshline.y(2), -substrate.height], [1.15, meshline.y(2), 0], [0 0 1], true);
+
+%setup a NF2FF box for the calculation of the far field
+start = [mesh.x(10)     mesh.y(10)     mesh.z(10)];
+stop  = [mesh.x(end-9) mesh.y(end-9) mesh.z(end-9)];
+[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop);
+
+% initialize the FDTD structure with excitation and open boundary conditions
+FDTD = InitFDTD( 'NrTs', 30000, 'EndCriteria', 1e-5, 'OverSampling', 20);
+FDTD = SetGaussExcite(FDTD, f_0, f_c );
+BC   = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'};
+FDTD = SetBoundaryCond(FDTD, BC );
+
+
+% remove old data, show structure, calculate new data
+[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
+[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder
+
+% write the data to the working directory
+WriteOpenEMS([Sim_Path '/' Sim_CSX], FDTD, CSX);
+% show the geometry for checking
+CSXGeomPlot([Sim_Path '/' Sim_CSX]);
+% run the simulation
+RunOpenEMS( Sim_Path, Sim_CSX);
+
+% plot amplitude and phase of the reflection coefficient
+freq  = linspace(f_0-f_c, f_0+f_c, 200);
+port = calcPort(port, Sim_Path, freq);
+s11 = port.uf.ref ./ port.uf.inc;
+s11phase = unwrap(arg(s11));
+figure %("visible", "off"); % use this to plot only into files at the end of this script
+ax = plotyy( freq/1e6, 20*log10(abs(s11)), freq/1e6, s11phase);
+grid on
+title( ['reflection coefficient ', suffix, ' S_{11}']);
+xlabel( 'frequency f / MHz' );
+ylabel( ax(1), 'reflection coefficient |S_{11}|' );
+ylabel(ax(2), 'S_{11} phase (rad)');
+
+% define an azimuthal trace around the monopole
+phi = [0] * pi / 180;
+theta = [-180:10:180] * pi / 180;
+
+% calculate the delay, the fidelity and the farfield
+[delay, fidelity, nf2ff] = DelayFidelity(nf2ff, port, Sim_Path, sin(tilt), cos(tilt), theta, phi, f_0, f_c, 'Mode', 1);
+
+%plot the gain at (close to) f_0
+f_0_nearest_ind = nthargout(2, @min, abs(nf2ff.freq -f_0));
+%turn directivity into gain
+nf2ff.Dmax(f_0_nearest_ind) *= nf2ff.Prad(f_0_nearest_ind) / calcPort(port, Sim_Path, nf2ff.freq(f_0_nearest_ind)).P_inc; 
+figure %("visible", "off");
+polarFF(nf2ff, 'xaxis', 'theta', 'freq_index', f_0_nearest_ind, 'logscale', [-4, 4]);
+title(["gain ", suffix, " / dBi"]);
+
+
+% We trick polarFF into plotting the delay in mm because
+% the axes of the vanilla polar plot can not be scaled.
+plotvar = delay * C0 * 1000; 
+maxplot = 80;
+minplot = 30;
+nf2ff.Dmax(1) = 10^(max(plotvar)/10);
+nf2ff.E_norm{1} = 10.^(plotvar/20);
+figure %("visible", "off");
+polarFF(nf2ff, 'xaxis', 'theta', 'logscale', [minplot, maxplot]);
+title(["delay ", suffix, " / mm"]);
+
+% The same for the fidelity in percent.
+plotvar = fidelity * 100;
+maxplot = 100;
+minplot = 98;
+nf2ff.Dmax(1) = 10^(max(plotvar)/10);
+nf2ff.E_norm{1} = 10.^(plotvar/20);
+figure %("visible", "off");
+polarFF(nf2ff, 'xaxis', 'theta', 'logscale', [minplot, maxplot]);
+title(["fidelity ", suffix, " / %"]);
+
+% save the plots in order to compare them after simulating the different channels
+print(1, ["s11_", suffix, ".png"]);
+print(2, ["farfield_", suffix, ".png"]);
+print(3, ["delay_mm_", suffix, ".png"]);
+print(4, ["fidelity_", suffix, ".png"]);
+return;

matlab/Tutorials/Rect_Waveguide.m

@@ -1,7 +1,7 @@
 %
 % Tutorials / Rect_Waveguide
 %
-% Describtion at:
+% Description at:
 % http://openems.de/index.php/Tutorial:_Rectangular_Waveguide
 %
 % Tested with
@@ -21,7 +21,7 @@ unit = 1e-6; %drawing unit in um
 % waveguide dimensions
 % WR42
 a = 10700;   %waveguide width
-b = 4300;    %waveguide heigth
+b = 4300;    %waveguide height
 length = 50000;
 
 % frequency range of interest

matlab/Tutorials/Simple_Patch_Antenna.m

@@ -1,20 +1,20 @@
+%% Simple Patch Antenna Tutorial
 %
-% Tutorials / simple patch antenna
-%
-% Describtion at:
-% http://openems.de/index.php/Tutorial:_Simple_Patch_Antenna
+% Description at:
+% <http://openems.de/index.php/Tutorial:_Simple_Patch_Antenna>
 %
 % Tested with
 %  - Matlab 2013a / Octave 4.0
-%  - openEMS v0.0.33
+%  - openEMS v0.0.35
 %
-% (C) 2010-2015 Thorsten Liebig <thorsten.liebig@uni-due.de>
+% (C) 2010-2017 Thorsten Liebig <thorsten.liebig@uni-due.de>
+%%
 
 close all
 clear
 clc
 
-%% setup the simulation
+%% Setup the Simulation
 physical_constants;
 unit = 1e-3; % all length in mm
 
@@ -38,7 +38,7 @@ feed.R = 50;     %feed resistance
 % size of the simulation box
 SimBox = [200 200 150];
 
-%% setup FDTD parameter & excitation function
+%% Setup FDTD Parameter & Excitation Function
 f0 = 2e9; % center frequency
 fc = 1e9; % 20 dB corner frequency
 FDTD = InitFDTD( 'NrTs', 30000 );
@@ -46,7 +46,7 @@ FDTD = SetGaussExcite( FDTD, f0, fc );
 BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % boundary conditions
 FDTD = SetBoundaryCond( FDTD, BC );
 
-%% setup CSXCAD geometry & mesh
+%% Setup CSXCAD Geometry & Mesh
 CSX = InitCSX();
 
 %initialize the mesh with the "air-box" dimensions
@@ -54,13 +54,13 @@ mesh.x = [-SimBox(1)/2 SimBox(1)/2];
 mesh.y = [-SimBox(2)/2 SimBox(2)/2];
 mesh.z = [-SimBox(3)/3 SimBox(3)*2/3];
 
-%% create patch
+% Create Patch
 CSX = AddMetal( CSX, 'patch' ); % create a perfect electric conductor (PEC)
 start = [-patch.width/2 -patch.length/2 substrate.thickness];
 stop  = [ patch.width/2  patch.length/2 substrate.thickness];
 CSX = AddBox(CSX,'patch',10,start,stop); % add a box-primitive to the metal property 'patch'
 
-%% create substrate
+% Create Substrate
 CSX = AddMaterial( CSX, 'substrate' );
 CSX = SetMaterialProperty( CSX, 'substrate', 'Epsilon', substrate.epsR, 'Kappa', substrate.kappa );
 start = [-substrate.width/2 -substrate.length/2 0];
@@ -70,18 +70,18 @@ CSX = AddBox( CSX, 'substrate', 0, start, stop );
 % add extra cells to discretize the substrate thickness
 mesh.z = [linspace(0,substrate.thickness,substrate.cells+1) mesh.z];
 
-%% create ground (same size as substrate)
+% Create Ground same size as substrate
 CSX = AddMetal( CSX, 'gnd' ); % create a perfect electric conductor (PEC)
 start(3)=0;
 stop(3) =0;
 CSX = AddBox(CSX,'gnd',10,start,stop);
 
-%% apply the excitation & resist as a current source
+% Apply the Excitation & Resist as a Current Source
 start = [feed.pos 0 0];
 stop  = [feed.pos 0 substrate.thickness];
 [CSX port] = AddLumpedPort(CSX, 5 ,1 ,feed.R, start, stop, [0 0 1], true);
 
-%% finalize the mesh
+% Finalize the Mesh
 % detect all edges except of the patch
 mesh = DetectEdges(CSX, mesh,'ExcludeProperty','patch');
 % detect and set a special 2D metal edge mesh for the patch
@@ -90,35 +90,41 @@ mesh = DetectEdges(CSX, mesh,'SetProperty','patch','2D_Metal_Edge_Res', c0/(f0+f
 mesh = SmoothMesh(mesh, c0/(f0+fc)/unit/20);
 CSX = DefineRectGrid(CSX, unit, mesh);
 
-%% add a nf2ff calc box; size is 3 cells away from MUR boundary condition
+CSX = AddDump(CSX,'Hf', 'DumpType', 11, 'Frequency',[2.4e9]);
+CSX = AddBox(CSX,'Hf',10,[-substrate.width -substrate.length -10*substrate.thickness],[substrate.width +substrate.length 10*substrate.thickness]); %assign box
+
+% add a nf2ff calc box; size is 3 cells away from MUR boundary condition
 start = [mesh.x(4)     mesh.y(4)     mesh.z(4)];
 stop  = [mesh.x(end-3) mesh.y(end-3) mesh.z(end-3)];
 [CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop);
 
-%% prepare simulation folder
+%% Prepare and Run Simulation
 Sim_Path = 'tmp_Patch_Ant';
 Sim_CSX = 'patch_ant.xml';
 
+% create an empty working directory
 [status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
 [status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder
 
-%% write openEMS compatible xml-file
+% write openEMS compatible xml-file
 WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
 
-%% show the structure
+% show the structure
 CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
 
-%% run openEMS
+% run openEMS
 RunOpenEMS( Sim_Path, Sim_CSX);
 
-%% postprocessing & do the plots
+%% Postprocessing & Plots
 freq = linspace( max([1e9,f0-fc]), f0+fc, 501 );
 port = calcPort(port, Sim_Path, freq);
 
-Zin = port.uf.tot ./ port.if.tot;
-s11 = port.uf.ref ./ port.uf.inc;
+%% Smith chart port reflection
+plotRefl(port, 'threshold', -10)
+title( 'reflection coefficient' );
 
 % plot feed point impedance
+Zin = port.uf.tot ./ port.if.tot;
 figure
 plot( freq/1e6, real(Zin), 'k-', 'Linewidth', 2 );
 hold on
@@ -130,6 +136,7 @@ ylabel( 'impedance Z_{in} / Ohm' );
 legend( 'real', 'imag' );
 
 % plot reflection coefficient S11
+s11 = port.uf.ref ./ port.uf.inc;
 figure
 plot( freq/1e6, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 );
 grid on
@@ -139,8 +146,8 @@ ylabel( 'reflection coefficient |S_{11}|' );
 
 drawnow
 
-%% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%find resonance frequncy from s11
+%% NFFF Plots
+%find resonance frequency from s11
 f_res_ind = find(s11==min(s11));
 f_res = freq(f_res_ind);
 
@@ -166,7 +173,7 @@ plotFFdB(nf2ff,'xaxis','theta','param',[1 2])
 
 drawnow
 
-%%
+% Show 3D pattern
 disp( 'calculating 3D far field pattern and dumping to vtk (use Paraview to visualize)...' );
 thetaRange = (0:2:180);
 phiRange = (0:2:360) - 180;

matlab/Tutorials/StripLine2MSL.m

@@ -0,0 +1,133 @@
+%
+% Stripline to Microstrip Line Transition
+%
+% Description at:
+% <http://openems.de/index.php/Tutorial:_Stripline_to_MSL_Transition>
+%
+% Tested with
+%  - Octave 4.0
+%  - openEMS v0.0.35
+%
+% (C) 2017 Thorsten Liebig <thorsten.liebig@gmx.de>
+
+close all
+clear
+clc
+
+%% Setup the Simulation
+physical_constants;
+unit = 1e-6; % specify everything in um
+
+line_length     = 15000; % line length of strip line and microstrip line
+substrate_width = 6000;
+air_spacer      = 4000;  % air spacer above the substrate
+
+msl_width               = 500;
+msl_substrate_thickness = 254;
+
+strip_width               = 500;
+strip_substrate_thickness = 512;
+
+connect_via_rad =  500/2;
+connect_via_gap = 1250/2;
+
+substrate_epr    = 3.66;
+substrate_kappa  = 1e-3 * 2*pi*2.45e9 * EPS0*substrate_epr; % substrate losses
+
+f_max = 10e9;
+resolution = 250;
+edge_res   = 25;
+feed_shift = 2500;
+meas_shift = 5000;
+
+%% Setup FDTD Parameters & Excitation Function
+FDTD = InitFDTD();
+FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2);
+BC   = {'PML_8' 'PML_8' 'PEC' 'PEC' 'PEC' 'MUR'};
+FDTD = SetBoundaryCond( FDTD, BC );
+
+%% Setup CSXCAD Geometry & Mesh
+CSX = InitCSX();
+edge_mesh  = [-1/3 2/3]*edge_res; % 1/3 - 2/3 rule for 2D metal edges
+
+mesh.x = SmoothMeshLines( [-connect_via_gap 0 connect_via_gap], 2*edge_res, 1.5 );
+mesh.x = SmoothMeshLines( [-line_length mesh.x line_length], resolution, 1.5);
+mesh.y = SmoothMeshLines( [0 msl_width/2+edge_mesh substrate_width/2], resolution/4 , 1.5);
+mesh.y = sort(unique([-mesh.y mesh.y]));
+mesh.z = SmoothMeshLines( [linspace(-strip_substrate_thickness,0,5) linspace(0,strip_substrate_thickness,5) linspace(strip_substrate_thickness,msl_substrate_thickness+strip_substrate_thickness,5) 2*strip_substrate_thickness+air_spacer] , resolution );
+CSX = DefineRectGrid( CSX, unit, mesh );
+
+% Create Substrate
+CSX = AddMaterial( CSX, 'RO4350B' );
+CSX = SetMaterialProperty( CSX, 'RO4350B', 'Epsilon', substrate_epr, 'Kappa', substrate_kappa );
+start = [mesh.x(1),   mesh.y(1),   -strip_substrate_thickness];
+stop  = [mesh.x(end), mesh.y(end), +strip_substrate_thickness+msl_substrate_thickness];
+CSX = AddBox( CSX, 'RO4350B', 0, start, stop );
+
+% Create a PEC called 'metal' and 'gnd'
+CSX = AddMetal( CSX, 'gnd' ); 
+CSX = AddMetal( CSX, 'metal' ); 
+
+% Create strip line port (incl. metal strip line)
+start = [-line_length -strip_width/2  0];
+stop  = [0            +strip_width/2  0];
+[CSX,port{1}] = AddStripLinePort( CSX, 100, 1, 'metal', start, stop, strip_substrate_thickness, 'x', [0 0 -1], 'ExcitePort', true, 'FeedShift', feed_shift, 'MeasPlaneShift', meas_shift );
+
+% Create MSL port on top
+start = [line_length  -strip_width/2 strip_substrate_thickness+msl_substrate_thickness];
+stop  = [0            +strip_width/2 strip_substrate_thickness];
+[CSX,port{2}] = AddMSLPort( CSX, 100, 2, 'metal', start, stop, 'x', [0 0 -1], 'MeasPlaneShift', meas_shift );
+
+% transitional via
+start = [0, 0, 0];
+stop  = [0, 0, strip_substrate_thickness+msl_substrate_thickness];
+CSX = AddCylinder(CSX, 'metal', 100, start, stop, connect_via_rad);
+
+% metal plane between strip line and MSL, including hole for transition
+p(1,1) = mesh.x(1);
+p(2,1) = mesh.y(1);
+p(1,2) = 0;
+p(2,2) = mesh.y(1);
+for a = linspace(-pi, pi, 21)
+  p(1,end+1) = connect_via_gap*sin(a);
+  p(2,end)   = connect_via_gap*cos(a);
+endfor
+p(1,end+1) = 0;
+p(2,end  ) = mesh.y(1);
+p(1,end+1) = mesh.x(end);
+p(2,end  ) = mesh.y(1);
+p(1,end+1) = mesh.x(end);
+p(2,end  ) = mesh.y(end);
+p(1,end+1) = mesh.x(1);
+p(2,end  ) = mesh.y(end);
+CSX = AddPolygon( CSX, 'gnd', 1, 'z', strip_substrate_thickness, p);
+
+%% Write/Show/Run the openEMS compatible xml-file
+Sim_Path = 'tmp';
+Sim_CSX = 'strip2msl.xml';
+
+[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
+[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder
+
+WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
+CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
+RunOpenEMS( Sim_Path, Sim_CSX );
+
+%% Post-Processing
+close all
+f = linspace( 0, f_max, 1601 );
+port = calcPort( port, Sim_Path, f, 'RefImpedance', 50);
+
+s11 = port{1}.uf.ref./ port{1}.uf.inc;
+s21 = port{2}.uf.ref./ port{1}.uf.inc;
+
+plot(f/1e9,20*log10(abs(s11)),'k-','LineWidth',2);
+hold on;
+grid on;
+plot(f/1e9,20*log10(abs(s21)),'r--','LineWidth',2);
+legend('S_{11}','S_{21}');
+ylabel('S-Parameter (dB)','FontSize',12);
+xlabel('frequency (GHz) \rightarrow','FontSize',12);
+ylim([-40 2]);
+
+