caiyu/openEMS
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

TFSF: some critical fixes & Tutorial updates

Browse Source
This commit is contained in:
Thorsten Liebig 2012-07-23 12:10:35 +02:00
parent df083a63cc
commit 3546fcc97d
5 changed files with 116 additions and 65 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

90
FDTD/extensions/operator_ext_tfsf.cpp
View File

@ -43,6 +43,8 @@ void Operator_Ext_TFST::Init()
m_CurrAmp[n][l][c]=NULL;
}
m_Frequency = 0.0;
m_PhVel = __C0__;
Operator_Extension::Init();
}
@ -95,6 +97,7 @@ bool Operator_Ext_TFST::BuildExtension()
return false;
}
m_PhVel = __C0__;
CSPropExcitation* elec=NULL;
CSProperties* prop=NULL;
CSPrimitives* prim=NULL;
@ -144,6 +147,20 @@ bool Operator_Ext_TFST::BuildExtension()
return false;
}
m_Frequency = elec->GetFrequency();
// if (m_Frequency<=0)
// m_Frequency = m_Op->GetExcitationSignal()->GetFrequencyOfInterest();
if (m_Frequency<=0)
m_PhVel=__C0__;
else
m_PhVel=m_Op->CalcNumericPhaseVelocity(m_Start,m_Stop,m_PropDir,m_Frequency);
if ((m_PhVel<0) || (m_PhVel>__C0__) || isnan(m_PhVel))
{
cerr << "Operator_Ext_TFST::BuildExtension: Warning, invalid phase velocity found, resetting to c0! " << endl;
m_PhVel = __C0__;
}
double origin[3];
unsigned int ui_origin[3];
for (int n=0;n<3;++n)
@ -163,11 +180,11 @@ bool Operator_Ext_TFST::BuildExtension()
if (m_IncLow[n])
{
ui_origin[n] = m_Start[n];
ui_origin[n] = m_Start[n]-1;
}
else
{
ui_origin[n] = m_Stop[n];
ui_origin[n] = m_Stop[n]+1;
}
origin[n] = m_Op->GetDiscLine(n,ui_origin[n]);
}
@ -244,16 +261,16 @@ bool Operator_Ext_TFST::BuildExtension()
if (m_ActiveDir[n][0])
{
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;
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/m_PhVel/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;
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/m_PhVel/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);
@ -268,16 +285,16 @@ bool Operator_Ext_TFST::BuildExtension()
{
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;
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/m_PhVel/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;
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/m_PhVel/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);
@ -287,20 +304,10 @@ bool Operator_Ext_TFST::BuildExtension()
m_CurrAmp[n][1][1][ui_pos]*=m_Op->GetIV(nPP,pos);
}
if (m_IncLow[n])
{
if (m_ActiveDir[n][0])
m_CurrAmp[n][0][0][ui_pos]*=-1;
if (m_ActiveDir[n][1])
m_CurrAmp[n][1][1][ui_pos]*=-1;
}
else
{
if (m_ActiveDir[n][0])
m_CurrAmp[n][0][1][ui_pos]*=-1;
if (m_ActiveDir[n][1])
m_CurrAmp[n][1][0][ui_pos]*=-1;
}
if (m_ActiveDir[n][0])
m_CurrAmp[n][0][0][ui_pos]*=-1;
if (m_ActiveDir[n][1])
m_CurrAmp[n][1][1][ui_pos]*=-1;
if (pos[nP]==m_Stop[nP])
{
@ -322,16 +329,16 @@ bool Operator_Ext_TFST::BuildExtension()
if (m_ActiveDir[n][0])
{
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;
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/m_PhVel/dT + 1.0;
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;
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/m_PhVel/dT + 1.0;
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);
@ -346,16 +353,16 @@ bool Operator_Ext_TFST::BuildExtension()
if (m_ActiveDir[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;
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/m_PhVel/dT + 1.0;
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;
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/m_PhVel/dT + 1.0;
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);
@ -365,20 +372,10 @@ bool Operator_Ext_TFST::BuildExtension()
m_VoltAmp[n][1][1][ui_pos]*=m_Op->GetVI(nPP,pos);
}
if (m_IncLow[n])
{
if (m_ActiveDir[n][1])
m_VoltAmp[n][1][0][ui_pos]*=-1;
if (m_ActiveDir[n][0])
m_VoltAmp[n][0][1][ui_pos]*=-1;
}
else
{
if (m_ActiveDir[n][0])
m_VoltAmp[n][0][0][ui_pos]*=-1;
if (m_ActiveDir[n][1])
m_VoltAmp[n][1][1][ui_pos]*=-1;
}
if (m_ActiveDir[n][1])
m_VoltAmp[n][1][0][ui_pos]*=-1;
if (m_ActiveDir[n][0])
m_VoltAmp[n][0][1][ui_pos]*=-1;
if (pos[nP]==m_Stop[nP])
{
@ -418,6 +415,7 @@ void Operator_Ext_TFST::ShowStat(ostream &ostr) const
Operator_Extension::ShowStat(ostr);
cout << "Active directions\t: " << "(" << m_ActiveDir[0][0] << "/" << m_ActiveDir[0][1] << ", " << m_ActiveDir[1][0] << "/" << m_ActiveDir[1][1] << ", " << m_ActiveDir[2][0] << "/" << m_ActiveDir[2][1] << ")" << endl;
cout << "Propagation direction\t: " << "(" << m_PropDir[0] << ", " << m_PropDir[1] << ", " << m_PropDir[2] << ")" << endl;
cout << "Rel. propagation speed\t: " << m_PhVel/__C0__ << "*c0 @ " << m_Frequency << " Hz" << endl;
cout << "E-field amplitude\t: " << "(" << m_E_Amp[0] << ", " << m_E_Amp[1] << ", " << m_E_Amp[2] << ")" << endl;
cout << "H-field 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;

3
FDTD/extensions/operator_ext_tfsf.h
View File

@ -60,6 +60,9 @@ protected:
double m_E_Amp[3];
double m_H_Amp[3];
double m_Frequency;
double m_PhVel;
unsigned int m_maxDelay;
// array setup [direction][low/high][component][ <mesh_position> ]

27
FDTD/operator.cpp
View File

@ -1740,3 +1740,30 @@ void Operator::DeleteExtension(Operator_Extension* op_ext)
}
}
}
double Operator::CalcNumericPhaseVelocity(unsigned int start[3], unsigned int stop[3], double propDir[3], float freq) const
{
double phv=__C0__;
double average_mesh_disc[3];
// double k=2*PI*freq/__C0__;
for (int n=0;n<3;++n)
{
average_mesh_disc[n] = fabs(GetDiscLine(n,start[n])-GetDiscLine(n,stop[n]))*GetGridDelta() / (fabs(stop[n]-start[n]));
}
for (int n=0;n<3;++n)
{
int nP = (n+1)%3;
int nPP = (n+2)%3;
if ((fabs(propDir[n])==1) && (propDir[nP]==0) && (propDir[nPP]==0))
{
double kx = asin(average_mesh_disc[0]/__C0__/dT*sin(2*PI*freq*dT/2))*2/average_mesh_disc[0];
return 2*PI*freq/kx;
}
}
cerr << "Operator::CalcNumericPhaseVelocity: Warning, propagation direction not in Cartesian direction, assuming phase velocity to be c0" << endl;
return phv;
}

2
FDTD/operator.h
View File

@ -147,6 +147,8 @@ public:
Operator_Ext_Excitation* GetExcitationExtension() const;
virtual double CalcNumericPhaseVelocity(unsigned int start[3], unsigned int stop[3], double propDir[3], float freq) const;
protected:
//! use New() for creating a new Operator
Operator();

59
matlab/Tutorials/RCS_Sphere.m
View File

@ -20,13 +20,16 @@ unit = 1e-3; % all length in mm
sphere.rad = 200;
inc_angle = 0 /180*pi; %incident angle (to x-axis) in rad
% size of the simulation box
SimBox = 1000;
PW_Box = 750;
%% setup FDTD parameter & excitation function
f_start = 200e6; % start frequency
f_start = 50e6; % start frequency
f_stop = 1000e6; % stop frequency
f0 = 500e6;
FDTD = InitFDTD( 30000 );
FDTD = SetGaussExcite( FDTD, 0.5*(f_start+f_stop), 0.5*(f_stop-f_start) );
@ -34,12 +37,12 @@ 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);
%create mesh
smooth_mesh = SmoothMeshLines([0 SimBox/2], max_res);
mesh.x = unique([-smooth_mesh smooth_mesh]);
mesh.y = mesh.x;
mesh.z = mesh.x;
@ -48,7 +51,7 @@ 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
k_dir = [cos(inc_angle) sin(inc_angle) 0]; % plane wave direction
E_dir = [0 0 1]; % plane wave polarization --> E_z
CSX = AddPlaneWaveExcite(CSX, 'plane_wave', k_dir, E_dir);
@ -92,26 +95,44 @@ RunOpenEMS( Sim_Path, Sim_CSX);
disp('Use Paraview to display the elctric fields dumped by openEMS');
%%
freq = 500e6;
EF = ReadUI( 'et', Sim_Path, freq ); % time domain/freq domain voltage
EF = ReadUI( 'et', Sim_Path, f0 ); % time domain/freq domain voltage
Pin = 0.5*norm(E_dir)^2/Z0 .* abs(EF.FD{1}.val).^2;
%%
nf2ff = CalcNF2FF(nf2ff, Sim_Path, freq, [-180:2:180]*pi/180, [0 90]*pi/180,'Mode',1);
RCS = 4*pi./Pin(1).*nf2ff.P_rad{1}(:,1);
polar(nf2ff.theta,RCS);
xlabel('z -->');
ylabel('x -->');
nf2ff = CalcNF2FF(nf2ff, Sim_Path, f0, pi/2, [-180:2:180]*pi/180, 'Mode',1);
RCS = 4*pi./Pin(1).*nf2ff.P_rad{1}(:);
polar(nf2ff.phi,RCS);
xlabel('x -->');
ylabel('y -->');
hold on
grid on
drawnow
%%
disp( 'calculating 3D far field pattern and dumping to vtk (use Paraview to visualize)...' );
thetaRange = (0:2:180);
phiRange = (0:2:360) - 180;
nf2ff = CalcNF2FF(nf2ff, Sim_Path, freq, thetaRange*pi/180, phiRange*pi/180,'Verbose',1,'Outfile','3D_Pattern.h5','Mode',1);
freq = linspace(f_start,f_stop,100);
EF = ReadUI( 'et', Sim_Path, freq ); % time domain/freq domain voltage
Pin = 0.5*norm(E_dir)^2/Z0 .* abs(EF.FD{1}.val).^2;
RCS = 4*pi./Pin(1).*nf2ff.P_rad{1};
DumpFF2VTK([Sim_Path '/3D_Pattern.vtk'],RCS,thetaRange,phiRange,1e-3);
nf2ff = CalcNF2FF(nf2ff, Sim_Path, freq, pi/2, pi+inc_angle, 'Mode',1);
for fn=1:numel(freq)
back_scat(fn) = 4*pi./Pin(fn).*nf2ff.P_rad{fn}(1);
end
disp('Use Paraview to display the 3D scattering pattern');
%%
figure
plot(freq/1e6,back_scat,'Linewidth',2);
grid on;
xlabel('frequency (MHz) \rightarrow');
ylabel('RCS (m^2) \rightarrow');
title('radar cross section');
%%
figure
lambda = c0./freq;
semilogy(sphere.rad*unit./lambda,back_scat/(pi*sphere.rad*unit*sphere.rad*unit),'Linewidth',2);
ylim([10^-2 10^1])
grid on;
xlabel('sphere radius / wavelength \rightarrow');
ylabel('RCS / (\pi a^2) \rightarrow');
title('normalized radar cross section');