1%
2% EXAMPLE / waveguide / circular waveguide
3%
4% This example demonstrates how to:
5%  - setup a circular waveguide
6%  - use analytic functions for waveguide excitations and voltage/current
7%  calculations
8%
9%
10% Tested with
11%  - Matlab 2009b
12%  - openEMS v0.0.17
13%
14% (C) 2010 Thorsten Liebig <thorsten.liebig@uni-due.de>
15
16close all
17clear
18clc
19
20%% switches & options...
21postprocessing_only = 0;
22use_pml = 0;            % use pml boundaries instead of mur
23openEMS_opts = '';
24% openEMS_opts = [openEMS_opts ' --disable-dumps'];
25
26%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
27numTS = 1e5;            %number of timesteps
28length = 1000;          %length of the waveguide
29unit = 1e-3;            %drawing unit used
30rad  = 300;             %radius of the circular waveguide
31mesh_res = [10 10 15];  %desired mesh resolution
32
33%excitation
34f0 = 350e6;             %center frequency
35f0_BW = 25e6;           %bandwidth: 10dB cut-off frequency
36
37physical_constants
38
39%% TE11 mode definitions (Pozar 3rd edition)
40p11 = 1.841;
41kc = p11 / rad /unit;
42k = 2*pi*f0/C0;
43fc = C0*kc/2/pi;
44beta = sqrt(k^2 - kc^2);
45n_eff = (beta/k);
46
47kc = kc*unit; %functions must be defined in drawing units
48func_Er = [ num2str(-1/kc^2) '/rho*cos(a)*j1('  num2str(kc) '*rho)'];
49func_Ea = [ num2str(1/kc) '*sin(a)*0.5*(j0('  num2str(kc) '*rho)-jn(2,'  num2str(kc) '*rho))'];
50func_Ex = ['(' func_Er '*cos(a) - ' func_Ea '*sin(a) )*(rho<' num2str(rad) ')'];
51func_Ey = ['(' func_Er '*sin(a) + ' func_Ea '*cos(a) )*(rho<' num2str(rad) ')'];
52
53func_Ha = [ num2str(-1/kc^2,'%14.13f') '/rho*cos(a)*j1('  num2str(kc,'%14.13f') '*rho)'];
54func_Hr = [ '-1*' num2str(1/kc,'%14.13f') '*sin(a)*0.5*(j0('  num2str(kc,'%14.13f') '*rho)-jn(2,'  num2str(kc,'%14.13f') '*rho))'];
55func_Hx = ['(' func_Hr '*cos(a) - ' func_Ha '*sin(a) )*(rho<' num2str(rad) ')'];
56func_Hy = ['(' func_Hr '*sin(a) + ' func_Ha '*cos(a) )*(rho<' num2str(rad) ')'];
57
58%% define files and path %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
59Sim_Path = 'tmp';
60Sim_CSX = 'Circ_WG.xml';
61
62if (postprocessing_only==0)
63    [status, message, messageid] = rmdir(Sim_Path,'s');
64    [status, message, messageid] = mkdir(Sim_Path);
65end
66
67%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
68FDTD = InitFDTD(numTS,1e-6,'OverSampling',5);
69FDTD = SetGaussExcite(FDTD,f0,f0_BW);
70BC = {'PEC','PEC','PEC','PEC','PEC','MUR'};
71if (use_pml>0)
72    BC = {'PEC','PEC','PEC','PEC','PEC','PML_8'};
73end
74FDTD = SetBoundaryCond(FDTD,BC,'MUR_PhaseVelocity',C0 / n_eff);
75
76%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77CSX = InitCSX();
78mesh.x = -mesh_res(1)/2-rad:mesh_res(1):rad+mesh_res(1)/2;
79mesh.y = -mesh_res(2)/2-rad:mesh_res(2):rad+mesh_res(2)/2;
80mesh.z = 0 : mesh_res(3) : length;
81CSX = DefineRectGrid(CSX, 1e-3,mesh);
82
83start = [0,0,0];
84stop = [0,0,length];
85
86%%% fill everything with copper, priority 0
87CSX = AddMetal(CSX,'copper');
88% CSX = SetMaterialProperty(CSX,'copper','Kappa',56e6);
89CSX = AddBox(CSX,'copper',0,[mesh.x(1) mesh.y(1) mesh.z(1)],[mesh.x(end) mesh.y(end) mesh.z(end)]);
90
91%%% cut out an air cylinder as circular waveguide... priority 5
92CSX = AddMaterial(CSX,'air');
93CSX = SetMaterialProperty(CSX,'air','Epsilon',1);
94CSX = AddCylinder(CSX,'air', 5 ,start,stop,rad);
95
96CSX = AddExcitation(CSX,'excite',0,[1 1 0]);
97weight{1} = func_Ex;
98weight{2} = func_Ey;
99weight{3} = 0;
100CSX = SetExcitationWeight(CSX, 'excite', weight );
101CSX = AddCylinder(CSX,'excite', 5 ,[0 0 -0.1],[0 0 0.1],rad);
102
103%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
104CSX = AddDump(CSX,'Et_','SubSampling','2,2,2','FileType',0,'DumpMode',2);
105start = [mesh.x(1) , 0 , mesh.z(1)];
106stop = [mesh.x(end), 0 , mesh.z(end)];
107CSX = AddBox(CSX,'Et_',0 , start,stop);
108
109CSX = AddDump(CSX,'Ht_','SubSampling','2,2,2','DumpType',1,'FileType',0,'DumpMode',2);
110CSX = AddBox(CSX,'Ht_',0,start,stop);
111
112%% define voltage calc boxes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
113%voltage calc
114start = [mesh.x(1)   mesh.y(1)   mesh.z(10)];
115stop  = [mesh.x(end) mesh.y(end) mesh.z(10)];
116CSX = AddProbe(CSX, 'ut1', 10, 1, [], 'ModeFunction',{func_Ex,func_Ey,0});
117CSX = AddBox(CSX,  'ut1',  0 ,start,stop);
118CSX = AddProbe(CSX,'it1', 11, 1, [], 'ModeFunction',{func_Hx,func_Hy,0});
119CSX = AddBox(CSX,'it1', 0 ,start,stop);
120
121start = [mesh.x(1)   mesh.y(1)   mesh.z(end-10)];
122stop  = [mesh.x(end) mesh.y(end) mesh.z(end-10)];
123CSX = AddProbe(CSX, 'ut2', 10, 1, [], 'ModeFunction',{func_Ex,func_Ey,0});
124CSX = AddBox(CSX,  'ut2',  0 ,start,stop);
125CSX = AddProbe(CSX,'it2', 11, 1, [], 'ModeFunction',{func_Hx,func_Hy,0});
126CSX = AddBox(CSX,'it2', 0 ,start,stop);
127
128port_dist = mesh.z(end-10) - mesh.z(10);
129
130%% Write openEMS
131if (postprocessing_only==0)
132    WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
133
134    RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts);
135end
136
137%% do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
138freq = linspace(f0-f0_BW,f0+f0_BW,201);
139U = ReadUI({'ut1','ut2'},[Sim_Path '/'],freq);
140I = ReadUI({'it1','it2'},[Sim_Path '/'],freq);
141Exc = ReadUI('et',Sim_Path,freq);
142
143k = 2*pi*freq/C0;
144kc = p11 / rad /unit;
145beta = sqrt(k.^2 - kc^2);
146
147ZL_a = Z0*k./beta ;
148
149uf1 = U.FD{1}.val./Exc.FD{1}.val;
150uf2 = U.FD{2}.val./Exc.FD{1}.val;
151if1 = I.FD{1}.val./Exc.FD{1}.val;
152if2 = I.FD{2}.val./Exc.FD{1}.val;
153
154uf1_inc = 0.5 * ( uf1 + if1 .* ZL_a );
155if1_inc = 0.5 * ( if1 + uf1 ./ ZL_a );
156uf2_inc = 0.5 * ( uf2 + if2 .* ZL_a );
157if2_inc = 0.5 * ( if2 + uf2 ./ ZL_a );
158
159uf1_ref = uf1 - uf1_inc;
160if1_ref = if1 - if1_inc;
161uf2_ref = uf2 - uf2_inc;
162if2_ref = if2 - if2_inc;
163
164% plot s-parameter
165figure
166s11 = uf1_ref./uf1_inc;
167s21 = uf2_inc./uf1_inc;
168plot(freq,20*log10(abs(s11)),'Linewidth',2);
169xlim([freq(1) freq(end)]);
170xlabel('frequency (Hz)')
171ylabel('s-para (dB)');
172% ylim([-40 5]);
173grid on;
174hold on;
175plot(freq,20*log10(abs(s21)),'r','Linewidth',2);
176legend('s11','s21','Location','SouthEast');
177
178% plot line-impedance comparison
179figure()
180ZL = uf1./if1;
181plot(freq,real(ZL),'Linewidth',2);
182xlim([freq(1) freq(end)]);
183xlabel('frequency (Hz)')
184ylabel('line-impedance (\Omega)');
185grid on;
186hold on;
187plot(freq,imag(ZL),'r--','Linewidth',2);
188plot(freq,ZL_a,'g-.','Linewidth',2);
189legend('\Re\{ZL\}','\Im\{ZL\}','ZL-analytic','Location','Best');
190
191% beta compare
192figure()
193da = angle(uf1_inc)-angle(uf2_inc);
194da = mod(da,2*pi);
195beta_12 = (da)/port_dist/unit;
196plot(freq,beta_12,'Linewidth',2);
197xlim([freq(1) freq(end)]);
198xlabel('frequency (Hz)');
199ylabel('\beta (m^{-1})');
200grid on;
201hold on;
202plot(freq,beta,'g--','Linewidth',2);
203legend('\beta-FDTD','\beta-analytic','Location','Best');
204
205%% visualize electric & magnetic fields
206disp('you will find vtk dump files in the simulation folder (tmp/)')
207disp('use paraview to visulaize them');
208