diff --git a/FDTD/extensions/engine_ext_tfsf.cpp b/FDTD/extensions/engine_ext_tfsf.cpp
new file mode 100644
index 0000000..db142ba
--- /dev/null
+++ b/FDTD/extensions/engine_ext_tfsf.cpp
@@ -0,0 +1,158 @@
+/*
+* Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+* 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 .
+*/
+
+#include "engine_ext_tfsf.h"
+#include "operator_ext_tfsf.h"
+#include "FDTD/engine_sse.h"
+
+Engine_Ext_TFSF::Engine_Ext_TFSF(Operator_Ext_TFST* op_ext) : Engine_Extension(op_ext)
+{
+ m_Op_TFSF = op_ext;
+ m_Priority = ENG_EXT_PRIO_TFSF;
+
+ m_DelayLookup = new unsigned int[m_Op_TFSF->m_maxDelay+1];
+}
+
+Engine_Ext_TFSF::~Engine_Ext_TFSF()
+{
+ delete[] m_DelayLookup;
+ m_DelayLookup = NULL;
+}
+
+void Engine_Ext_TFSF::DoPostVoltageUpdates()
+{
+ unsigned int numTS = m_Eng->GetNumberOfTimesteps();
+ unsigned int length = m_Op_TFSF->m_Exc->GetLength();
+
+ for (unsigned int n=0;n<=m_Op_TFSF->m_maxDelay;++n)
+ {
+ if ( numTS < n )
+ m_DelayLookup[n]=0;
+ else if ( numTS-n > length)
+ m_DelayLookup[n]=0;
+ else
+ m_DelayLookup[n] = numTS - n;
+ }
+
+ //get the current signal since an H-field is added ...
+ FDTD_FLOAT* signal = m_Op_TFSF->m_Exc->GetCurrentSignal();
+
+ int nP,nPP;
+ unsigned int ui_pos;
+ unsigned int pos[3];
+ for (int n=0;n<3;++n)
+ {
+ nP = (n+1)%3;
+ nPP = (n+2)%3;
+ pos[nP] = m_Op_TFSF->m_Start[nP];
+
+ ui_pos = 0;
+ for (unsigned int i=0;im_numLines[nP];++i)
+ {
+ pos[nPP] = m_Op_TFSF->m_Start[nPP];
+ for (unsigned int j=0;jm_numLines[nPP];++j)
+ {
+ // current updates
+ pos[n] = m_Op_TFSF->m_Start[n];
+
+ m_Eng->SetVolt(nP,pos, m_Eng->GetVolt(nP,pos)
+ + (1.0-m_Op_TFSF->m_VoltDelayDelta[n][0][0][ui_pos])*m_Op_TFSF->m_VoltAmp[n][0][0][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_VoltDelay[n][0][0][ui_pos]]]
+ + m_Op_TFSF->m_VoltDelayDelta[n][0][0][ui_pos] *m_Op_TFSF->m_VoltAmp[n][0][0][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_VoltDelay[n][0][0][ui_pos]]] );
+
+ m_Eng->SetVolt(nPP,pos, m_Eng->GetVolt(nPP,pos)
+ + (1.0-m_Op_TFSF->m_VoltDelayDelta[n][0][1][ui_pos])*m_Op_TFSF->m_VoltAmp[n][0][1][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_VoltDelay[n][0][1][ui_pos]]]
+ + m_Op_TFSF->m_VoltDelayDelta[n][0][1][ui_pos] *m_Op_TFSF->m_VoltAmp[n][0][1][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_VoltDelay[n][0][1][ui_pos]]] );
+
+ pos[n] = m_Op_TFSF->m_Stop[n];
+
+ m_Eng->SetVolt(nP,pos, m_Eng->GetVolt(nP,pos)
+ + (1.0-m_Op_TFSF->m_VoltDelayDelta[n][1][0][ui_pos])*m_Op_TFSF->m_VoltAmp[n][1][0][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_VoltDelay[n][1][0][ui_pos]]]
+ + m_Op_TFSF->m_VoltDelayDelta[n][1][0][ui_pos] *m_Op_TFSF->m_VoltAmp[n][1][0][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_VoltDelay[n][1][0][ui_pos]]] );
+
+ m_Eng->SetVolt(nPP,pos, m_Eng->GetVolt(nPP,pos)
+ + (1.0-m_Op_TFSF->m_VoltDelayDelta[n][1][1][ui_pos])*m_Op_TFSF->m_VoltAmp[n][1][1][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_VoltDelay[n][1][1][ui_pos]]]
+ + m_Op_TFSF->m_VoltDelayDelta[n][1][1][ui_pos] *m_Op_TFSF->m_VoltAmp[n][1][1][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_VoltDelay[n][1][1][ui_pos]]] );
+
+ ++pos[nPP];
+ ++ui_pos;
+ }
+ ++pos[nP];
+ }
+ }
+}
+
+void Engine_Ext_TFSF::DoPostCurrentUpdates()
+{
+ unsigned int numTS = m_Eng->GetNumberOfTimesteps();
+ unsigned int length = m_Op_TFSF->m_Exc->GetLength();
+
+ for (unsigned int n=0;nm_maxDelay;++n)
+ {
+ if ( numTS < n )
+ m_DelayLookup[n]=0;
+ else if ( numTS-n > length)
+ m_DelayLookup[n]=0;
+ else
+ m_DelayLookup[n] = numTS - n;
+ }
+
+ //get the current signal since an E-field is added ...
+ FDTD_FLOAT* signal = m_Op_TFSF->m_Exc->GetVoltageSignal();
+
+ int nP,nPP;
+ unsigned int ui_pos;
+ unsigned int pos[3];
+ for (int n=0;n<3;++n)
+ {
+ nP = (n+1)%3;
+ nPP = (n+2)%3;
+ pos[nP] = m_Op_TFSF->m_Start[nP];
+
+ ui_pos = 0;
+ for (unsigned int i=0;im_numLines[nP];++i)
+ {
+ pos[nPP] = m_Op_TFSF->m_Start[nPP];
+ for (unsigned int j=0;jm_numLines[nPP];++j)
+ {
+ // current updates
+ pos[n] = m_Op_TFSF->m_Start[n]-1;
+
+ m_Eng->SetCurr(nP,pos, m_Eng->GetCurr(nP,pos)
+ + (1.0-m_Op_TFSF->m_CurrDelayDelta[n][0][0][ui_pos])*m_Op_TFSF->m_CurrAmp[n][0][0][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_CurrDelay[n][0][0][ui_pos]]]
+ + m_Op_TFSF->m_CurrDelayDelta[n][0][0][ui_pos] *m_Op_TFSF->m_CurrAmp[n][0][0][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_CurrDelay[n][0][0][ui_pos]]] );
+
+ m_Eng->SetCurr(nPP,pos, m_Eng->GetCurr(nPP,pos)
+ + (1.0-m_Op_TFSF->m_CurrDelayDelta[n][0][1][ui_pos])*m_Op_TFSF->m_CurrAmp[n][0][1][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_CurrDelay[n][0][1][ui_pos]]]
+ + m_Op_TFSF->m_CurrDelayDelta[n][0][1][ui_pos] *m_Op_TFSF->m_CurrAmp[n][0][1][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_CurrDelay[n][0][1][ui_pos]]] );
+
+ pos[n] = m_Op_TFSF->m_Stop[n];
+
+ m_Eng->SetCurr(nP,pos, m_Eng->GetCurr(nP,pos)
+ + (1.0-m_Op_TFSF->m_CurrDelayDelta[n][1][0][ui_pos])*m_Op_TFSF->m_CurrAmp[n][1][0][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_CurrDelay[n][1][0][ui_pos]]]
+ + m_Op_TFSF->m_CurrDelayDelta[n][1][0][ui_pos] *m_Op_TFSF->m_CurrAmp[n][1][0][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_CurrDelay[n][1][0][ui_pos]]] );
+
+ m_Eng->SetCurr(nPP,pos, m_Eng->GetCurr(nPP,pos)
+ + (1.0-m_Op_TFSF->m_CurrDelayDelta[n][1][1][ui_pos])*m_Op_TFSF->m_CurrAmp[n][1][1][ui_pos]*signal[m_DelayLookup[ m_Op_TFSF->m_CurrDelay[n][1][1][ui_pos]]]
+ + m_Op_TFSF->m_CurrDelayDelta[n][1][1][ui_pos] *m_Op_TFSF->m_CurrAmp[n][1][1][ui_pos]*signal[m_DelayLookup[1+m_Op_TFSF->m_CurrDelay[n][1][1][ui_pos]]] );
+
+ ++pos[nPP];
+ ++ui_pos;
+ }
+ ++pos[nP];
+ }
+ }
+}
diff --git a/FDTD/extensions/engine_ext_tfsf.h b/FDTD/extensions/engine_ext_tfsf.h
new file mode 100644
index 0000000..e630379
--- /dev/null
+++ b/FDTD/extensions/engine_ext_tfsf.h
@@ -0,0 +1,40 @@
+/*
+* Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+* 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 .
+*/
+
+#ifndef ENGINE_EXT_TFSF_H
+#define ENGINE_EXT_TFSF_H
+
+#include "engine_extension.h"
+
+class Operator_Ext_TFST;
+
+class Engine_Ext_TFSF : public Engine_Extension
+{
+public:
+ Engine_Ext_TFSF(Operator_Ext_TFST* op_ext);
+ virtual ~Engine_Ext_TFSF();
+
+ virtual void DoPostVoltageUpdates();
+ virtual void DoPostCurrentUpdates();
+
+protected:
+ Operator_Ext_TFST* m_Op_TFSF;
+
+ unsigned int* m_DelayLookup;
+};
+
+#endif // ENGINE_EXT_TFSF_H
diff --git a/FDTD/extensions/engine_extension.h b/FDTD/extensions/engine_extension.h
index 27f7ea5..01d3194 100644
--- a/FDTD/extensions/engine_extension.h
+++ b/FDTD/extensions/engine_extension.h
@@ -23,6 +23,7 @@
// priority definitions for some important extensions
#define ENG_EXT_PRIO_UPML +1e6 //unaxial pml extension priority
#define ENG_EXT_PRIO_CYLINDER +1e5 //cylindrial extension priority
+#define ENG_EXT_PRIO_TFSF +5e4 //total-field/scattered-field extension priority
#define ENG_EXT_PRIO_EXCITATION -1000 //excitation priority
#define ENG_EXT_PRIO_CYLINDERMULTIGRID -3000 //cylindrial multi-grid extension priority
diff --git a/FDTD/extensions/operator_ext_tfsf.cpp b/FDTD/extensions/operator_ext_tfsf.cpp
new file mode 100644
index 0000000..ea3aa0f
--- /dev/null
+++ b/FDTD/extensions/operator_ext_tfsf.cpp
@@ -0,0 +1,381 @@
+/*
+* Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+* 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 .
+*/
+
+#include "operator_ext_tfsf.h"
+#include "engine_ext_tfsf.h"
+#include
+
+Operator_Ext_TFST::Operator_Ext_TFST(Operator* op) : Operator_Extension(op)
+{
+ Init();
+}
+
+Operator_Ext_TFST::~Operator_Ext_TFST()
+{
+ Reset();
+}
+
+void Operator_Ext_TFST::Init()
+{
+ for (int n=0;n<3;++n)
+ for (int l=0;l<2;++l)
+ for (int c=0;c<2;++c)
+ {
+ m_VoltDelay[n][l][c]=NULL;
+ m_VoltDelayDelta[n][l][c]=NULL;
+ m_VoltAmp[n][l][c]=NULL;
+ m_CurrDelay[n][l][c]=NULL;
+ m_CurrDelayDelta[n][l][c]=NULL;
+ m_CurrAmp[n][l][c]=NULL;
+ }
+
+ Operator_Extension::Init();
+}
+
+void Operator_Ext_TFST::Reset()
+{
+ for (int n=0;n<3;++n)
+ for (int l=0;l<2;++l)
+ for (int c=0;c<2;++c)
+ {
+ delete[] m_VoltDelay[n][l][c];
+ m_VoltDelay[n][l][c]=NULL;
+ delete[] m_VoltDelayDelta[n][l][c];
+ m_VoltDelayDelta[n][l][c]=NULL;
+ delete[] m_VoltAmp[n][l][c];
+ m_VoltAmp[n][l][c]=NULL;
+ delete[] m_CurrDelay[n][l][c];
+ m_CurrDelay[n][l][c]=NULL;
+ delete[] m_CurrDelayDelta[n][l][c];
+ m_CurrDelayDelta[n][l][c]=NULL;
+ delete[] m_CurrAmp[n][l][c];
+ m_CurrAmp[n][l][c]=NULL;
+ }
+ Operator_Extension::Reset();
+}
+
+Operator_Extension* Operator_Ext_TFST::Clone(Operator* op)
+{
+ UNUSED(op);
+ return NULL;
+}
+
+bool Operator_Ext_TFST::BuildExtension()
+{
+ m_Exc = m_Op->GetExcitationSignal();
+ double dT = m_Op->GetTimestep();
+ if (dT==0)
+ return false;
+ if (m_Exc==0)
+ return false;
+
+ Reset();
+ ContinuousStructure* CSX = m_Op->GetGeometryCSX();
+
+ vector vec_prop = CSX->GetPropertyByType(CSProperties::EXCITATION);
+
+ if (vec_prop.size()==0)
+ {
+ cerr << "Operator_Ext_TFST::BuildExtension: Warning, no excitation properties found" << endl;
+ SetActive(false);
+ return false;
+ }
+
+ CSPropExcitation* elec=NULL;
+ CSProperties* prop=NULL;
+ CSPrimitives* prim=NULL;
+ CSPrimBox* box=NULL;
+ for (size_t p=0; pToExcitation();
+ if (elec->GetExcitType()!=10)
+ continue;
+ if (prop->GetQtyPrimitives()!=1)
+ {
+ cerr << "Operator_Ext_TFST::BuildExtension: Warning, plane wave excitation found with more (or less) than one primitive, skipping..." << endl;
+ continue;
+ }
+ prim = prop->GetPrimitive(0);
+ if (prim->GetType()!=CSPrimitives::BOX)
+ {
+ cerr << "Operator_Ext_TFST::BuildExtension: Warning, plane wave excitation found with false non-Box primitive, skipping..." << endl;
+ continue;
+ }
+ box = prim->ToBox();
+ if (box==NULL) //sanity check, should not happen!
+ {
+ SetActive(false);
+ return false;
+ }
+
+ // found a plane-wave excite with exactly one box
+ for (int n=0;n<3;++n)
+ m_PropDir[n] = elec->GetPropagationDir(n);
+ double dir_norm = sqrt(m_PropDir[0]*m_PropDir[0]+m_PropDir[1]*m_PropDir[1]+m_PropDir[2]*m_PropDir[2]);
+ if (dir_norm==0)
+ {
+ cerr << "Operator_Ext_TFST::BuildExtension: Warning, plane wave direction is zero, skipping..." << endl;
+ SetActive(false);
+ return false;
+ }
+
+ //make it a unit vector
+ m_PropDir[0]/=dir_norm;m_PropDir[1]/=dir_norm;m_PropDir[2]/=dir_norm;
+
+ if (m_Op->SnapBox2Mesh(box->GetStartCoord()->GetNativeCoords(), box->GetStopCoord()->GetNativeCoords(), m_Start, m_Stop)!=3)
+ {
+ cerr << "Operator_Ext_TFST::BuildExtension: Warning, plane wave box dimension is invalid, skipping..." << endl;
+ SetActive(false);
+ return false;
+ }
+
+ double origin[3];
+ unsigned int ui_origin[3];
+ for (int n=0;n<3;++n)
+ {
+ m_E_Amp[n] = elec->GetExcitation(n);
+ m_numLines[n] = m_Stop[n]-m_Start[n]+1;
+ m_IncLow[n] = m_PropDir[n]>=0;
+
+ if (m_IncLow[n])
+ {
+ ui_origin[n] = m_Start[n];
+ }
+ else
+ {
+ ui_origin[n] = m_Stop[n];
+ }
+ origin[n] = m_Op->GetDiscLine(n,ui_origin[n]);
+ }
+
+ double dotEk = (m_E_Amp[0]*m_PropDir[0] + m_E_Amp[1]*m_PropDir[1] + m_E_Amp[2]*m_PropDir[2]);
+ double angle = acos( dotEk / (m_E_Amp[0]*m_E_Amp[0] + m_E_Amp[1]*m_E_Amp[1] + m_E_Amp[2]*m_E_Amp[2]) ) / PI * 180;
+
+ if (angle==0)
+ {
+ cerr << "Operator_Ext_TFST::BuildExtension: Warning, plane wave direction and polarization is identical, skipping..." << endl;
+ SetActive(false);
+ return false;
+ }
+ if (angle!=90)
+ {
+ cerr << "Operator_Ext_TFST::BuildExtension: Warning, angle between propagation direction and polarization is not 90deg, resetting E-polarization to : (";
+ for (int n=0;n<3;++n)
+ m_E_Amp[n]-=m_PropDir[n]*dotEk;
+ cerr << m_E_Amp[0] << "," << m_E_Amp[1] << "," << m_E_Amp[2] << ")" << endl;
+ }
+
+ int nP,nPP;
+ for (int n=0;n<3;++n)
+ {
+ nP = (n+1)%3;
+ nPP = (n+2)%3;
+ m_H_Amp[n] = m_PropDir[nP]*m_E_Amp[nPP] - m_PropDir[nPP]*m_E_Amp[nP];
+ m_H_Amp[n] /= __Z0__;
+ }
+
+ double coord[3];
+ double unit = m_Op->GetGridDelta();
+ double delay;
+ double dist;
+ unsigned int pos[3];
+ unsigned int numP, ui_pos;
+ m_maxDelay = 0;
+ for (int n=0;n<3;++n)
+ {
+ nP = (n+1)%3;
+ nPP = (n+2)%3;
+ pos[n] = 0; //unused
+ pos[nP] = m_Start[nP];
+
+ numP = m_numLines[nP]*m_numLines[nPP];
+
+ for (int l=0;l<2;++l)
+ for (int c=0;c<2;++c)
+ {
+ m_VoltDelay[n][l][c]=new unsigned int[numP];
+ m_VoltDelayDelta[n][l][c]=new FDTD_FLOAT[numP];
+ m_VoltAmp[n][l][c]=new FDTD_FLOAT[numP];
+
+ m_CurrDelay[n][l][c]=new unsigned int[numP];
+ m_CurrDelayDelta[n][l][c]=new FDTD_FLOAT[numP];
+ m_CurrAmp[n][l][c]=new FDTD_FLOAT[numP];
+ }
+
+ ui_pos = 0;
+ for (unsigned int i=0;iGetYeeCoords(nP,pos,coord,false);
+ dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
+ delay = dist*unit/__C0__/dT;
+ m_maxDelay = max((unsigned int)delay,m_maxDelay);
+ m_CurrDelay[n][0][1][ui_pos] = floor(delay);
+ m_CurrDelayDelta[n][0][1][ui_pos] = delay - floor(delay);
+ m_CurrAmp[n][0][1][ui_pos] = m_E_Amp[nP]*m_Op->GetEdgeLength(nP,pos);
+
+ m_Op->GetYeeCoords(nPP,pos,coord,false);
+ dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
+ delay = dist*unit/__C0__/dT;
+ m_maxDelay = max((unsigned int)delay,m_maxDelay);
+ m_CurrDelay[n][0][0][ui_pos] = floor(delay);
+ m_CurrDelayDelta[n][0][0][ui_pos] = delay - floor(delay);
+ m_CurrAmp[n][0][0][ui_pos] = m_E_Amp[nPP]*m_Op->GetEdgeLength(nPP,pos);
+
+ --pos[n];
+ m_CurrAmp[n][0][0][ui_pos]*=m_Op->GetIV(nP,pos);
+ m_CurrAmp[n][0][1][ui_pos]*=m_Op->GetIV(nPP,pos);
+
+ pos[n] = m_Stop[n];
+ m_Op->GetYeeCoords(nP,pos,coord,false);
+ dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
+ delay = dist*unit/__C0__/dT;
+ m_maxDelay = max((unsigned int)delay,m_maxDelay);
+ m_CurrDelay[n][1][1][ui_pos] = floor(delay);
+ m_CurrDelayDelta[n][1][1][ui_pos] = delay - floor(delay);
+ m_CurrAmp[n][1][1][ui_pos] = m_E_Amp[nP]*m_Op->GetEdgeLength(nP,pos);
+
+ m_Op->GetYeeCoords(nPP,pos,coord,false);
+ dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
+ delay = dist*unit/__C0__/dT;
+ m_maxDelay = max((unsigned int)delay,m_maxDelay);
+ m_CurrDelay[n][1][0][ui_pos] = floor(delay);
+ m_CurrDelayDelta[n][1][0][ui_pos] = delay - floor(delay);
+ m_CurrAmp[n][1][0][ui_pos] = m_E_Amp[nPP]*m_Op->GetEdgeLength(nPP,pos);
+
+ m_CurrAmp[n][1][0][ui_pos]*=m_Op->GetIV(nP,pos);
+ m_CurrAmp[n][1][1][ui_pos]*=m_Op->GetIV(nPP,pos);
+
+ if (m_IncLow[n])
+ {
+ m_CurrAmp[n][0][0][ui_pos]*=-1;
+ m_CurrAmp[n][1][1][ui_pos]*=-1;
+ }
+ else
+ {
+ m_CurrAmp[n][0][1][ui_pos]*=-1;
+ m_CurrAmp[n][1][0][ui_pos]*=-1;
+ }
+
+ if (pos[nP]==m_Stop[nP])
+ {
+ m_CurrAmp[n][0][1][ui_pos]=0;
+ m_CurrAmp[n][1][1][ui_pos]=0;
+ }
+ if (pos[nPP]==m_Stop[nPP])
+ {
+ m_CurrAmp[n][0][0][ui_pos]=0;
+ m_CurrAmp[n][1][0][ui_pos]=0;
+ }
+
+ // voltage updates
+ pos[n] = m_Start[n]-1;
+ m_Op->GetYeeCoords(nP,pos,coord,true);
+ dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
+ delay = dist*unit/__C0__/dT;
+ m_maxDelay = max((unsigned int)delay,m_maxDelay);
+ m_VoltDelay[n][0][1][ui_pos] = floor(delay);
+ m_VoltDelayDelta[n][0][1][ui_pos] = delay - floor(delay);
+ m_VoltAmp[n][0][1][ui_pos] = m_H_Amp[nP]*m_Op->GetEdgeLength(nP,pos,true);
+
+ m_Op->GetYeeCoords(nPP,pos,coord,true);
+ dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
+ delay = dist*unit/__C0__/dT;
+ m_maxDelay = max((unsigned int)delay,m_maxDelay);
+ m_VoltDelay[n][0][0][ui_pos] = floor(delay);
+ m_VoltDelayDelta[n][0][0][ui_pos] = delay - floor(delay);
+ m_VoltAmp[n][0][0][ui_pos] = m_H_Amp[nPP]*m_Op->GetEdgeLength(nPP,pos,true);
+
+ ++pos[n];
+ m_VoltAmp[n][0][0][ui_pos]*=m_Op->GetVI(nP,pos);
+ m_VoltAmp[n][0][1][ui_pos]*=m_Op->GetVI(nPP,pos);
+
+ pos[n] = m_Stop[n];
+ m_Op->GetYeeCoords(nP,pos,coord,true);
+ dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
+ delay = dist*unit/__C0__/dT;
+ m_maxDelay = max((unsigned int)delay,m_maxDelay);
+ m_VoltDelay[n][1][1][ui_pos] = floor(delay);
+ m_VoltDelayDelta[n][1][1][ui_pos] = delay - floor(delay);
+ m_VoltAmp[n][1][1][ui_pos] = m_H_Amp[nP]*m_Op->GetEdgeLength(nP,pos,true);
+
+ m_Op->GetYeeCoords(nPP,pos,coord,true);
+ dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
+ delay = dist*unit/__C0__/dT;
+ m_maxDelay = max((unsigned int)delay,m_maxDelay);
+ m_VoltDelay[n][1][0][ui_pos] = floor(delay);
+ m_VoltDelayDelta[n][1][0][ui_pos] = delay - floor(delay);
+ m_VoltAmp[n][1][0][ui_pos] = m_H_Amp[nPP]*m_Op->GetEdgeLength(nPP,pos,true);
+
+ m_VoltAmp[n][1][0][ui_pos]*=m_Op->GetVI(nP,pos);
+ m_VoltAmp[n][1][1][ui_pos]*=m_Op->GetVI(nPP,pos);
+
+ if (m_IncLow[n])
+ {
+ m_VoltAmp[n][1][0][ui_pos]*=-1;
+ m_VoltAmp[n][0][1][ui_pos]*=-1;
+ }
+ else
+ {
+ m_VoltAmp[n][0][0][ui_pos]*=-1;
+ m_VoltAmp[n][1][1][ui_pos]*=-1;
+ }
+
+ if (pos[nP]==m_Stop[nP])
+ {
+ m_VoltAmp[n][0][0][ui_pos]=0;
+ m_VoltAmp[n][1][0][ui_pos]=0;
+ }
+ if (pos[nPP]==m_Stop[nPP])
+ {
+ m_VoltAmp[n][0][1][ui_pos]=0;
+ m_VoltAmp[n][1][1][ui_pos]=0;
+ }
+
+ ++pos[nPP];
+ ++ui_pos;
+ }
+ ++pos[nP];
+ }
+ }
+ ++m_maxDelay;
+ return true;
+ }
+ SetActive(false);
+ return false;
+}
+
+Engine_Extension* Operator_Ext_TFST::CreateEngineExtention()
+{
+ return new Engine_Ext_TFSF(this);
+}
+
+void Operator_Ext_TFST::ShowStat(ostream &ostr) const
+{
+ Operator_Extension::ShowStat(ostr);
+ cout << "Propagation direction\t: " << "(" << m_PropDir[0] << ", " << m_PropDir[1] << ", " << m_PropDir[2] << ")" << endl;
+ cout << "E-field amplitude\t: " << "(" << m_E_Amp[0] << ", " << m_E_Amp[1] << ", " << m_E_Amp[2] << ")" << endl;
+ cout << "m_H_Amp amplitude\t: " << "(" << m_H_Amp[0]*__Z0__ << ", " << m_H_Amp[1]*__Z0__ << ", " << m_H_Amp[2]*__Z0__ << ")/Z0" << endl;
+ cout << "Box Dimensions\t\t: " << m_numLines[0] << " x " << m_numLines[1] << " x " << m_numLines[2] << endl;
+ cout << "Max. Delay (TS)\t\t: " << m_maxDelay << endl;
+ cout << "Memory usage (est.)\t: ~" << m_numLines[0] * m_numLines[1] * m_numLines[2] * 6 * 4 * 4 / 1024 << " kiB" << endl;
+}
diff --git a/FDTD/extensions/operator_ext_tfsf.h b/FDTD/extensions/operator_ext_tfsf.h
new file mode 100644
index 0000000..4f42b64
--- /dev/null
+++ b/FDTD/extensions/operator_ext_tfsf.h
@@ -0,0 +1,75 @@
+/*
+* Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+* 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 .
+*/
+
+#ifndef OPERATOR_EXT_TFSF_H
+#define OPERATOR_EXT_TFSF_H
+
+#include "operator_extension.h"
+#include "FDTD/operator.h"
+#include "tools/constants.h"
+
+class Excitation;
+
+class Operator_Ext_TFST : public Operator_Extension
+{
+ friend class Engine_Ext_TFSF;
+public:
+ Operator_Ext_TFST(Operator* op);
+ ~Operator_Ext_TFST();
+
+ 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 false;}
+ virtual bool IsCylindricalMultiGridSave(bool child) const {UNUSED(child); return false;}
+
+ virtual string GetExtensionName() const {return string("Total-Field/Scattered-Field Extension");}
+
+ virtual void ShowStat(ostream &ostr) const;
+
+ virtual void Init();
+ virtual void Reset();
+
+protected:
+ Excitation* m_Exc;
+
+ bool m_IncLow[3];
+ unsigned int m_Start[3];
+ unsigned int m_Stop[3];
+ unsigned int m_numLines[3];
+
+ double m_PropDir[3];
+ double m_E_Amp[3];
+ double m_H_Amp[3];
+
+ unsigned int m_maxDelay;
+
+ // array setup [direction][low/high][component][ ]
+ unsigned int* m_VoltDelay[3][2][2];
+ FDTD_FLOAT* m_VoltDelayDelta[3][2][2];
+ FDTD_FLOAT* m_VoltAmp[3][2][2];
+
+ unsigned int* m_CurrDelay[3][2][2];
+ FDTD_FLOAT* m_CurrDelayDelta[3][2][2];
+ FDTD_FLOAT* m_CurrAmp[3][2][2];
+
+};
+
+#endif // OPERATOR_EXT_TFSF_H
diff --git a/matlab/Tutorials/RCS_Sphere.m b/matlab/Tutorials/RCS_Sphere.m
new file mode 100644
index 0000000..ce8d0f6
--- /dev/null
+++ b/matlab/Tutorials/RCS_Sphere.m
@@ -0,0 +1,89 @@
+%
+% Tutorials / radar cross section of a metal sphere
+%
+% This tutorial is unfinished!
+%
+% Describtion at:
+% http://openems.de/index.php/Tutorial:_RCS_Sphere
+%
+% Tested with
+% - Matlab 2011a / Octave 3.4.3
+% - openEMS v0.0.29
+%
+% (C) 2012 Thorsten Liebig
+
+close all
+clear
+clc
+
+%% setup the simulation
+physical_constants;
+unit = 1e-3; % all length in mm
+
+sphere.rad = 160;
+
+% size of the simulation box
+SimBox = 1000;
+PW_Box = 500;
+
+%% setup FDTD parameter & excitation function
+f_start = 200e6; % start frequency
+f_stop = 1000e6; % stop frequency
+
+FDTD = InitFDTD( 30000 );
+FDTD = SetGaussExcite( FDTD, 0.5*(f_start+f_stop), 0.5*(f_stop-f_start) );
+BC = [1 1 1 1 1 1]*3; % set boundary conditions
+FDTD = SetBoundaryCond( FDTD, BC );
+
+%% setup CSXCAD geometry & mesh
+% currently, openEMS cannot automatically generate a mesh
+max_res = c0 / f_stop / unit / 20; % cell size: lambda/20
+CSX = InitCSX();
+
+%create fixed lines for the simulation box, substrate and port
+mesh.x = SmoothMeshLines([-SimBox/2 0 SimBox/2], max_res);
+mesh.y = mesh.x;
+mesh.z = mesh.x;
+
+%% create metal sphere
+CSX = AddMetal( CSX, 'sphere' ); % create a perfect electric conductor (PEC)
+CSX = AddSphere(CSX,'sphere',10,[0 0 0],sphere.rad);
+
+%% plane wave excitation
+k_dir = [1 0 0]; % plane wave direction --> x-direction
+E_dir = [0 0 1]; % plane wave polarization --> E_z
+
+CSX = AddPlaneWaveExcite(CSX, 'plane_wave', k_dir, E_dir);
+start = [-PW_Box/2 -PW_Box/2 -PW_Box/2];
+stop = -start;
+CSX = AddBox(CSX, 'plane_wave', 0, start, stop);
+
+%% dump boxes
+CSX = AddDump(CSX, 'Et');
+start = [mesh.x(1) mesh.y(1) 0];
+stop = [mesh.x(end) mesh.y(end) 0];
+CSX = AddBox(CSX, 'Et', 0, start, stop);
+
+% add 8 lines in all direction as pml spacing
+mesh = AddPML(mesh,8);
+
+CSX = DefineRectGrid( CSX, unit, mesh );
+
+%% prepare simulation folder
+Sim_Path = 'Sphere_RCS';
+Sim_CSX = 'Sphere_RCS.xml';
+
+[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
+WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
+
+%% show the structure
+CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
+
+%% run openEMS
+RunOpenEMS( Sim_Path, Sim_CSX);
+
+%%
+disp('Use Paraview to display the elctric fields dumped by openEMS');
diff --git a/openEMS.pro b/openEMS.pro
index b4b2ac6..e40bf93 100644
--- a/openEMS.pro
+++ b/openEMS.pro
@@ -123,7 +123,9 @@ SOURCES += FDTD/extensions/engine_extension.cpp \
FDTD/extensions/operator_ext_cylinder.cpp \
FDTD/extensions/engine_ext_cylinder.cpp \
FDTD/extensions/operator_ext_excitation.cpp \
- FDTD/extensions/engine_ext_excitation.cpp
+ FDTD/extensions/engine_ext_excitation.cpp \
+ FDTD/extensions/operator_ext_tfsf.cpp \
+ FDTD/extensions/engine_ext_tfsf.cpp
# Common source files
SOURCES += Common/operator_base.cpp \
@@ -188,7 +190,9 @@ HEADERS += FDTD/extensions/operator_extension.h \
FDTD/extensions/operator_ext_upml.h \
FDTD/extensions/engine_ext_upml.h \
FDTD/extensions/operator_ext_excitation.h \
- FDTD/extensions/engine_ext_excitation.h
+ FDTD/extensions/engine_ext_excitation.h \
+ FDTD/extensions/operator_ext_tfsf.h \
+ FDTD/extensions/engine_ext_tfsf.h
# Common header files
HEADERS += Common/operator_base.h \
diff --git a/openems.cpp b/openems.cpp
index 109814a..6eeb683 100644
--- a/openems.cpp
+++ b/openems.cpp
@@ -24,6 +24,7 @@
#include "FDTD/engine_multithread.h"
#include "FDTD/operator_multithread.h"
#include "FDTD/extensions/operator_ext_excitation.h"
+#include "FDTD/extensions/operator_ext_tfsf.h"
#include "FDTD/extensions/operator_ext_mur_abc.h"
#include "FDTD/extensions/operator_ext_pml_sf.h"
#include "FDTD/extensions/operator_ext_upml.h"
@@ -631,6 +632,7 @@ int openEMS::SetupFDTD(const char* file)
m_Exc = new Excitation();
FDTD_Op->SetExcitationSignal(m_Exc);
FDTD_Op->AddExtension(new Operator_Ext_Excitation(FDTD_Op));
+ FDTD_Op->AddExtension(new Operator_Ext_TFST(FDTD_Op));
if (FDTD_Op->SetGeometryCSX(m_CSX)==false) return(2);