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

matlab examples: cleaned up MSL2; updated MSL port

Browse Source
This commit is contained in:
Sebastian Held 2010-10-19 12:33:03 +02:00
parent da3e55a65c
commit 244b92e529
4 changed files with 100 additions and 120 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
matlab/AddMSLPort.m
View File

@ -1,5 +1,5 @@
function [CSX,port] = AddMSLPort( CSX, portnr, materialname, start, stop, dir, evec, refplaneshift, excitename ) function [CSX,port] = AddMSLPort( CSX, portnr, materialname, start, stop, dir, evec, refplaneshift, excitename )
% [CSX,port] = AddMSLPort( CSX, portnr, materialname, start, stop, dir, evec, excitename ) % [CSX,port] = AddMSLPort( CSX, portnr, materialname, start, stop, dir, evec, refplaneshift, excitename )
% %
% CSX: CSX-object created by InitCSX() % CSX: CSX-object created by InitCSX()
% portnr: (integer) number of the port % portnr: (integer) number of the port
@ -152,6 +152,7 @@ port.direction = direction;
% port.idx_height = idx_height; % port.idx_height = idx_height;
port.excite = 0; port.excite = 0;
port.refplaneshift = 0; port.refplaneshift = 0;
port.measplanepos = abs(v2_start(idx_prop) - start(idx_prop));
if (nargin >= 8) && (~isempty(refplaneshift)) if (nargin >= 8) && (~isempty(refplaneshift))
% refplaneshift counts from start of port % refplaneshift counts from start of port

13
matlab/calcMSLPort.m
View File

@ -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) % ref_shift: (optional) reference plane shift measured from start of port (in drawing units)
% %
% output: % output:
% S11: reflection coefficient % S11: reflection coefficient (normalized to ZL)
% beta: propagation constant % beta: propagation constant
% ZL: characteristic line impedance % ZL: characteristic line impedance
% %
@ -21,7 +21,7 @@ function [S11,beta,ZL] = calcMSLPort( portstruct, SimDir, f, ref_shift )
% %
% openEMS matlab interface % openEMS matlab interface
% ----------------------- % -----------------------
% Sebastian Held <sebastian.held@uni-due.de> % (C) 2010 Sebastian Held <sebastian.held@uni-due.de>
% See also AddMSLPort % See also AddMSLPort
% check % check
@ -90,15 +90,10 @@ S11 = (-1i * dEt + Et .* temp) ./ (Et .* temp + 1i * dEt); % solution 1
% determine ZL % determine ZL
ZL = sqrt(Et .* dEt ./ (Ht .* dHt)); ZL = sqrt(Et .* dEt ./ (Ht .* dHt));
% reference plane shift % reference plane shift (lossless)
if (nargin > 3) if (nargin > 3)
% renormalize the shift to the measurement plane % renormalize the shift to the measurement plane
if (portstruct.stop(portstruct.idx_prop) - portstruct.start(portstruct.idx_prop) > 0) ref_shift = ref_shift - portstruct.measplanepos;
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.drawingunit; ref_shift = ref_shift * portstruct.drawingunit;
S11 = S11 .* exp(2i*real(beta)*ref_shift); S11 = S11 .* exp(2i*real(beta)*ref_shift);
S11_corrected = S11_corrected .* exp(2i*real(beta)*ref_shift); S11_corrected = S11_corrected .* exp(2i*real(beta)*ref_shift);

2
matlab/examples/microstrip/MSL.m
View File

@ -180,8 +180,6 @@ xlabel( 'frequency f / MHz' );
ylabel( 'characteristic impedance Z / Ohm' ); ylabel( 'characteristic impedance Z / Ohm' );
legend( 'real', 'imag' ); legend( 'real', 'imag' );
%% visualize electric and magnetic fields %% visualize electric and magnetic fields
% you will find vtk dump files in the simulation folder (tmp/) % you will find vtk dump files in the simulation folder (tmp/)
% use paraview to visualize them % use paraview to visualize them

202
matlab/examples/microstrip/MSL2.m
View File

@ -1,139 +1,114 @@
% %
% microstrip line example % EXAMPLE / microstrip / MSL2
%
% this example shows how to use a MSL port
% %
% This example shows how to use the MSL-port.
% The MSL is excited at the center of the computational volume. The % 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 % boundary at xmin is an absorbing boundary (Mur) and at xmax an electric
% wall. The reflection coefficient at this wall is S11 = -1. % 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 <sebastian.held@uni-due.de>
close all close all
clear clear
clc clc
physical_constants %% switches
postproc_only = 0;
postprocessing_only = 0;
%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
drawingunits = 1e-6; % specify everything in um physical_constants;
MSL_length = 10000; unit = 1e-6; % specify everything in um
MSL_width = 1000; MSL_length = 10000;
MSL_width = 1000;
substrate_thickness = 254; substrate_thickness = 254;
substrate_epr = 3.66;
mesh_res = [200 0 0]; % mesh_res = [200 0 0];
max_timesteps = 20000;
min_decrement = 1e-6; %% prepare simulation folder
f_max = 8e9; 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%% %% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%
max_timesteps = 20000;
min_decrement = 1e-6;
f_max = 7e9;
FDTD = InitFDTD( max_timesteps, min_decrement, 'OverSampling', 10 ); FDTD = InitFDTD( max_timesteps, min_decrement, 'OverSampling', 10 );
FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2 ); 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 ); FDTD = SetBoundaryCond( FDTD, BC );
%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CSX = InitCSX(); CSX = InitCSX();
mesh.x = -MSL_length : mesh_res(1) : MSL_length; resolution = c0/(f_max*sqrt(substrate_epr))/unit /50; % resolution of lambda/50
mesh.y = linspace(-MSL_width/2,MSL_width/2,10); % discretize the width of the MSL with 10 cells mesh.x = SmoothMeshLines( [-MSL_length MSL_length], resolution );
temp1 = linspace(-4*MSL_width,mesh.y(1),20); mesh.y = SmoothMeshLines( [-4*MSL_width -MSL_width/2 MSL_width/2 4*MSL_width], resolution );
temp2 = linspace(mesh.y(end),4*MSL_width,20); mesh.z = SmoothMeshLines( [linspace(0,substrate_thickness,5) 10*substrate_thickness], resolution );
mesh.y = [temp1(1:end-1), mesh.y, temp2(2:end)]; % add coarser discretization CSX = DefineRectGrid( CSX, unit, mesh );
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' );
%% substrate %% 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]; start = [mesh.x(1), mesh.y(1), 0];
stop = [mesh.x(end), mesh.y(end), substrate_thickness]; stop = [mesh.x(end), mesh.y(end), substrate_thickness];
CSX = AddBox( CSX, 'RO4350B', 0, start, stop ); CSX = AddBox( CSX, 'RO4350B', 0, start, stop );
%% MSL port %% MSL port
CSX = AddExcitation( CSX, 'excite', 0, [0 0 1]); CSX = AddMetal( CSX, 'PEC' );
portstart = [ 0, -MSL_width/2, substrate_thickness]; portstart = [ 0, -MSL_width/2, substrate_thickness];
portstop = [ MSL_length, MSL_width/2, 0]; portstop = [ MSL_length, MSL_width/2, 0];
[CSX,portstruct] = AddMSLPort( CSX, 1, 'PEC', portstart, portstop, [1 0 0], [0 0 1], 'excite' ); [CSX,portstruct] = AddMSLPort( CSX, 1, 'PEC', portstart, portstop, [1 0 0], [0 0 1], [], 'excite' );
%% MSL %% MSL
start = [-MSL_length, -MSL_width/2, substrate_thickness]; start = [-MSL_length, -MSL_width/2, substrate_thickness];
stop = [ 0, 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... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% define dump boxes
CSX = AddDump(CSX,'Et_','DumpType',0,'DumpMode',0); start = [mesh.x(1), mesh.y(1), substrate_thickness/2];
start = [mesh.x(1) , mesh.y(1), substrate_thickness/2];
stop = [mesh.x(end), mesh.y(end), substrate_thickness/2]; stop = [mesh.x(end), mesh.y(end), substrate_thickness/2];
CSX = AddBox(CSX,'Et_',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',0); CSX = AddDump( CSX, 'Ht_', 'DumpType', 1,'DumpMode', 2 ); % cell interpolated
CSX = AddBox(CSX,'Ht_',0,start,stop); CSX = AddBox( CSX, 'Ht_', 0, start, stop );
%% write openEMS compatible xml-file
%% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
%% cd to working dir and run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% show the structure
savePath = pwd; CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
cd(Sim_Path); %cd to working dir
args = [Sim_CSX ' ' openEMS_opts];
if ~postprocessing_only
invoke_openEMS(args);
end
cd(savePath);
%% run openEMS
openEMS_opts = '';
%% postproc & do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% openEMS_opts = [openEMS_opts ' --engine=fastest'];
U = ReadUI({'port_ut1A','port_ut1B','et'},'tmp/'); % openEMS_opts = [openEMS_opts ' --debug-material'];
I = ReadUI({'port_it1A','port_it1B'},'tmp/'); % openEMS_opts = [openEMS_opts ' --debug-boxes'];
delta_t_2 = I.TD{1}.t(1) - U.TD{1}.t(1); % half time-step (s) % openEMS_opts = [openEMS_opts ' --debug-PEC'];
if ~postproc_only
% create finer frequency resolution RunOpenEMS( Sim_Path, Sim_CSX, openEMS_opts );
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
end 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; % Z = (U.FD{1}.val+U.FD{2}.val)/2 ./ I.FD{1}.val;
% plot( f*1e-9, [real(Z);imag(Z)],'Linewidth',2); % plot( f*1e-9, [real(Z);imag(Z)],'Linewidth',2);
@ -143,23 +118,30 @@ f = U.FD{1}.f;
% legend( {'real','imaginary'}, 'location', 'northwest' ) % legend( {'real','imaginary'}, 'location', 'northwest' )
% title( 'line impedance (will fail in case of reflections!)' ); % title( 'line impedance (will fail in case of reflections!)' );
% figure figure
% plotyy(U.TD{1}.t/1e-6,[U.TD{1}.val;U.TD{2}.val],U.TD{1}.t/1e-6,et); 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)'); xlabel( 'time (us)' );
% ylabel('amplitude (V)'); ylabel( 'amplitude (V)' );
% grid on; grid on
% title( 'Time domain voltage probes and excitation signal' ); title( 'Time domain voltage probes and excitation signal' );
% legend( {'ut1A','ut1B','excitation'} );
% figure % now make the y-axis symmetric to y=0 (align zeros of y1 and y2)
% plot(I.TD{1}.t/1e-6,[I.TD{1}.val;I.TD{2}.val]); y1 = ylim(ax(1));
% xlabel('time (us)'); y2 = ylim(ax(2));
% ylabel('amplitude (A)'); ylim( ax(1), [-max(abs(y1)) max(abs(y1))] );
% grid on; ylim( ax(2), [-max(abs(y2)) max(abs(y2))] );
% title( 'Time domain current probes' );
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 ); [S11,beta,ZL] = calcMSLPort( portstruct, Sim_Path, f );
% attention! the reflection coefficient S11 is normalized to ZL!
figure figure
plot( sin(0:0.01:2*pi), cos(0:0.01:2*pi), 'Color', [.7 .7 .7] ); 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 figure
plotyy( f/1e9, 20*log10(abs(S11)), f/1e9, angle(S11)/pi*180 ); 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)' ); xlabel( 'frequency (GHz)' );
ylabel( '|S11| (dB)' );
title( 'Reflection coefficient S11 at the measurement plane' ); title( 'Reflection coefficient S11 at the measurement plane' );
figure figure
@ -198,7 +181,6 @@ ylabel('impedance (Ohm)');
grid on; grid on;
legend( {'real','imaginary'}, 'location', 'northeast' ) legend( {'real','imaginary'}, 'location', 'northeast' )
title( 'Characteristic line impedance ZL' ); title( 'Characteristic line impedance ZL' );
ylim( [-2*mean(real(ZL)) 2*mean(real(ZL))] );
% reference plane shift (to the end of the port) % reference plane shift (to the end of the port)
ref_shift = abs(portstop(1) - portstart(1)); ref_shift = abs(portstop(1) - portstart(1));
@ -219,3 +201,7 @@ plot( S11, 'k' );
plot( real(S11(1)), imag(S11(1)), '*r' ); plot( real(S11(1)), imag(S11(1)), '*r' );
axis equal axis equal
title( 'Reflection coefficient S11 at the reference plane (at the electric wall)' ); 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