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

waveguide examples update

Browse Source
This commit is contained in:
Thorsten Liebig 2011-01-02 11:32:00 +01:00
parent 7d09021dd4
commit dce988b52f
3 changed files with 294 additions and 158 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

192
matlab/examples/waveguide/Coax.m
View File

@ -1,43 +1,59 @@
%
% EXAMPLE / waveguide / coaxial cable
%
% This example demonstrates how to:
% - setup a coaxial waveguide
% - use analytic functions for waveguide excitations and voltage/current
% calculations
%
%
% Tested with
% - Matlab 2009b
% - openEMS v0.0.17
%
% (C) 2010 Thorsten Liebig <thorsten.liebig@uni-due.de>
close all close all
clear clear
clc clc
%% switches & options...
postprocessing_only = 0;
use_pml = 0; % use pml boundaries instead of mur
openEMS_opts = '';
% openEMS_opts = [openEMS_opts ' --disable-dumps'];
%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
abs_length = 250; numTS = 5000; %number of timesteps
length = 1000; length = 1000; %length of the waveguide
coax_rad_i = 100; unit = 1e-3; %drawing unit used
coax_rad_ai = 230; coax_rad_i = 100; %inner radius
coax_rad_aa = 240; coax_rad_ai = 230; %inner radius of outer cladding
mesh_res = [5 5 5]; coax_rad_aa = 240; %outer radius of outer cladding
mesh_res = [5 5 5]; %desired mesh resolution
EPS0 = 8.85418781762e-12; physical_constants;
MUE0 = 1.256637062e-6;
C0 = 1/sqrt(EPS0*MUE0);
Z0 = sqrt(MUE0/EPS0);
%excitation
f0 = 0.5e9; f0 = 0.5e9;
epsR = 1; epsR = 1;
%% define openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% create sim path %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
openEMS_opts = '';
% openEMS_opts = [openEMS_opts ' --debug-material'];
% openEMS_opts = [openEMS_opts ' --debug-operator'];
% openEMS_opts = [openEMS_opts ' --disable-dumps --engine=fastest'];
openEMS_opts = [openEMS_opts ' --engine=fastest'];
Sim_Path = 'tmp'; Sim_Path = 'tmp';
Sim_CSX = 'coax.xml'; Sim_CSX = 'coax.xml';
if (exist(Sim_Path,'dir')) if (postprocessing_only==0)
rmdir(Sim_Path,'s'); [status, message, messageid] = rmdir(Sim_Path,'s');
[status, message, messageid] = mkdir(Sim_Path);
end end
mkdir(Sim_Path);
%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FDTD = InitFDTD(5e5,1e-5); FDTD = InitFDTD(numTS,1e-5);
FDTD = SetGaussExcite(FDTD,f0,f0); FDTD = SetGaussExcite(FDTD,f0,f0);
BC = [0 0 0 0 0 0]; %electric walls only BC = {'PEC','PEC','PEC','PEC','PEC','MUR'};
if (use_pml>0)
BC = {'PEC','PEC','PEC','PEC','PEC','PML_8'};
end
FDTD = SetBoundaryCond(FDTD,BC); FDTD = SetBoundaryCond(FDTD,BC);
%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -45,16 +61,7 @@ CSX = InitCSX();
mesh.x = -2.5*mesh_res(1)-coax_rad_aa : mesh_res(1) : coax_rad_aa+2.5*mesh_res(1); mesh.x = -2.5*mesh_res(1)-coax_rad_aa : mesh_res(1) : coax_rad_aa+2.5*mesh_res(1);
mesh.y = mesh.x; mesh.y = mesh.x;
mesh.z = 0 : mesh_res(3) : length; mesh.z = 0 : mesh_res(3) : length;
CSX = DefineRectGrid(CSX, 1e-3,mesh); CSX = DefineRectGrid(CSX, unit, mesh);
%% fake pml %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
finalKappa = 0.3/abs_length^4;
finalSigma = finalKappa*MUE0/EPS0;
CSX = AddMaterial(CSX,'pml');
CSX = SetMaterialProperty(CSX,'pml','Kappa',finalKappa);
CSX = SetMaterialProperty(CSX,'pml','Sigma',finalSigma);
CSX = SetMaterialWeight(CSX,'pml','Kappa',['pow(abs(z)-' num2str(length-abs_length) ',4)']);
CSX = SetMaterialWeight(CSX,'pml','Sigma',['pow(abs(z)-' num2str(length-abs_length) ',4)']);
%%% coax %%% coax
CSX = AddMaterial(CSX,'copper'); CSX = AddMaterial(CSX,'copper');
@ -62,11 +69,10 @@ CSX = SetMaterialProperty(CSX,'copper','Kappa',56e6);
start = [0, 0 , 0];stop = [0, 0 , length]; start = [0, 0 , 0];stop = [0, 0 , length];
CSX = AddCylinder(CSX,'copper',0 ,start,stop,coax_rad_i); CSX = AddCylinder(CSX,'copper',0 ,start,stop,coax_rad_i);
CSX = AddCylindricalShell(CSX,'copper',0 ,start,stop,0.5*(coax_rad_aa+coax_rad_ai),(coax_rad_aa-coax_rad_ai)); CSX = AddCylindricalShell(CSX,'copper',0 ,start,stop,0.5*(coax_rad_aa+coax_rad_ai),(coax_rad_aa-coax_rad_ai));
start(3) = length-abs_length;
CSX = AddCylindricalShell(CSX,'pml',0 ,start,stop,0.5*(coax_rad_i+coax_rad_ai),(coax_rad_ai-coax_rad_i));
%% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
start(3) = 0; stop(3)=mesh_res(1)/2; start = [0,0,-0.1];
stop = [0,0, 0.1];
CSX = AddExcitation(CSX,'excite',0,[1 1 0]); CSX = AddExcitation(CSX,'excite',0,[1 1 0]);
weight{1} = '(x)/(x*x+y*y)'; weight{1} = '(x)/(x*x+y*y)';
weight{2} = 'y/pow(rho,2)'; weight{2} = 'y/pow(rho,2)';
@ -83,37 +89,101 @@ CSX = AddBox(CSX,'Et_',0 , start,stop);
CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',2); CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',2);
CSX = AddBox(CSX,'Ht_',0,start,stop); CSX = AddBox(CSX,'Ht_',0,start,stop);
%% define voltage calc boxes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%voltage calc %voltage calc
CSX = AddProbe(CSX,'ut1_1',0); CSX = AddProbe(CSX,'ut1',0);
start = [ coax_rad_i 0 length/2 ];stop = [ coax_rad_ai 0 length/2 ]; start = [ coax_rad_i 0 mesh.z(10) ];
CSX = AddBox(CSX,'ut1_1', 0 ,start,stop); stop = [ coax_rad_ai 0 mesh.z(10) ];
CSX = AddProbe(CSX,'ut1_2',0); CSX = AddBox(CSX,'ut1', 0 ,start,stop);
start = [ coax_rad_i 0 length/2+mesh_res(3) ];stop = [ coax_rad_ai 0 length/2+mesh_res(3) ]; CSX = AddProbe(CSX,'ut2',0);
CSX = AddBox(CSX,'ut1_2', 0 ,start,stop); start = [ coax_rad_i 0 mesh.z(end-10)];
stop = [ coax_rad_ai 0 mesh.z(end-10)];
CSX = AddBox(CSX,'ut2', 0 ,start,stop);
%current calc %current calc, for each position there are two currents, which will get
CSX = AddProbe(CSX,'it1',1); %averaged to match the voltage position in between (!Yee grid!)
CSX = AddProbe(CSX,'it1a',1);
mid = 0.5*(coax_rad_i+coax_rad_ai); mid = 0.5*(coax_rad_i+coax_rad_ai);
start = [ -mid -mid length/2 ];stop = [ mid mid length/2 ]; start = [ -mid -mid mesh.z(9) ];
CSX = AddBox(CSX,'it1', 0 ,start,stop); stop = [ mid mid mesh.z(9) ];
CSX = AddBox(CSX,'it1a', 0 ,start,stop);
CSX = AddProbe(CSX,'it1b',1);
start = [ -mid -mid mesh.z(10) ];
stop = [ mid mid mesh.z(10) ];
CSX = AddBox(CSX,'it1b', 0 ,start,stop);
%Write openEMS compatoble xml-file CSX = AddProbe(CSX,'it2a',1);
WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); start = [ -mid -mid mesh.z(end-11) ];
stop = [ mid mid mesh.z(end-11) ];
CSX = AddBox(CSX,'it2a', 0 ,start,stop);
CSX = AddProbe(CSX,'it2b',1);
start = [ -mid -mid mesh.z(end-10) ];
stop = [ mid mid mesh.z(end-10) ];
CSX = AddBox(CSX,'it2b', 0 ,start,stop);
%% cd to working dir and run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Write openEMS
RunOpenEMS(Sim_Path,Sim_CSX,openEMS_opts); if (postprocessing_only==0)
WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts);
end
%% postproc & do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%
UI = ReadUI({'ut1_1','ut1_2','it1'},Sim_Path); freq = linspace(0,2*f0,201);
u_f = (UI.FD{1}.val + UI.FD{2}.val)/2; %averaging voltages to fit current U = ReadUI({'ut1','ut2'},[Sim_Path '/'],freq);
i_f = UI.FD{3}.val; I = ReadUI({'it1a','it1b','it2a','it2b'},[Sim_Path '/'],freq);
Exc = ReadUI('et',Sim_Path,freq);
ZL = Z0/2/pi/sqrt(epsR)*log(coax_rad_ai/coax_rad_i); %analytic line-impedance of a coax %% plot voltages
plot(UI.FD{1}.f,ZL*ones(size(u_f)),'g'); figure
plot(U.TD{1}.t, U.TD{1}.val,'Linewidth',2);
hold on; hold on;
grid on; grid on;
Z = u_f./i_f; plot(U.TD{2}.t, U.TD{2}.val,'r--','Linewidth',2);
plot(UI.FD{1}.f,real(Z),'Linewidth',2); xlabel('time (s)')
plot(UI.FD{1}.f,imag(Z),'r','Linewidth',2); ylabel('voltage (V)')
xlim([0 2*f0]); legend('u_1(t)','u_2(t)')
legend('Z_L','\Re\{Z\}','\Im\{Z\}','Location','Best');
%% calculate incoming and reflected voltages & currents
ZL_a = ones(size(freq))*Z0/2/pi/sqrt(epsR)*log(coax_rad_ai/coax_rad_i); %analytic line-impedance of a coax
uf1 = U.FD{1}.val./Exc.FD{1}.val;
uf2 = U.FD{2}.val./Exc.FD{1}.val;
if1 = 0.5*(I.FD{1}.val+I.FD{2}.val)./Exc.FD{1}.val;
if2 = 0.5*(I.FD{3}.val+I.FD{4}.val)./Exc.FD{1}.val;
uf1_inc = 0.5 * ( uf1 + if1 .* ZL_a );
if1_inc = 0.5 * ( if1 + uf1 ./ ZL_a );
uf2_inc = 0.5 * ( uf2 + if2 .* ZL_a );
if2_inc = 0.5 * ( if2 + uf2 ./ ZL_a );
uf1_ref = uf1 - uf1_inc;
if1_ref = if1 - if1_inc;
uf2_ref = uf2 - uf2_inc;
if2_ref = if2 - if2_inc;
% plot s-parameter
figure
s11 = uf1_ref./uf1_inc;
s21 = uf2_inc./uf1_inc;
plot(freq,20*log10(abs(s11)),'Linewidth',2);
xlim([freq(1) freq(end)]);
xlabel('frequency (Hz)')
ylabel('s-para (dB)');
% ylim([-40 5]);
grid on;
hold on;
plot(freq,20*log10(abs(s21)),'r','Linewidth',2);
legend('s11','s21','Location','SouthEast');
% plot line-impedance comparison
figure()
ZL = uf1./if1;
plot(freq,real(ZL),'Linewidth',2);
xlim([freq(1) freq(end)]);
xlabel('frequency (Hz)')
ylabel('line-impedance (\Omega)');
grid on;
hold on;
plot(freq,imag(ZL),'r--','Linewidth',2);
plot(freq,ZL_a,'g-.','Linewidth',2);
legend('\Re\{ZL\}','\Im\{ZL\}','ZL-analytic','Location','Best');

222
matlab/examples/waveguide/Coax_CylinderCoords.m
View File

@ -1,71 +1,72 @@
%
% EXAMPLE / waveguide / coaxial cable using cylindrical coordinates
%
% This example demonstrates how to:
% - use cylindrical coordinates
% - setup a coaxial waveguide
% - use analytic functions for waveguide excitations and voltage/current
% calculations
%
%
% Tested with
% - Matlab 2009b
% - openEMS v0.0.17
%
% (C) 2010 Thorsten Liebig <thorsten.liebig@uni-due.de>
close all close all
clear clear
clc clc
%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% switches & options...
EPS0 = 8.85418781762e-12; postprocessing_only = 0;
MUE0 = 1.256637062e-6; use_pml = 0; % use pml boundaries instead of mur
C0 = 1/sqrt(EPS0*MUE0); openEMS_opts = '';
Z0 = sqrt(MUE0/EPS0); % openEMS_opts = [openEMS_opts ' --disable-dumps'];
%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
numTS = 1e5; %number of timesteps
length = 1000; %length of the waveguide
unit = 1e-3; %drawing unit used
coax_rad_i = 100; %inner radius
coax_rad_a = 230; %outer radius
mesh_res = [10 nan 10]; %desired mesh resolution
N_alpha = 71; %mesh lines in azimuth direction
physical_constants;
%excitation
f0 = 0.5e9; f0 = 0.5e9;
epsR = 1; epsR = 1;
kappa = 0;
length = 3000;
port_dist = 1500;
rad_i = 100;
rad_a = 230;
max_mesh = 10;
max_alpha = max_mesh;
N_alpha = ceil(rad_a * 2*pi / max_alpha);
mesh_res = [max_mesh 2*pi/N_alpha max_mesh];
%% define and openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
openEMS_opts = '';
openEMS_opts = [openEMS_opts ' --disable-dumps'];
% openEMS_opts = [openEMS_opts ' --debug-operator'];
% openEMS_opts = [openEMS_opts ' --debug-material'];
%% create sim path %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Sim_Path = 'tmp'; Sim_Path = 'tmp';
Sim_CSX = 'coax.xml'; Sim_CSX = 'coax.xml';
if (exist(Sim_Path,'dir')) if (postprocessing_only==0)
rmdir(Sim_Path,'s'); [status, message, messageid] = rmdir(Sim_Path,'s');
[status, message, messageid] = mkdir(Sim_Path);
end end
mkdir(Sim_Path);
%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FDTD = InitCylindricalFDTD(1e5,1e-6,'OverSampling',10); FDTD = InitCylindricalFDTD(numTS,1e-5);
FDTD = SetGaussExcite(FDTD,f0,f0); FDTD = SetGaussExcite(FDTD,f0,f0);
BC = [0 0 1 1 0 0]; BC = {'PEC','PEC','PEC','PEC','PEC','MUR'};
if (use_pml>0)
BC = {'PEC','PEC','PEC','PEC','PEC','PML_8'};
end
FDTD = SetBoundaryCond(FDTD,BC); FDTD = SetBoundaryCond(FDTD,BC);
%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CSX = InitCSX(); CSX = InitCSX('CoordSystem',1);
mesh.x = rad_i : mesh_res(1) : rad_a; mesh.x = coax_rad_i : mesh_res(1) : coax_rad_a;
mesh.y = linspace(0,2*pi,N_alpha); mesh.y = linspace(0,2*pi,N_alpha);
% mesh.y = mesh.y + mesh_res(2)/2;
mesh.z = 0 : mesh_res(3) : length; mesh.z = 0 : mesh_res(3) : length;
CSX = DefineRectGrid(CSX, 1e-3,mesh); CSX = DefineRectGrid(CSX, unit, mesh);
%% fake pml %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
abs_length = 30*(mesh.z(2)-mesh.z(1))
finalKappa = 0.3;
finalSigma = finalKappa*MUE0/EPS0/epsR;
CSX = AddMaterial(CSX,'pml');
CSX = SetMaterialProperty(CSX,'pml','Kappa',finalKappa+kappa,'Epsilon',epsR);
CSX = SetMaterialProperty(CSX,'pml','Sigma',finalSigma);
CSX = SetMaterialWeight(CSX,'pml','Kappa',['pow(abs(z)-' num2str(length-abs_length) ',4)/' num2str(abs_length^4)]);
CSX = SetMaterialWeight(CSX,'pml','Sigma',['pow(abs(z)-' num2str(length-abs_length) ',4)/' num2str(abs_length^4)]);
start = [rad_i mesh.y(1) length-abs_length];
stop = [rad_a mesh.y(end) length];
CSX = AddBox(CSX,'pml',0 ,start,stop);
%% material %% material
CSX = AddMaterial(CSX,'fill'); CSX = AddMaterial(CSX,'fill');
CSX = SetMaterialProperty(CSX,'fill','Epsilon',epsR,'Kappa',kappa); CSX = SetMaterialProperty(CSX,'fill','Epsilon',epsR);
start = [mesh.x(1) mesh.y(1) 0]; start = [mesh.x(1) mesh.y(1) 0];
stop = [mesh.x(end) mesh.y(end) length]; stop = [mesh.x(end) mesh.y(end) length];
CSX = AddBox(CSX,'fill',0 ,start,stop); CSX = AddBox(CSX,'fill',0 ,start,stop);
@ -76,8 +77,8 @@ weight{1} = '1/rho';
weight{2} = 0; weight{2} = 0;
weight{3} = 0; weight{3} = 0;
CSX = SetExcitationWeight(CSX, 'excite', weight ); CSX = SetExcitationWeight(CSX, 'excite', weight );
start = [rad_i mesh.y(1) 0]; start = [coax_rad_i mesh.y(1) 0];
stop = [rad_a mesh.y(end) 0]; stop = [coax_rad_a mesh.y(end) 0];
CSX = AddBox(CSX,'excite',0 ,start,stop); CSX = AddBox(CSX,'excite',0 ,start,stop);
%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -89,48 +90,101 @@ CSX = AddBox(CSX,'Et_',0 , start,stop);
CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',0); CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',0);
CSX = AddBox(CSX,'Ht_',0,start,stop); CSX = AddBox(CSX,'Ht_',0,start,stop);
% voltage calc (take a voltage average to be at the same spot as the %% define voltage calc boxes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% current calculation) %voltage calc
CSX = AddProbe(CSX,'ut1_1',0); CSX = AddProbe(CSX,'ut1',0);
start = [ rad_i 0 port_dist ];stop = [ rad_a 0 port_dist ]; start = [ coax_rad_i 0 mesh.z(10) ];
CSX = AddBox(CSX,'ut1_1', 0 ,start,stop); stop = [ coax_rad_a 0 mesh.z(10) ];
CSX = AddProbe(CSX,'ut1_2',0); CSX = AddBox(CSX,'ut1', 0 ,start,stop);
start = [ rad_i 0 port_dist+mesh_res(3) ];stop = [ rad_a 0 port_dist+mesh_res(3) ]; CSX = AddProbe(CSX,'ut2',0);
CSX = AddBox(CSX,'ut1_2', 0 ,start,stop); start = [ coax_rad_i 0 mesh.z(end-10)];
stop = [ coax_rad_a 0 mesh.z(end-10)];
CSX = AddBox(CSX,'ut2', 0 ,start,stop);
%current calc, for each position there are two currents, which will get
%averaged to match the voltage position in between (!Yee grid!)
CSX = AddProbe(CSX,'it1a',1);
mid = 0.5*(coax_rad_i+coax_rad_a);
start = [ 0 mesh.z(1) mesh.z(9) ];
stop = [ mid mesh.z(end) mesh.z(9) ];
CSX = AddBox(CSX,'it1a', 0 ,start,stop);
CSX = AddProbe(CSX,'it1b',1);
start = [ 0 mesh.z(1) mesh.z(10) ];
stop = [ mid mesh.z(end) mesh.z(10) ];
CSX = AddBox(CSX,'it1b', 0 ,start,stop);
CSX = AddProbe(CSX,'ut_ex',0); CSX = AddProbe(CSX,'it2a',1);
start = [ rad_i 0 0 ];stop = [ rad_a 0 0 ]; start = [ 0 mesh.z(1) mesh.z(end-11) ];
CSX = AddBox(CSX,'ut_ex', 0 ,start,stop); stop = [ mid mesh.z(end) mesh.z(end-11) ];
CSX = AddBox(CSX,'it2a', 0 ,start,stop);
CSX = AddProbe(CSX,'it2b',1);
start = [ 0 mesh.z(1) mesh.z(end-10) ];
stop = [ mid mesh.z(end) mesh.z(end-10) ];
CSX = AddBox(CSX,'it2b', 0 ,start,stop);
% current calc %% Write openEMS
CSX = AddProbe(CSX,'it1',1); if (postprocessing_only==0)
mid = 0.5*(rad_i+rad_a); WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
start = [ 0 mesh.y(1) port_dist ];stop = [ mid mesh.y(end) port_dist ]; RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts);
CSX = AddBox(CSX,'it1', 0 ,start,stop); end
%% run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%
WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); freq = linspace(0,2*f0,201);
RunOpenEMS(Sim_Path,Sim_CSX,openEMS_opts); U = ReadUI({'ut1','ut2'},[Sim_Path '/'],freq);
I = ReadUI({'it1a','it1b','it2a','it2b'},[Sim_Path '/'],freq);
%% postproc & do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Exc = ReadUI('et',Sim_Path,freq);
UI = ReadUI({'ut1_1','ut1_2','it1'},Sim_Path);
plot(UI.TD{1}.t,UI.TD{1}.val)
UI_ex = ReadUI({'ut_ex'},'tmp/');
hold on;
plot(UI_ex.TD{1}.t,UI_ex.TD{1}.val,'r--');
u_f = (UI.FD{1}.val + UI.FD{2}.val)/2; %averaging voltages to fit current
i_f = UI.FD{3}.val;
%% plot voltages
figure figure
ZL = Z0/2/pi/sqrt(epsR)*log(rad_a/rad_i); %analytic line-impedance of a coax plot(U.TD{1}.t, U.TD{1}.val,'Linewidth',2);
plot(UI.FD{1}.f,ZL*ones(size(u_f)),'g');
hold on; hold on;
grid on; grid on;
Z = u_f./i_f; plot(U.TD{2}.t, U.TD{2}.val,'r--','Linewidth',2);
plot(UI.FD{1}.f,real(Z),'Linewidth',2); xlabel('time (s)')
plot(UI.FD{1}.f,imag(Z),'r','Linewidth',2); ylabel('voltage (V)')
xlim([0 2*f0]); legend('u_1(t)','u_2(t)')
legend('Z_L','\Re\{Z\}','\Im\{Z\}','Location','Best');
%% calculate incoming and reflected voltages & currents
ZL_a = ones(size(freq))*Z0/2/pi/sqrt(epsR)*log(coax_rad_a/coax_rad_i); %analytic line-impedance of a coax
uf1 = U.FD{1}.val./Exc.FD{1}.val;
uf2 = U.FD{2}.val./Exc.FD{1}.val;
if1 = 0.5*(I.FD{1}.val+I.FD{2}.val)./Exc.FD{1}.val;
if2 = 0.5*(I.FD{3}.val+I.FD{4}.val)./Exc.FD{1}.val;
uf1_inc = 0.5 * ( uf1 + if1 .* ZL_a );
if1_inc = 0.5 * ( if1 + uf1 ./ ZL_a );
uf2_inc = 0.5 * ( uf2 + if2 .* ZL_a );
if2_inc = 0.5 * ( if2 + uf2 ./ ZL_a );
uf1_ref = uf1 - uf1_inc;
if1_ref = if1 - if1_inc;
uf2_ref = uf2 - uf2_inc;
if2_ref = if2 - if2_inc;
% plot s-parameter
figure
s11 = uf1_ref./uf1_inc;
s21 = uf2_inc./uf1_inc;
plot(freq,20*log10(abs(s11)),'Linewidth',2);
xlim([freq(1) freq(end)]);
xlabel('frequency (Hz)')
ylabel('s-para (dB)');
% ylim([-40 5]);
grid on;
hold on;
plot(freq,20*log10(abs(s21)),'r','Linewidth',2);
legend('s11','s21','Location','SouthEast');
% plot line-impedance comparison
figure()
ZL = uf1./if1;
plot(freq,real(ZL),'Linewidth',2);
xlim([freq(1) freq(end)]);
xlabel('frequency (Hz)')
ylabel('line-impedance (\Omega)');
grid on;
hold on;
plot(freq,imag(ZL),'r--','Linewidth',2);
plot(freq,ZL_a,'g-.','Linewidth',2);
legend('\Re\{ZL\}','\Im\{ZL\}','ZL-analytic','Location','Best');

38
matlab/examples/waveguide/Rect_Waveguide.m
View File

@ -18,7 +18,7 @@ clear
clc clc
%% switches %% switches
postproc_only = 1; postproc_only = 0;
%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
physical_constants; physical_constants;
@ -26,7 +26,7 @@ unit = 1e-3; %drawing unit in mm
numTS = 50000; %max. number of timesteps numTS = 50000; %max. number of timesteps
% waveguide dimensions % waveguide dimensions
length = 2000; length = 1000;
a = 1000; %waveguide width a = 1000; %waveguide width
b = 600; %waveguide heigth b = 600; %waveguide heigth
@ -79,10 +79,8 @@ Sim_Path = 'tmp';
Sim_CSX = 'rect_wg.xml'; Sim_CSX = 'rect_wg.xml';
if (postproc_only==0) if (postproc_only==0)
if (exist(Sim_Path,'dir')) [status, message, messageid] = rmdir(Sim_Path,'s');
rmdir(Sim_Path,'s'); [status, message, messageid] = mkdir(Sim_Path);
end
mkdir(Sim_Path);
end end
%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -124,6 +122,8 @@ CSX = AddProbe(CSX, 'ut2', 10, 1, [], 'ModeFunction',{func_Ex,func_Ey,0});
CSX = AddBox(CSX, 'ut2', 0 ,start,stop); CSX = AddBox(CSX, 'ut2', 0 ,start,stop);
CSX = AddProbe(CSX,'it2', 11, 1, [], 'ModeFunction',{func_Hx,func_Hy,0}); CSX = AddProbe(CSX,'it2', 11, 1, [], 'ModeFunction',{func_Hx,func_Hy,0});
CSX = AddBox(CSX,'it2', 0 ,start,stop); CSX = AddBox(CSX,'it2', 0 ,start,stop);
port_dist = mesh.z(end-15) - mesh.z(15);
%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CSX = AddDump(CSX,'Et','FileType',1,'SubSampling','4,4,2'); CSX = AddDump(CSX,'Et','FileType',1,'SubSampling','4,4,2');
@ -188,6 +188,20 @@ xlabel('freq (Hz)');
xlim([freq(1) freq(end)]); xlim([freq(1) freq(end)]);
legend('\Re(Z_L)','\Im(Z_L)','Z_L analytic','Location','Best'); legend('\Re(Z_L)','\Im(Z_L)','Z_L analytic','Location','Best');
%% beta compare
figure()
da = unwrap(angle(uf1_inc./uf2_inc)) ;
% da = mod(da,2*pi);
beta_12 = (da)/port_dist/unit;
plot(freq,beta_12,'Linewidth',2);
xlim([freq(1) freq(end)]);
xlabel('frequency (Hz)');
ylabel('\beta (m^{-1})');
grid on;
hold on;
plot(freq,beta,'g--','Linewidth',2);
legend('\beta-FDTD','\beta-analytic','Location','Best');
%% Plot the field dumps %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Plot the field dumps %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dump_file = [Sim_Path '/Et.h5']; dump_file = [Sim_Path '/Et.h5'];
figure() figure()
@ -200,17 +214,15 @@ PlotHDF5FieldData(dump_file, PlotArgs)
%% dump frequency to vtk %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% dump frequency to vtk %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cleanup and create dump folder % cleanup and create dump folder
vtk_path = [Sim_Path '/vtk']; vtk_path = [Sim_Path '/vtk'];
if exist(vtk_path,'dir') [status, message, messageid] = rmdir(vtk_path,'s');
rmdir(vtk_path,'s'); [status, message, messageid] = mkdir(vtk_path);
end
mkdir(vtk_path);
disp('Dumping to vtk files... this may take a minute...') disp('Dumping to vtk files... this may take a minute...')
% define interpolation mesh % define interpolation mesh
mesh_interp{1}=mesh.x * unit; mesh_interp{1}=mesh.x * unit;
mesh_interp{2}=b/2 * unit; mesh_interp{2}=b/2 * unit;
mesh_interp{3}=mesh.z * unit; mesh_interp{3}=mesh.z * unit;
[field_FD mesh_FD] = ReadHDF5Dump(dump_file,'Interpolation',mesh_interp,'Frequency',vtk_dump_freq); [field mesh_FD] = ReadHDF5Dump(dump_file,'Interpolation',mesh_interp,'Frequency',vtk_dump_freq);
% dump animated phase to vtk % dump animated phase to vtk
for n=1:numel(vtk_dump_freq) for n=1:numel(vtk_dump_freq)
@ -218,11 +230,11 @@ for n=1:numel(vtk_dump_freq)
phase = phase(1:end-1); phase = phase(1:end-1);
for ph = phase for ph = phase
filename = [vtk_path '/E_xz_f=' num2str(vtk_dump_freq(n)) '_p' num2str(ph) '.vtk']; filename = [vtk_path '/E_xz_f=' num2str(vtk_dump_freq(n)) '_p' num2str(ph) '.vtk'];
Dump2VTK(filename,real(field_FD.values{n}.*exp(1j*ph/180*pi)),mesh_FD,'E-Field'); Dump2VTK(filename,real(field.FD.values{n}.*exp(1j*ph/180*pi)),mesh_FD,'E-Field');
end end
filename = [vtk_path '/E_xz_f=' num2str(vtk_dump_freq(n)) '_mag.vtk']; filename = [vtk_path '/E_xz_f=' num2str(vtk_dump_freq(n)) '_mag.vtk'];
Dump2VTK(filename,abs(field_FD.values{n}),mesh_FD,'E-Field'); Dump2VTK(filename,abs(field.FD.values{n}),mesh_FD,'E-Field');
end end
disp('done... you can open and visualize the vtk-files using Paraview (www.paraview.org)!') disp('done... you can open and visualize the vtk-files using Paraview (www.paraview.org)!')