From dce988b52f940bfe5db48eb91f6e7e6efc467636 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Sun, 2 Jan 2011 11:32:00 +0100 Subject: [PATCH] waveguide examples update --- matlab/examples/waveguide/Coax.m | 192 ++++++++++----- .../examples/waveguide/Coax_CylinderCoords.m | 222 +++++++++++------- matlab/examples/waveguide/Rect_Waveguide.m | 38 ++- 3 files changed, 294 insertions(+), 158 deletions(-) diff --git a/matlab/examples/waveguide/Coax.m b/matlab/examples/waveguide/Coax.m index 629d063..d7ce8a7 100644 --- a/matlab/examples/waveguide/Coax.m +++ b/matlab/examples/waveguide/Coax.m @@ -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 + close all clear 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -abs_length = 250; -length = 1000; -coax_rad_i = 100; -coax_rad_ai = 230; -coax_rad_aa = 240; -mesh_res = [5 5 5]; +numTS = 5000; %number of timesteps +length = 1000; %length of the waveguide +unit = 1e-3; %drawing unit used +coax_rad_i = 100; %inner radius +coax_rad_ai = 230; %inner radius of outer cladding +coax_rad_aa = 240; %outer radius of outer cladding +mesh_res = [5 5 5]; %desired mesh resolution -EPS0 = 8.85418781762e-12; -MUE0 = 1.256637062e-6; -C0 = 1/sqrt(EPS0*MUE0); -Z0 = sqrt(MUE0/EPS0); +physical_constants; +%excitation f0 = 0.5e9; epsR = 1; -%% define openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -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']; - +%% create sim path %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Sim_Path = 'tmp'; Sim_CSX = 'coax.xml'; -if (exist(Sim_Path,'dir')) - rmdir(Sim_Path,'s'); +if (postprocessing_only==0) + [status, message, messageid] = rmdir(Sim_Path,'s'); + [status, message, messageid] = mkdir(Sim_Path); end -mkdir(Sim_Path); %% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% -FDTD = InitFDTD(5e5,1e-5); +FDTD = InitFDTD(numTS,1e-5); 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); %% 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.y = mesh.x; mesh.z = 0 : mesh_res(3) : length; -CSX = DefineRectGrid(CSX, 1e-3,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)']); +CSX = DefineRectGrid(CSX, unit, mesh); %%% coax CSX = AddMaterial(CSX,'copper'); @@ -62,11 +69,10 @@ CSX = SetMaterialProperty(CSX,'copper','Kappa',56e6); start = [0, 0 , 0];stop = [0, 0 , length]; 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)); -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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -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]); weight{1} = '(x)/(x*x+y*y)'; 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 = AddBox(CSX,'Ht_',0,start,stop); +%% define voltage calc boxes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %voltage calc -CSX = AddProbe(CSX,'ut1_1',0); -start = [ coax_rad_i 0 length/2 ];stop = [ coax_rad_ai 0 length/2 ]; -CSX = AddBox(CSX,'ut1_1', 0 ,start,stop); -CSX = AddProbe(CSX,'ut1_2',0); -start = [ coax_rad_i 0 length/2+mesh_res(3) ];stop = [ coax_rad_ai 0 length/2+mesh_res(3) ]; -CSX = AddBox(CSX,'ut1_2', 0 ,start,stop); +CSX = AddProbe(CSX,'ut1',0); +start = [ coax_rad_i 0 mesh.z(10) ]; +stop = [ coax_rad_ai 0 mesh.z(10) ]; +CSX = AddBox(CSX,'ut1', 0 ,start,stop); +CSX = AddProbe(CSX,'ut2',0); +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 -CSX = AddProbe(CSX,'it1',1); +%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_ai); -start = [ -mid -mid length/2 ];stop = [ mid mid length/2 ]; -CSX = AddBox(CSX,'it1', 0 ,start,stop); +start = [ -mid -mid mesh.z(9) ]; +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 -WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); +CSX = AddProbe(CSX,'it2a',1); +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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -RunOpenEMS(Sim_Path,Sim_CSX,openEMS_opts); +%% Write openEMS +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); -u_f = (UI.FD{1}.val + UI.FD{2}.val)/2; %averaging voltages to fit current -i_f = UI.FD{3}.val; +%% +freq = linspace(0,2*f0,201); +U = ReadUI({'ut1','ut2'},[Sim_Path '/'],freq); +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(UI.FD{1}.f,ZL*ones(size(u_f)),'g'); +%% plot voltages +figure +plot(U.TD{1}.t, U.TD{1}.val,'Linewidth',2); hold on; grid on; -Z = u_f./i_f; -plot(UI.FD{1}.f,real(Z),'Linewidth',2); -plot(UI.FD{1}.f,imag(Z),'r','Linewidth',2); -xlim([0 2*f0]); -legend('Z_L','\Re\{Z\}','\Im\{Z\}','Location','Best'); +plot(U.TD{2}.t, U.TD{2}.val,'r--','Linewidth',2); +xlabel('time (s)') +ylabel('voltage (V)') +legend('u_1(t)','u_2(t)') + +%% 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'); diff --git a/matlab/examples/waveguide/Coax_CylinderCoords.m b/matlab/examples/waveguide/Coax_CylinderCoords.m index 70f80f4..5e4606c 100644 --- a/matlab/examples/waveguide/Coax_CylinderCoords.m +++ b/matlab/examples/waveguide/Coax_CylinderCoords.m @@ -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 + close all clear clc -%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -EPS0 = 8.85418781762e-12; -MUE0 = 1.256637062e-6; -C0 = 1/sqrt(EPS0*MUE0); -Z0 = sqrt(MUE0/EPS0); +%% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +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; 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_CSX = 'coax.xml'; -if (exist(Sim_Path,'dir')) - rmdir(Sim_Path,'s'); +if (postprocessing_only==0) + [status, message, messageid] = rmdir(Sim_Path,'s'); + [status, message, messageid] = mkdir(Sim_Path); end -mkdir(Sim_Path); %% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% -FDTD = InitCylindricalFDTD(1e5,1e-6,'OverSampling',10); +FDTD = InitCylindricalFDTD(numTS,1e-5); 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); %% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -CSX = InitCSX(); -mesh.x = rad_i : mesh_res(1) : rad_a; +CSX = InitCSX('CoordSystem',1); +mesh.x = coax_rad_i : mesh_res(1) : coax_rad_a; mesh.y = linspace(0,2*pi,N_alpha); -% mesh.y = mesh.y + mesh_res(2)/2; mesh.z = 0 : mesh_res(3) : length; -CSX = DefineRectGrid(CSX, 1e-3,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); +CSX = DefineRectGrid(CSX, unit, mesh); %% material 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]; stop = [mesh.x(end) mesh.y(end) length]; CSX = AddBox(CSX,'fill',0 ,start,stop); @@ -76,8 +77,8 @@ weight{1} = '1/rho'; weight{2} = 0; weight{3} = 0; CSX = SetExcitationWeight(CSX, 'excite', weight ); -start = [rad_i mesh.y(1) 0]; -stop = [rad_a mesh.y(end) 0]; +start = [coax_rad_i mesh.y(1) 0]; +stop = [coax_rad_a mesh.y(end) 0]; CSX = AddBox(CSX,'excite',0 ,start,stop); %% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -89,48 +90,101 @@ CSX = AddBox(CSX,'Et_',0 , start,stop); CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',0); CSX = AddBox(CSX,'Ht_',0,start,stop); -% voltage calc (take a voltage average to be at the same spot as the -% current calculation) -CSX = AddProbe(CSX,'ut1_1',0); -start = [ rad_i 0 port_dist ];stop = [ rad_a 0 port_dist ]; -CSX = AddBox(CSX,'ut1_1', 0 ,start,stop); -CSX = AddProbe(CSX,'ut1_2',0); -start = [ rad_i 0 port_dist+mesh_res(3) ];stop = [ rad_a 0 port_dist+mesh_res(3) ]; -CSX = AddBox(CSX,'ut1_2', 0 ,start,stop); +%% define voltage calc boxes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%voltage calc +CSX = AddProbe(CSX,'ut1',0); +start = [ coax_rad_i 0 mesh.z(10) ]; +stop = [ coax_rad_a 0 mesh.z(10) ]; +CSX = AddBox(CSX,'ut1', 0 ,start,stop); +CSX = AddProbe(CSX,'ut2',0); +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); -start = [ rad_i 0 0 ];stop = [ rad_a 0 0 ]; -CSX = AddBox(CSX,'ut_ex', 0 ,start,stop); +CSX = AddProbe(CSX,'it2a',1); +start = [ 0 mesh.z(1) mesh.z(end-11) ]; +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 -CSX = AddProbe(CSX,'it1',1); -mid = 0.5*(rad_i+rad_a); -start = [ 0 mesh.y(1) port_dist ];stop = [ mid mesh.y(end) port_dist ]; -CSX = AddBox(CSX,'it1', 0 ,start,stop); +%% Write openEMS +if (postprocessing_only==0) + WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts); +end -%% run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); -RunOpenEMS(Sim_Path,Sim_CSX,openEMS_opts); - -%% postproc & do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -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; +%% +freq = linspace(0,2*f0,201); +U = ReadUI({'ut1','ut2'},[Sim_Path '/'],freq); +I = ReadUI({'it1a','it1b','it2a','it2b'},[Sim_Path '/'],freq); +Exc = ReadUI('et',Sim_Path,freq); +%% plot voltages figure -ZL = Z0/2/pi/sqrt(epsR)*log(rad_a/rad_i); %analytic line-impedance of a coax -plot(UI.FD{1}.f,ZL*ones(size(u_f)),'g'); +plot(U.TD{1}.t, U.TD{1}.val,'Linewidth',2); hold on; grid on; -Z = u_f./i_f; -plot(UI.FD{1}.f,real(Z),'Linewidth',2); -plot(UI.FD{1}.f,imag(Z),'r','Linewidth',2); -xlim([0 2*f0]); -legend('Z_L','\Re\{Z\}','\Im\{Z\}','Location','Best'); +plot(U.TD{2}.t, U.TD{2}.val,'r--','Linewidth',2); +xlabel('time (s)') +ylabel('voltage (V)') +legend('u_1(t)','u_2(t)') + +%% 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'); diff --git a/matlab/examples/waveguide/Rect_Waveguide.m b/matlab/examples/waveguide/Rect_Waveguide.m index 212fb01..a9601ad 100644 --- a/matlab/examples/waveguide/Rect_Waveguide.m +++ b/matlab/examples/waveguide/Rect_Waveguide.m @@ -18,7 +18,7 @@ clear clc %% switches -postproc_only = 1; +postproc_only = 0; %% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% physical_constants; @@ -26,7 +26,7 @@ unit = 1e-3; %drawing unit in mm numTS = 50000; %max. number of timesteps % waveguide dimensions -length = 2000; +length = 1000; a = 1000; %waveguide width b = 600; %waveguide heigth @@ -79,10 +79,8 @@ Sim_Path = 'tmp'; Sim_CSX = 'rect_wg.xml'; if (postproc_only==0) - if (exist(Sim_Path,'dir')) - rmdir(Sim_Path,'s'); - end - mkdir(Sim_Path); + [status, message, messageid] = rmdir(Sim_Path,'s'); + [status, message, messageid] = mkdir(Sim_Path); end %% 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 = AddProbe(CSX,'it2', 11, 1, [], 'ModeFunction',{func_Hx,func_Hy,0}); CSX = AddBox(CSX,'it2', 0 ,start,stop); + +port_dist = mesh.z(end-15) - mesh.z(15); %% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CSX = AddDump(CSX,'Et','FileType',1,'SubSampling','4,4,2'); @@ -188,6 +188,20 @@ xlabel('freq (Hz)'); xlim([freq(1) freq(end)]); 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% dump_file = [Sim_Path '/Et.h5']; figure() @@ -200,17 +214,15 @@ PlotHDF5FieldData(dump_file, PlotArgs) %% dump frequency to vtk %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % cleanup and create dump folder vtk_path = [Sim_Path '/vtk']; -if exist(vtk_path,'dir') - rmdir(vtk_path,'s'); -end -mkdir(vtk_path); +[status, message, messageid] = rmdir(vtk_path,'s'); +[status, message, messageid] = mkdir(vtk_path); disp('Dumping to vtk files... this may take a minute...') % define interpolation mesh mesh_interp{1}=mesh.x * unit; mesh_interp{2}=b/2 * 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 for n=1:numel(vtk_dump_freq) @@ -218,11 +230,11 @@ for n=1:numel(vtk_dump_freq) phase = phase(1:end-1); for ph = phase 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 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 disp('done... you can open and visualize the vtk-files using Paraview (www.paraview.org)!')