From be8a3fbc51eb95d6a1b06f7731c9c972b33a1a02 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Wed, 18 Jul 2012 13:12:25 +0200 Subject: [PATCH] new extension: total-field/scattered-field excitation See new matlab tutorial: RCS_Sphere.m Signed-off-by: Thorsten Liebig --- FDTD/extensions/engine_ext_tfsf.cpp | 158 +++++++++++ FDTD/extensions/engine_ext_tfsf.h | 40 +++ FDTD/extensions/engine_extension.h | 1 + FDTD/extensions/operator_ext_tfsf.cpp | 381 ++++++++++++++++++++++++++ FDTD/extensions/operator_ext_tfsf.h | 75 +++++ matlab/Tutorials/RCS_Sphere.m | 89 ++++++ openEMS.pro | 8 +- openems.cpp | 2 + 8 files changed, 752 insertions(+), 2 deletions(-) create mode 100644 FDTD/extensions/engine_ext_tfsf.cpp create mode 100644 FDTD/extensions/engine_ext_tfsf.h create mode 100644 FDTD/extensions/operator_ext_tfsf.cpp create mode 100644 FDTD/extensions/operator_ext_tfsf.h create mode 100644 matlab/Tutorials/RCS_Sphere.m 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);