diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 4891f8c..ed8904a 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -393,15 +393,51 @@ void Operator::DumpPEC2File( string filename ) double scaling = 1; #endif - for (pos[0]=0; pos[0]SetPos(pos[0],pos[1],pos[2]); if (EC_C[n][i]>0) { - GetVV(n,pos[0],pos[1],pos[2]) = (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]); - GetVI(n,pos[0],pos[1],pos[2]) = (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]); + SetVV(n,pos[0],pos[1],pos[2], (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]) ); + SetVI(n,pos[0],pos[1],pos[2], (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]) ); } else { - GetVV(n,pos[0],pos[1],pos[2]) = 0; - GetVI(n,pos[0],pos[1],pos[2]) = 0; + SetVV(n,pos[0],pos[1],pos[2], 0 ); + SetVI(n,pos[0],pos[1],pos[2], 0 ); } if (EC_L[n][i]>0) { - GetII(n,pos[0],pos[1],pos[2]) = (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]); - GetIV(n,pos[0],pos[1],pos[2]) = (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]); + SetII(n,pos[0],pos[1],pos[2], (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]) ); + SetIV(n,pos[0],pos[1],pos[2], (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]) ); } else { - GetII(n,pos[0],pos[1],pos[2]) = 0; - GetIV(n,pos[0],pos[1],pos[2]) = 0; + SetII(n,pos[0],pos[1],pos[2], 0 ); + SetIV(n,pos[0],pos[1],pos[2], 0 ); } } @@ -632,19 +668,19 @@ void Operator::ApplyElectricBC(bool* dirs) for (pos[nPP]=0;pos[nPP]GetExcitType()==1) //hard excite { - GetVV(n,pos[0],pos[1],pos[2]) = 0; - GetVI(n,pos[0],pos[1],pos[2]) = 0; + SetVV(n,pos[0],pos[1],pos[2], 0 ); + SetVI(n,pos[0],pos[1],pos[2], 0 ); } } } @@ -1203,8 +1239,8 @@ bool Operator::CalcFieldExcitation() } if (elec->GetExcitType()==3) //hard excite { - GetII(n,pos[0],pos[1],pos[2]) = 0; - GetIV(n,pos[0],pos[1],pos[2]) = 0; + SetII(n,pos[0],pos[1],pos[2], 0 ); + SetIV(n,pos[0],pos[1],pos[2], 0 ); } } } @@ -1266,8 +1302,8 @@ bool Operator::CalcFieldExcitation() } if (elec->GetExcitType()==1) //hard excite { - GetVV(n,pos[0],pos[1],pos[2]) = 0; - GetVI(n,pos[0],pos[1],pos[2]) = 0; + SetVV(n,pos[0],pos[1],pos[2], 0 ); + SetVI(n,pos[0],pos[1],pos[2], 0 ); } } } @@ -1321,8 +1357,8 @@ void Operator::CalcPEC_Range(unsigned int startX, unsigned int stopX, unsigned i { if (prop->GetType()==CSProperties::METAL) //set to PEC { - GetVV(n,pos[0],pos[1],pos[2]) = 0; - GetVI(n,pos[0],pos[1],pos[2]) = 0; + SetVV(n,pos[0],pos[1],pos[2], 0 ); + SetVI(n,pos[0],pos[1],pos[2], 0 ); ++counter[n]; // cerr << "CartOperator::CalcPEC: PEC found at " << pos[0] << " ; " << pos[1] << " ; " << pos[2] << endl; } @@ -1359,8 +1395,8 @@ void Operator::CalcPEC_Curves() for (size_t t=0;tsetupExcitation(Excite,maxTS);}; - inline virtual FDTD_FLOAT& GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vv[n][x][y][z]; } - inline virtual FDTD_FLOAT& GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vi[n][x][y][z]; } + // the next four functions need to be reimplemented in a derived class + inline virtual FDTD_FLOAT GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vv[n][x][y][z]; } + inline virtual FDTD_FLOAT GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vi[n][x][y][z]; } + inline virtual FDTD_FLOAT GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return ii[n][x][y][z]; } + inline virtual FDTD_FLOAT GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return iv[n][x][y][z]; } - inline virtual FDTD_FLOAT& GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return ii[n][x][y][z]; } - inline virtual FDTD_FLOAT& GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return iv[n][x][y][z]; } + // the next four functions need to be reimplemented in a derived class + inline virtual void SetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { vv[n][x][y][z] = value; } + inline virtual void SetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { vi[n][x][y][z] = value; } + inline virtual void SetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { ii[n][x][y][z] = value; } + inline virtual void SetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { iv[n][x][y][z] = value; } virtual void SetBoundaryCondition(int* BCs) {for (int n=0;n<6;++n) m_BC[n]=BCs[n];} virtual void ApplyElectricBC(bool* dirs); //applied by default to all boundaries diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp index ccf9b79..385b29a 100644 --- a/FDTD/operator_cylinder.cpp +++ b/FDTD/operator_cylinder.cpp @@ -207,8 +207,8 @@ void Operator_Cylinder::ApplyElectricBC(bool* dirs) { for (pos[2]=0;pos[2]GetVV(2,0,0,pos[2]) = 1; + m_Op->SetVV(2,0,0,pos[2], 1); vv_R0[pos[2]] = (1-dT*G/2/C)/(1+dT*G/2/C); vi_R0[pos[2]] = (dT/C)/(1+dT*G/2/C); } diff --git a/FDTD/operator_ext_upml.cpp b/FDTD/operator_ext_upml.cpp index 0bfc575..f85f693 100644 --- a/FDTD/operator_ext_upml.cpp +++ b/FDTD/operator_ext_upml.cpp @@ -378,8 +378,8 @@ bool Operator_Ext_UPML::BuildExtension() //modify the original operator to perform eq. (7.85) by the main engine (EC-FDTD: equation is multiplied by delta_n) //the engine extension will replace the original voltages with the "voltage flux" (volt*eps0) prior to the voltage updates //after the updates are done the extension will calculate the new voltages (see below) and place them back into the main field domain - m_Op->GetVV(n,pos[0],pos[1],pos[2]) = (2*__EPS0__ - kappa_v[nP]*dT) / (2*__EPS0__ + kappa_v[nP]*dT); - m_Op->GetVI(n,pos[0],pos[1],pos[2]) = (2*__EPS0__*dT) / (2*__EPS0__ + kappa_v[nP]*dT) * m_Op->GetEdgeLength(n,pos) / m_Op->GetEdgeArea(n,pos); + m_Op->SetVV(n,pos[0],pos[1],pos[2], (2*__EPS0__ - kappa_v[nP]*dT) / (2*__EPS0__ + kappa_v[nP]*dT) ); + m_Op->SetVI(n,pos[0],pos[1],pos[2], (2*__EPS0__*dT) / (2*__EPS0__ + kappa_v[nP]*dT) * m_Op->GetEdgeLength(n,pos) / m_Op->GetEdgeArea(n,pos) ); //operators needed by eq. (7.88) to calculate new voltages from old voltages and old and new "voltage fluxes" @@ -392,7 +392,7 @@ bool Operator_Ext_UPML::BuildExtension() { //disable upml GetVV(n,loc_pos) = m_Op->GetVV(n,pos[0],pos[1],pos[2]); - m_Op->GetVV(n,pos[0],pos[1],pos[2]) = 0; + m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 ); GetVVFO(n,loc_pos) = 0; GetVVFN(n,loc_pos) = 1; } @@ -405,8 +405,8 @@ bool Operator_Ext_UPML::BuildExtension() //modify the original operator to perform eq. (7.89) by the main engine (EC-FDTD: equation is multiplied by delta_n) //the engine extension will replace the original currents with the "current flux" (curr*mu0) prior to the current updates //after the updates are done the extension will calculate the new currents (see below) and place them back into the main field domain - m_Op->GetII(n,pos[0],pos[1],pos[2]) = (2*__EPS0__ - kappa_i[nP]*dT) / (2*__EPS0__ + kappa_i[nP]*dT); - m_Op->GetIV(n,pos[0],pos[1],pos[2]) = (2*__EPS0__*dT) / (2*__EPS0__ + kappa_i[nP]*dT) * m_Op->GetEdgeLength(n,pos,true) / m_Op->GetEdgeArea(n,pos,true); + m_Op->SetII(n,pos[0],pos[1],pos[2], (2*__EPS0__ - kappa_i[nP]*dT) / (2*__EPS0__ + kappa_i[nP]*dT) ); + m_Op->SetIV(n,pos[0],pos[1],pos[2], (2*__EPS0__*dT) / (2*__EPS0__ + kappa_i[nP]*dT) * m_Op->GetEdgeLength(n,pos,true) / m_Op->GetEdgeArea(n,pos,true) ); //operators needed by eq. (7.90) to calculate new currents from old currents and old and new "current fluxes" GetII(n,loc_pos) = (2*__EPS0__ - kappa_i[nPP]*dT) / (2*__EPS0__ + kappa_i[nPP]*dT); @@ -418,7 +418,7 @@ bool Operator_Ext_UPML::BuildExtension() { //disable upml GetII(n,loc_pos) = m_Op->GetII(n,pos[0],pos[1],pos[2]); - m_Op->GetII(n,pos[0],pos[1],pos[2]) = 0; + m_Op->SetII(n,pos[0],pos[1],pos[2], 0 ); GetIIFO(n,loc_pos) = 0; GetIIFN(n,loc_pos) = 1; } diff --git a/FDTD/operator_sse.h b/FDTD/operator_sse.h index 7d9687a..d86d089 100644 --- a/FDTD/operator_sse.h +++ b/FDTD/operator_sse.h @@ -32,11 +32,15 @@ public: virtual void DumpOperator2File(string filename); - 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]; } + 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]; } + inline virtual FDTD_FLOAT GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_ii[n][x][y][z%numVectors].f[z/numVectors]; } + inline virtual FDTD_FLOAT GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_iv[n][x][y][z%numVectors].f[z/numVectors]; } - inline virtual FDTD_FLOAT& GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_ii[n][x][y][z%numVectors].f[z/numVectors]; } - inline virtual FDTD_FLOAT& GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_iv[n][x][y][z%numVectors].f[z/numVectors]; } + inline virtual void SetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { f4_vv[n][x][y][z%numVectors].f[z/numVectors] = value; } + inline virtual void SetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { f4_vi[n][x][y][z%numVectors].f[z/numVectors] = value; } + inline virtual void SetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { f4_ii[n][x][y][z%numVectors].f[z/numVectors] = value; } + inline virtual void SetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { f4_iv[n][x][y][z%numVectors].f[z/numVectors] = value; } protected: //! use New() for creating a new Operator diff --git a/FDTD/operator_sse_compressed.h b/FDTD/operator_sse_compressed.h index 2928a33..6aeb29e 100644 --- a/FDTD/operator_sse_compressed.h +++ b/FDTD/operator_sse_compressed.h @@ -49,11 +49,15 @@ public: virtual int CalcECOperator(); - inline virtual FDTD_FLOAT& GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { 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 ) { 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);} + 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);} + inline virtual FDTD_FLOAT GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_ii_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetII(n,x,y,z);} + inline virtual FDTD_FLOAT GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_iv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetIV(n,x,y,z);} - inline virtual FDTD_FLOAT& GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { if (m_Use_Compression) return f4_ii_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetII(n,x,y,z);} - inline virtual FDTD_FLOAT& GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { if (m_Use_Compression) return f4_iv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetIV(n,x,y,z);} + inline virtual void SetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { if (m_Use_Compression) f4_vv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors] = value; else Operator_sse::SetVV(n,x,y,z,value);} + inline virtual void SetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { if (m_Use_Compression) f4_vi_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors] = value; else Operator_sse::SetVI(n,x,y,z,value);} + inline virtual void SetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { if (m_Use_Compression) f4_ii_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors] = value; else Operator_sse::SetII(n,x,y,z,value);} + inline virtual void SetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { if (m_Use_Compression) f4_iv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors] = value; else Operator_sse::SetIV(n,x,y,z,value);} virtual void ShowStat() const; diff --git a/FDTD/process_efield.cpp b/FDTD/process_efield.cpp index 464c5e3..b1c189c 100644 --- a/FDTD/process_efield.cpp +++ b/FDTD/process_efield.cpp @@ -111,7 +111,7 @@ int ProcessEField::Process() field *= m_weight; for (size_t n=0;n typedef std::complex double_complex; -#define II double_complex(0.0,1.0) +#define _I double_complex(0.0,1.0) #include #include diff --git a/FDTD/processintegral.cpp b/FDTD/processintegral.cpp index a690b65..cc2321a 100644 --- a/FDTD/processintegral.cpp +++ b/FDTD/processintegral.cpp @@ -77,7 +77,7 @@ int ProcessIntegral::Process() double T = time; for (size_t n=0;n= 8) && (~isempty(refplaneshift)) % refplaneshift counts from start of port diff --git a/matlab/calcMSLPort.m b/matlab/calcMSLPort.m index d76b3cd..5f6cc40 100644 --- a/matlab/calcMSLPort.m +++ b/matlab/calcMSLPort.m @@ -12,7 +12,7 @@ function [S11,beta,ZL] = calcMSLPort( portstruct, SimDir, f, ref_shift ) % ref_shift: (optional) reference plane shift measured from start of port (in drawing units) % % output: -% S11: reflection coefficient +% S11: reflection coefficient (normalized to ZL) % beta: propagation constant % ZL: characteristic line impedance % @@ -21,7 +21,7 @@ function [S11,beta,ZL] = calcMSLPort( portstruct, SimDir, f, ref_shift ) % % openEMS matlab interface % ----------------------- -% Sebastian Held +% (C) 2010 Sebastian Held % See also AddMSLPort % check @@ -90,15 +90,10 @@ S11 = (-1i * dEt + Et .* temp) ./ (Et .* temp + 1i * dEt); % solution 1 % determine ZL ZL = sqrt(Et .* dEt ./ (Ht .* dHt)); -% reference plane shift +% reference plane shift (lossless) if (nargin > 3) % renormalize the shift to the measurement plane - if (portstruct.stop(portstruct.idx_prop) - portstruct.start(portstruct.idx_prop) > 0) - dir = +1; - else - dir = -1; - end - ref_shift = ref_shift - dir*(portstruct.v2_start(portstruct.idx_prop) - portstruct.start(portstruct.idx_prop)); + ref_shift = ref_shift - portstruct.measplanepos; ref_shift = ref_shift * portstruct.drawingunit; S11 = S11 .* exp(2i*real(beta)*ref_shift); S11_corrected = S11_corrected .* exp(2i*real(beta)*ref_shift); diff --git a/matlab/examples/microstrip/MSL.m b/matlab/examples/microstrip/MSL.m index 6f44529..b6ec0b3 100644 --- a/matlab/examples/microstrip/MSL.m +++ b/matlab/examples/microstrip/MSL.m @@ -180,8 +180,6 @@ xlabel( 'frequency f / MHz' ); ylabel( 'characteristic impedance Z / Ohm' ); legend( 'real', 'imag' ); - - %% visualize electric and magnetic fields % you will find vtk dump files in the simulation folder (tmp/) % use paraview to visualize them diff --git a/matlab/examples/microstrip/MSL2.m b/matlab/examples/microstrip/MSL2.m index fe4daaa..adbb817 100644 --- a/matlab/examples/microstrip/MSL2.m +++ b/matlab/examples/microstrip/MSL2.m @@ -1,139 +1,114 @@ % -% microstrip line example -% -% this example shows how to use a MSL port +% EXAMPLE / microstrip / MSL2 % +% This example shows how to use the MSL-port. % The MSL is excited at the center of the computational volume. The % boundary at xmin is an absorbing boundary (Mur) and at xmax an electric % wall. The reflection coefficient at this wall is S11 = -1. +% Direction of propagation is x. % - +% This example demonstrates: +% - simple microstrip geometry (made of PEC) +% - MSL port +% - MSL analysis +% +% +% Tested with +% - Matlab 2009b +% - Octave 3.3.52 +% - openEMS v0.0.14 +% +% (C) 2010 Sebastian Held close all clear clc -physical_constants - - -postprocessing_only = 0; - +%% switches +postproc_only = 0; %% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -drawingunits = 1e-6; % specify everything in um -MSL_length = 10000; -MSL_width = 1000; +physical_constants; +unit = 1e-6; % specify everything in um +MSL_length = 10000; +MSL_width = 1000; substrate_thickness = 254; +substrate_epr = 3.66; -mesh_res = [200 0 0]; -max_timesteps = 20000; -min_decrement = 1e-6; -f_max = 8e9; +% mesh_res = [200 0 0]; + +%% prepare simulation folder +Sim_Path = 'tmp'; +Sim_CSX = 'msl2.xml'; +if ~postproc_only + [status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory + [status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder +end %% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%% +max_timesteps = 20000; +min_decrement = 1e-6; +f_max = 7e9; FDTD = InitFDTD( max_timesteps, min_decrement, 'OverSampling', 10 ); FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2 ); -BC = [2 0 0 0 0 1]; +BC = {'MUR' 'PEC' 'PEC' 'PEC' 'PEC' 'PMC'}; FDTD = SetBoundaryCond( FDTD, BC ); %% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CSX = InitCSX(); -mesh.x = -MSL_length : mesh_res(1) : MSL_length; -mesh.y = linspace(-MSL_width/2,MSL_width/2,10); % discretize the width of the MSL with 10 cells -temp1 = linspace(-4*MSL_width,mesh.y(1),20); -temp2 = linspace(mesh.y(end),4*MSL_width,20); -mesh.y = [temp1(1:end-1), mesh.y, temp2(2:end)]; % add coarser discretization -mesh.z = linspace(0,substrate_thickness,5); % discretize the substrate with 5 cells -temp1 = linspace(substrate_thickness,2*substrate_thickness,5); -mesh.z = [mesh.z temp1(2:end)]; % add same space above the strip -temp1 = linspace(2*substrate_thickness,5*substrate_thickness,10); -mesh.z = [mesh.z temp1(2:end)]; % coarser discretization -CSX = DefineRectGrid( CSX, drawingunits, mesh ); - -%% Material definitions -CSX = AddMetal( CSX, 'PEC' ); -CSX = AddMaterial( CSX, 'RO4350B' ); +resolution = c0/(f_max*sqrt(substrate_epr))/unit /50; % resolution of lambda/50 +mesh.x = SmoothMeshLines( [-MSL_length MSL_length], resolution ); +mesh.y = SmoothMeshLines( [-4*MSL_width -MSL_width/2 MSL_width/2 4*MSL_width], resolution ); +mesh.z = SmoothMeshLines( [linspace(0,substrate_thickness,5) 10*substrate_thickness], resolution ); +CSX = DefineRectGrid( CSX, unit, mesh ); %% substrate -CSX = SetMaterialProperty( CSX, 'RO4350B', 'Epsilon', 3.66 ); +CSX = AddMaterial( CSX, 'RO4350B' ); +CSX = SetMaterialProperty( CSX, 'RO4350B', 'Epsilon', substrate_epr ); start = [mesh.x(1), mesh.y(1), 0]; stop = [mesh.x(end), mesh.y(end), substrate_thickness]; CSX = AddBox( CSX, 'RO4350B', 0, start, stop ); %% MSL port -CSX = AddExcitation( CSX, 'excite', 0, [0 0 1]); +CSX = AddMetal( CSX, 'PEC' ); portstart = [ 0, -MSL_width/2, substrate_thickness]; -portstop = [ MSL_length, MSL_width/2, 0]; -[CSX,portstruct] = AddMSLPort( CSX, 1, 'PEC', portstart, portstop, [1 0 0], [0 0 1], 'excite' ); +portstop = [ MSL_length, MSL_width/2, 0]; +[CSX,portstruct] = AddMSLPort( CSX, 1, 'PEC', portstart, portstop, [1 0 0], [0 0 1], [], 'excite' ); %% MSL start = [-MSL_length, -MSL_width/2, substrate_thickness]; stop = [ 0, MSL_width/2, substrate_thickness]; -CSX = AddBox( CSX, 'PEC', 0, start, stop ); +CSX = AddBox( CSX, 'PEC', 999, start, stop ); % priority needs to be higher than -%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -CSX = AddDump(CSX,'Et_','DumpType',0,'DumpMode',0); -start = [mesh.x(1) , mesh.y(1), substrate_thickness/2]; +%% define dump boxes +start = [mesh.x(1), mesh.y(1), substrate_thickness/2]; stop = [mesh.x(end), mesh.y(end), substrate_thickness/2]; -CSX = AddBox(CSX,'Et_',0 , start,stop); - -CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',0); -CSX = AddBox(CSX,'Ht_',0,start,stop); +CSX = AddDump( CSX, 'Et_', 'DumpType', 0,'DumpMode', 2 ); % cell interpolated +CSX = AddBox( CSX, 'Et_', 0, start, stop ); +CSX = AddDump( CSX, 'Ht_', 'DumpType', 1,'DumpMode', 2 ); % cell interpolated +CSX = AddBox( CSX, 'Ht_', 0, start, stop ); - -%% define openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -openEMS_opts = ''; -openEMS_opts = [openEMS_opts ' --disable-dumps']; -% openEMS_opts = [openEMS_opts ' --debug-material']; -% openEMS_opts = [openEMS_opts ' --debug-operator']; -% openEMS_opts = [openEMS_opts ' --debug-boxes']; -% openEMS_opts = [openEMS_opts ' --engine=sse-compressed']; -% openEMS_opts = [openEMS_opts ' --engine=multithreaded']; -openEMS_opts = [openEMS_opts ' --engine=fastest']; - -Sim_Path = 'tmp'; -Sim_CSX = 'MSL2.xml'; - -if ~postprocessing_only - rmdir(Sim_Path,'s'); -end -mkdir(Sim_Path); - -%% Write openEMS compatible xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% write openEMS compatible xml-file WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); -%% cd to working dir and run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -savePath = pwd; -cd(Sim_Path); %cd to working dir -args = [Sim_CSX ' ' openEMS_opts]; -if ~postprocessing_only - invoke_openEMS(args); -end -cd(savePath); +%% show the structure +CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); - - -%% postproc & do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -U = ReadUI({'port_ut1A','port_ut1B','et'},'tmp/'); -I = ReadUI({'port_it1A','port_it1B'},'tmp/'); -delta_t_2 = I.TD{1}.t(1) - U.TD{1}.t(1); % half time-step (s) - -% create finer frequency resolution -f = linspace( 0, f_max, 1601 ); -for n=1:numel(U.FD) - U.FD{n}.f = f; - U.FD{n}.val = DFT_time2freq( U.TD{n}.t, U.TD{n}.val, f ); -end -for n=1:numel(I.FD) - I.FD{n}.f = f; - I.FD{n}.val = DFT_time2freq( I.TD{n}.t, I.TD{n}.val, f ); - I.FD{n}.val = I.FD{n}.val .* exp(-1i*2*pi*I.FD{n}.f*delta_t_2); % compensate half time-step advance of H-field +%% run openEMS +openEMS_opts = ''; +openEMS_opts = [openEMS_opts ' --engine=fastest']; +% openEMS_opts = [openEMS_opts ' --debug-material']; +% openEMS_opts = [openEMS_opts ' --debug-boxes']; +% openEMS_opts = [openEMS_opts ' --debug-PEC']; +if ~postproc_only + RunOpenEMS( Sim_Path, Sim_CSX, openEMS_opts ); end -% interpolate et to the time spacing of the voltage probes -et = interp1( U.TD{3}.t, U.TD{3}.val, U.TD{1}.t ); -f = U.FD{1}.f; +%% postprocess +f = linspace( 1e6, f_max, 1601 ); +U = ReadUI( {'port_ut1A','port_ut1B','et'}, 'tmp/', f ); +I = ReadUI( {'port_it1A','port_it1B'}, 'tmp/', f ); % Z = (U.FD{1}.val+U.FD{2}.val)/2 ./ I.FD{1}.val; % plot( f*1e-9, [real(Z);imag(Z)],'Linewidth',2); @@ -143,23 +118,30 @@ f = U.FD{1}.f; % legend( {'real','imaginary'}, 'location', 'northwest' ) % title( 'line impedance (will fail in case of reflections!)' ); -% figure -% plotyy(U.TD{1}.t/1e-6,[U.TD{1}.val;U.TD{2}.val],U.TD{1}.t/1e-6,et); -% xlabel('time (us)'); -% ylabel('amplitude (V)'); -% grid on; -% title( 'Time domain voltage probes and excitation signal' ); -% -% figure -% plot(I.TD{1}.t/1e-6,[I.TD{1}.val;I.TD{2}.val]); -% xlabel('time (us)'); -% ylabel('amplitude (A)'); -% grid on; -% title( 'Time domain current probes' ); +figure +ax = plotyy( U.TD{1}.t/1e-6, [U.TD{1}.val;U.TD{2}.val], U.TD{3}.t/1e-6, U.TD{3}.val ); +xlabel( 'time (us)' ); +ylabel( 'amplitude (V)' ); +grid on +title( 'Time domain voltage probes and excitation signal' ); +legend( {'ut1A','ut1B','excitation'} ); +% now make the y-axis symmetric to y=0 (align zeros of y1 and y2) +y1 = ylim(ax(1)); +y2 = ylim(ax(2)); +ylim( ax(1), [-max(abs(y1)) max(abs(y1))] ); +ylim( ax(2), [-max(abs(y2)) max(abs(y2))] ); +figure +plot( I.TD{1}.t/1e-6, [I.TD{1}.val;I.TD{2}.val] ); +xlabel( 'time (us)' ); +ylabel( 'amplitude (A)' ); +grid on +title( 'Time domain current probes' ); +legend( {'it1A','it1B'} ); -%% port analysis +% port analysis [S11,beta,ZL] = calcMSLPort( portstruct, Sim_Path, f ); +% attention! the reflection coefficient S11 is normalized to ZL! figure plot( sin(0:0.01:2*pi), cos(0:0.01:2*pi), 'Color', [.7 .7 .7] ); @@ -180,8 +162,9 @@ title( 'Reflection coefficient S11 at the measurement plane' ); figure plotyy( f/1e9, 20*log10(abs(S11)), f/1e9, angle(S11)/pi*180 ); -legend( {'abs(S11)', 'angle(S11)'} ); +legend( {'|S11|', 'angle(S11)'} ); xlabel( 'frequency (GHz)' ); +ylabel( '|S11| (dB)' ); title( 'Reflection coefficient S11 at the measurement plane' ); figure @@ -198,7 +181,6 @@ ylabel('impedance (Ohm)'); grid on; legend( {'real','imaginary'}, 'location', 'northeast' ) title( 'Characteristic line impedance ZL' ); -ylim( [-2*mean(real(ZL)) 2*mean(real(ZL))] ); % reference plane shift (to the end of the port) ref_shift = abs(portstop(1) - portstart(1)); @@ -219,3 +201,7 @@ plot( S11, 'k' ); plot( real(S11(1)), imag(S11(1)), '*r' ); axis equal title( 'Reflection coefficient S11 at the reference plane (at the electric wall)' ); + +%% visualize electric and magnetic fields +% you will find vtk dump files in the simulation folder (tmp/) +% use paraview to visualize them diff --git a/paraview/view_debug-PEC.cpd b/paraview/view_debug-PEC.cpd new file mode 100644 index 0000000..a521f66 --- /dev/null +++ b/paraview/view_debug-PEC.cpd @@ -0,0 +1,617 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +