Tutorials: use new waveguide ports

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>

matlab/Tutorials/Circ_Waveguide.m

@@ -6,9 +6,9 @@
 %
 % Tested with
 %  - Matlab 2011a / Octave 3.4.3
-%  - openEMS v0.0.26
+%  - openEMS v0.0.31
 %
-% (C) 2010-2012 Thorsten Liebig <thorsten.liebig@gmx.de>
+% (C) 2010-2013 Thorsten Liebig <thorsten.liebig@gmx.de>
 
 close all
 clear
@@ -28,33 +28,8 @@ f_stop  =  500e6;
 
 mesh_res = [10 2*pi/49.999 10]; %targeted mesh resolution
 
-%% mode functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% by David M. Pozar, Microwave Engineering, third edition
-freq = linspace(f_start,f_stop,201);
-p11 = 1.841;
-kc = p11 / rad /unit;
-k = 2*pi*freq/C0;
-fc = C0*kc/2/pi;
-beta = sqrt(k.^2 - kc^2);
-n_eff = (beta/k);
-ZL_a = k * Z0 ./ beta;    %analytic waveguide impedance
-
-% TE_11 mode profile E- and H-field
-kc = kc*unit; %functions must be defined in drawing units
-func_Er = [ num2str(-1/kc^2,15) '/rho*cos(a)*j1('  num2str(kc,15) '*rho)'];
-func_Ea = [ num2str(1/kc,15) '*sin(a)*0.5*(j0('  num2str(kc,15) '*rho)-jn(2,'  num2str(kc,15) '*rho))'];
-
-func_Ha = [ num2str(-1/kc^2,15) '/rho*cos(a)*j1('  num2str(kc,15) '*rho)'];
-func_Hr = [ '-1*' num2str(1/kc,15) '*sin(a)*0.5*(j0('  num2str(kc,15) '*rho)-jn(2,'  num2str(kc,15) '*rho))'];
-
-disp([' Cutoff frequencies for this mode and wavguide is: ' num2str(fc/1e6) ' MHz']);
-
-if (f_start<fc)
-    warning('openEMS:example','f_start is smaller than the cutoff-frequency, this may result in a long simulation... ');
-end
-
 %% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-FDTD = InitFDTD(1e6,1e-5,'CoordSystem',1);
+FDTD = InitFDTD('EndCriteria',1e-4,'CoordSystem',1);
 FDTD = SetGaussExcite(FDTD,0.5*(f_start+f_stop),0.5*(f_stop-f_start));
 
 % boundary conditions
@@ -68,35 +43,14 @@ mesh.a = SmoothMeshLines([0 2*pi], mesh_res(2)); % mesh in aziumthal dir.
 mesh.z = SmoothMeshLines([0 length], mesh_res(3));
 CSX = DefineRectGrid(CSX, unit,mesh);
 
-%% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% apply the waveguide port %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% ra-mode profile excitation located directly on top of pml (first 8 z-lines)
+start=[mesh.r(1)   mesh.a(1)   mesh.z(8)];
-CSX = AddExcitation(CSX,'excite',0,[1 1 0]);
+stop =[mesh.r(end) mesh.a(end) mesh.z(15)];
-weight{1} = func_Er;
+[CSX, port{1}] = AddCircWaveGuidePort( CSX, 0, 1, start, stop, rad*unit, 'TE11', 0, 1);
-weight{2} = func_Ea;
-weight{3} = 0;
-CSX = SetExcitationWeight(CSX,'excite',weight);
-start=[mesh.r(1)   mesh.a(1)   mesh.z(8) ];
-stop =[mesh.r(end) mesh.a(end) mesh.z(8) ];
-CSX = AddBox(CSX,'excite',0 ,start,stop);
 
-%% voltage and current definitions using the mode matching probes %%%%%%%%%
+start=[mesh.r(1)   mesh.a(1)   mesh.z(end-13)];
-%port 1
+stop =[mesh.r(end) mesh.a(end) mesh.z(end-14)];
-start = [mesh.r(1)   mesh.a(1)   mesh.z(15)];
+[CSX, port{2}] = AddCircWaveGuidePort( CSX, 0, 2, start, stop, rad*unit, 'TE11');
-stop  = [mesh.r(end) mesh.a(end) mesh.z(15)];
-CSX = AddProbe(CSX, 'ut1', 10, 'ModeFunction',{func_Er,func_Ea,0});
-CSX = AddBox(CSX,  'ut1',  0 ,start,stop);
-CSX = AddProbe(CSX,'it1', 11, 'ModeFunction',{func_Hr,func_Ha,0});
-CSX = AddBox(CSX,'it1', 0 ,start,stop);
-
-%port 2
-start(3) = mesh.z(end-14);
-stop(3) = mesh.z(end-14);
-CSX = AddProbe(CSX, 'ut2', 10, 'ModeFunction',{func_Er,func_Ea,0});
-CSX = AddBox(CSX,  'ut2',  0 ,start,stop);
-CSX = AddProbe(CSX,'it2', 11, 'ModeFunction',{func_Hr,func_Ha,0});
-CSX = AddBox(CSX,'it2', 0 ,start,stop);
-
-port_dist = mesh.z(end-14) - mesh.z(15);
 
 %% define dump box... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 CSX = AddDump(CSX,'Et','FileType',1,'SubSampling','4,4,4');
@@ -106,7 +60,7 @@ CSX = AddBox(CSX,'Et',0 , start,stop);
 
 %% Write openEMS compatoble xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Sim_Path = 'tmp';
-Sim_CSX = 'rect_wg.xml';
+Sim_CSX = 'circ_wg.xml';
 
 [status, message, messageid] = rmdir(Sim_Path,'s');
 [status, message, messageid] = mkdir(Sim_Path);
@@ -116,29 +70,16 @@ WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
 RunOpenEMS(Sim_Path, Sim_CSX)
 
 %% postproc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-U = ReadUI({'ut1','ut2'},[Sim_Path '/'],freq);
+freq = linspace(f_start,f_stop,201);
-I = ReadUI({'it1','it2'},[Sim_Path '/'],freq);
+port = calcPort( port, Sim_Path, freq);
-Exc = ReadUI('et',Sim_Path,freq);
 
-uf1 = U.FD{1}.val./Exc.FD{1}.val;
+s11 = port{1}.uf.ref./ port{1}.uf.inc;
-uf2 = U.FD{2}.val./Exc.FD{1}.val;
+s21 = port{2}.uf.ref./ port{1}.uf.inc;
-if1 = I.FD{1}.val./Exc.FD{1}.val;
+ZL = port{1}.uf.tot./port{1}.if.tot;
-if2 = I.FD{2}.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_inc - if1; 
-uf2_ref = uf2 - uf2_inc;
-if2_ref = if2_inc - if2;
 
 %% plot s-parameter %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 figure
-s11 = uf1_ref./uf1_inc;
-s21 = uf2_inc./uf1_inc;
 plot(freq*1e-6,20*log10(abs(s11)),'k-','Linewidth',2);
 xlim([freq(1) freq(end)]*1e-6);
 grid on;
@@ -151,14 +92,14 @@ xlabel('frequency (MHz) \rightarrow','FontSize',12);
 
 %% compare analytic and numerical wave-impedance %%%%%%%%%%%%%%%%%%%%%%%%%%
 figure
-ZL = uf1./if1;
 plot(freq*1e-6,real(ZL),'Linewidth',2);
 hold on;
 grid on;
 plot(freq*1e-6,imag(ZL),'r--','Linewidth',2);
-plot(freq*1e-6,ZL_a,'g-.','Linewidth',2);
+plot(freq*1e-6,port{1}.ZL,'g-.','Linewidth',2);
 ylabel('ZL (\Omega)','FontSize',12);
 xlabel('frequency (MHz) \rightarrow','FontSize',12);
 xlim([freq(1) freq(end)]*1e-6);
 l = legend('\Re(Z_L)','\Im(Z_L)','Z_L analytic','Location','Best');
 set(l,'FontSize',12);
+

matlab/Tutorials/Conical_Horn_Antenna.m

@@ -5,8 +5,8 @@
 % http://openems.de/index.php/Tutorial:_Conical_Horn_Antenna
 %
 % Tested with
-%  - Matlab 2011a / Octave 3.4.3
+%  - Matlab 2011a / Octave 3.6.3
-%  - openEMS v0.0.27
+%  - openEMS v0.0.31
 %
 % (C) 2011,2012 Thorsten Liebig <thorsten.liebig@uni-due.de>
 
@@ -40,35 +40,6 @@ f_stop  =  20e9;
 % frequency of interest
 f0 = 15e9;
 
-%% mode functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% by David M. Pozar, Microwave Engineering, third edition, page 113
-freq = linspace(f_start,f_stop,201);
-
-p11 = 1.841;
-kc = p11 / horn.radius /unit;
-k = 2*pi*freq/C0;
-fc = C0*kc/2/pi;
-beta = sqrt(k.^2 - kc^2);
-ZL_a = k * Z0 ./ beta;    %analytic waveguide impedance
-
-% mode profile E- and H-field
-kc = kc*unit;
-func_Er = [ num2str(-1/kc^2,'%14.13f') '/rho*cos(a)*j1('  num2str(kc,'%14.13f') '*rho)'];
-func_Ea = [ num2str(1/kc,'%14.13f') '*sin(a)*0.5*(j0('  num2str(kc,'%14.13f') '*rho)-jn(2,'  num2str(kc,'%14.13f') '*rho))'];
-func_Ex = ['(' func_Er '*cos(a) - ' func_Ea '*sin(a) ) * (rho<' num2str(horn.radius) ')'];
-func_Ey = ['(' func_Er '*sin(a) + ' func_Ea '*cos(a) ) * (rho<' num2str(horn.radius) ')'];
-
-func_Ha = [ num2str(-1/kc^2,'%14.13f') '/rho*cos(a)*j1('  num2str(kc,'%14.13f') '*rho)'];
-func_Hr = [ '-1*' num2str(1/kc,'%14.13f') '*sin(a)*0.5*(j0('  num2str(kc,'%14.13f') '*rho)-jn(2,'  num2str(kc,'%14.13f') '*rho))'];
-func_Hx = ['(' func_Hr '*cos(a) - ' func_Ha '*sin(a) ) * (rho<' num2str(horn.radius) ')'];
-func_Hy = ['(' func_Hr '*sin(a) + ' func_Ha '*cos(a) ) * (rho<' num2str(horn.radius) ')'];
-
-disp([' Cutoff frequencies for this mode and wavguide is: ' num2str(fc/1e9) ' GHz']);
-
-if (f_start<fc)
-    warning('openEMS:example','f_start is smaller than the cutoff-frequency, this may result in a long simulation... ');
-end
-
 %% setup FDTD parameter & excitation function
 FDTD = InitFDTD( 30000, 1e-4 );
 FDTD = SetGaussExcite(FDTD,0.5*(f_start+f_stop),0.5*(f_stop-f_start));
@@ -112,31 +83,17 @@ CSX = AddRotPoly(CSX,'Conical_Horn',10,0,2,p);
 % horn aperture
 A = pi*((horn.radius + sin(horn.angle)*horn.length)*unit)^2;
 
-% %% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% xy-mode profile excitation located directly on top of pml (first 8 z-lines)
+start=[-horn.radius -horn.radius mesh.z(10) ];
-CSX = AddExcitation(CSX,'excite',0,[1 1 0]);
+stop =[+horn.radius +horn.radius mesh.z(1)+horn.feed_length/2 ];
-weight{1} = func_Ex;
+[CSX, port] = AddCircWaveGuidePort( CSX, 0, 1, start, stop, horn.radius*unit, 'TE11', 0, 1);
-weight{2} = func_Ey;
-weight{3} = 0;
-CSX = SetExcitationWeight(CSX,'excite',weight);
-start=[0 0 mesh.z(8)-0.1 ];
-stop =[0 0 mesh.z(8)+0.1 ];
-CSX = AddCylinder(CSX,'excite',0 ,start,stop,horn.radius);
 
+%%
 CSX = AddDump(CSX,'Exc_dump');
-start=[-horn.radius -horn.radius mesh.z(8)-0.1 ];
+start=[-horn.radius -horn.radius mesh.z(8)];
-stop =[+horn.radius +horn.radius mesh.z(8)+0.1 ];
+stop =[+horn.radius +horn.radius mesh.z(8)];
 CSX = AddBox(CSX,'Exc_dump',0,start,stop);
 
-%% voltage and current definitions using the mode matching probes %%%%%%%%%
-%port 1
-start = [-horn.radius -horn.radius mesh.z(1)+horn.feed_length/2];
-stop  = [ horn.radius  horn.radius mesh.z(1)+horn.feed_length/2];
-CSX = AddProbe(CSX, 'ut1', 10, 1, [], 'ModeFunction',{func_Ex,func_Ey,0});
-CSX = AddBox(CSX,  'ut1',  0 ,start,stop);
-CSX = AddProbe(CSX,'it1', 11, 1, [], 'ModeFunction',{func_Hx,func_Hy,0});
-CSX = AddBox(CSX,'it1', 0 ,start,stop);
-
 %% nf2ff calc
 start = [mesh.x(9) mesh.y(9) mesh.z(9)];
 stop  = [mesh.x(end-8) mesh.y(end-8) mesh.z(end-8)];
@@ -159,16 +116,17 @@ CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
 RunOpenEMS( Sim_Path, Sim_CSX);
 
 %% postprocessing & do the plots
-U = ReadUI( 'ut1', Sim_Path, freq ); % time domain/freq domain voltage
+freq = linspace(f_start,f_stop,201);
-I = ReadUI( 'it1', Sim_Path, freq ); % time domain/freq domain current (half time step is corrected)
+
+port = calcPort(port, Sim_Path, freq);
+
+Zin = port.uf.tot ./ port.if.tot;
+s11 = port.uf.ref ./ port.uf.inc;
+P_in = 0.5 * port.uf.inc .* conj( port.if.inc ); % antenna feed power
+
 
 % plot reflection coefficient S11
 figure
-uf_inc = 0.5*(U.FD{1}.val + I.FD{1}.val .* ZL_a);
-if_inc = 0.5*(I.FD{1}.val + U.FD{1}.val ./ ZL_a);
-uf_ref = U.FD{1}.val - uf_inc;
-if_ref = if_inc - I.FD{1}.val;
-s11 = uf_ref ./ uf_inc;
 plot( freq/1e9, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 );
 ylim([-60 0]);
 grid on
@@ -176,8 +134,6 @@ title( 'reflection coefficient S_{11}' );
 xlabel( 'frequency f / GHz' );
 ylabel( 'reflection coefficient |S_{11}|' );
 
-P_in = 0.5*uf_inc .* conj( if_inc ); % antenna feed power
-
 drawnow
 
 %% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matlab/Tutorials/Horn_Antenna.m

@@ -43,37 +43,12 @@ f_stop  =  20e9;
 f0 = 15e9;
 
 %waveguide TE-mode definition
-m = 1;
+TE_mode = 'TE10';
-n = 0;
-
-%% mode functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% by David M. Pozar, Microwave Engineering, third edition, page 113
-freq = linspace(f_start,f_stop,201);
 a = horn.width;
 b = horn.height;
-k = 2*pi*freq/c0;
-kc = sqrt((m*pi/a/unit)^2 + (n*pi/b/unit)^2);
-fc = c0*kc/2/pi;          %cut-off frequency
-beta = sqrt(k.^2 - kc^2); %waveguide phase-constant
-ZL_a = k * Z0 ./ beta;    %analytic waveguide impedance
-
-% mode profile E- and H-field
-x_pos = ['(x-' num2str(a/2) ')'];
-y_pos = ['(y-' num2str(b/2) ')'];
-func_Ex = [num2str( n/b/unit) '*cos(' num2str(m*pi/a) '*' x_pos ')*sin('  num2str(n*pi/b) '*' y_pos ')'];
-func_Ey = [num2str(-m/a/unit) '*sin(' num2str(m*pi/a) '*' x_pos ')*cos('  num2str(n*pi/b) '*' y_pos ')'];
-
-func_Hx = [num2str(m/a/unit) '*sin(' num2str(m*pi/a) '*' x_pos ')*cos('  num2str(n*pi/b) '*' y_pos ')'];
-func_Hy = [num2str(n/b/unit) '*cos(' num2str(m*pi/a) '*' x_pos ')*sin('  num2str(n*pi/b) '*' y_pos ')'];
-
-disp([' Cutoff frequencies for this mode and wavguide is: ' num2str(fc/1e9) ' GHz']);
-
-if (f_start<fc)
-    warning('openEMS:example','f_start is smaller than the cutoff-frequency, this may result in a long simulation... ');
-end
 
 %% setup FDTD parameter & excitation function
-FDTD = InitFDTD( 30000, 1e-4 );
+FDTD = InitFDTD('EndCriteria', 1e-4);
 FDTD = SetGaussExcite(FDTD,0.5*(f_start+f_stop),0.5*(f_stop-f_start));
 BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'}; % boundary conditions
 FDTD = SetBoundaryCond( FDTD, BC );
@@ -138,25 +113,10 @@ CSX = AddLinPoly( CSX, 'horn', 10, 0, -horn.thickness/2, p, horn.thickness, 'Tra
 % horn aperture
 A = (a + 2*sin(horn.angle(1))*horn.length)*unit * (b + 2*sin(horn.angle(2))*horn.length)*unit;
 
-% %% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% xy-mode profile excitation located directly on top of pml (first 8 z-lines)
-CSX = AddExcitation(CSX,'excite',0,[1 1 0]);
-weight{1} = func_Ex;
-weight{2} = func_Ey;
-weight{3} = 0;
-CSX = SetExcitationWeight(CSX,'excite',weight);
 start=[-a/2 -b/2 mesh.z(8) ];
-stop =[ a/2  b/2 mesh.z(8) ];
+stop =[ a/2  b/2 mesh.z(1)+horn.feed_length/2 ];
-CSX = AddBox(CSX,'excite',0 ,start,stop);
+[CSX, port] = AddRectWaveGuidePort( CSX, 0, 1, start, stop, 2, a*unit, b*unit, TE_mode, 1);
-
-%% voltage and current definitions using the mode matching probes %%%%%%%%%
-%port 1
-start(3) = mesh.z(1)+horn.feed_length/2;
-stop(3)  = start(3);
-CSX = AddProbe(CSX, 'ut1', 10, 1, [], 'ModeFunction',{func_Ex,func_Ey,0});
-CSX = AddBox(CSX,  'ut1',  0 ,start,stop);
-CSX = AddProbe(CSX,'it1', 11, 1, [], 'ModeFunction',{func_Hx,func_Hy,0});
-CSX = AddBox(CSX,'it1', 0 ,start,stop);
 
 %% nf2ff calc
 start = [mesh.x(9) mesh.y(9) mesh.z(9)];
@@ -180,16 +140,14 @@ CSXGeomPlot([Sim_Path '/' Sim_CSX]);
 RunOpenEMS(Sim_Path, Sim_CSX);
 
 %% postprocessing & do the plots
-U = ReadUI( 'ut1', Sim_Path, freq ); % time domain/freq domain voltage
+freq = linspace(f_start,f_stop,201);
-I = ReadUI( 'it1', Sim_Path, freq ); % time domain/freq domain current (half time step is corrected)
+
+port = calcPort(port, Sim_Path, freq);
+
+Zin = port.uf.tot ./ port.if.tot;
+s11 = port.uf.ref ./ port.uf.inc;
+P_in = 0.5 * port.uf.inc .* conj( port.if.inc ); % antenna feed power
 
-% plot reflection coefficient S11
-figure
-uf_inc = 0.5*(U.FD{1}.val + I.FD{1}.val .* ZL_a);
-if_inc = 0.5*(I.FD{1}.val + U.FD{1}.val ./ ZL_a);
-uf_ref = U.FD{1}.val - uf_inc;
-if_ref = if_inc - I.FD{1}.val;
-s11 = uf_ref ./ uf_inc;
 plot( freq/1e9, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 );
 ylim([-60 0]);
 grid on
@@ -197,8 +155,6 @@ title( 'reflection coefficient S_{11}' );
 xlabel( 'frequency f / GHz' );
 ylabel( 'reflection coefficient |S_{11}|' );
 
-P_in = 0.5*uf_inc .* conj( if_inc ); % antenna feed power
-
 drawnow
 
 %% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

matlab/Tutorials/Rect_Waveguide.m

@@ -6,9 +6,9 @@
 %
 % Tested with
 %  - Matlab 2011a / Octave 3.4.3
-%  - openEMS v0.0.26
+%  - openEMS v0.0.31
 %
-% (C) 2010-2012 Thorsten Liebig <thorsten.liebig@gmx.de>
+% (C) 2010-2013 Thorsten Liebig <thorsten.liebig@gmx.de>
 
 close all
 clear
@@ -17,7 +17,6 @@ clc
 %% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 physical_constants;
 unit = 1e-3; %drawing unit in mm
-numTS = 50000; %max. number of timesteps
 
 % waveguide dimensions
 length = 5000;
@@ -29,35 +28,12 @@ f_start =  300e6;
 f_stop  =  500e6;
 
 %waveguide TE-mode definition
-m = 1;
+TE_mode = 'TE11';
-n = 1;
 
 mesh_res = [10 10 10]; %targeted mesh resolution
 
-%% mode functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% by David M. Pozar, Microwave Engineering, third edition, page 113
-freq = linspace(f_start,f_stop,201);
-k = 2*pi*freq/c0;
-kc = sqrt((m*pi/a/unit)^2 + (n*pi/b/unit)^2);
-fc = c0*kc/2/pi;          %cut-off frequency
-beta = sqrt(k.^2 - kc^2); %waveguide phase-constant
-ZL_a = k * Z0 ./ beta;    %analytic waveguide impedance
-
-% mode profile E- and H-field
-func_Ex = [num2str( n/b/unit) '*cos(' num2str(m*pi/a) '*x)*sin('  num2str(n*pi/b) '*y)'];
-func_Ey = [num2str(-m/a/unit) '*sin(' num2str(m*pi/a) '*x)*cos('  num2str(n*pi/b) '*y)'];
-
-func_Hx = [num2str(m/a/unit) '*sin(' num2str(m*pi/a) '*x)*cos('  num2str(n*pi/b) '*y)'];
-func_Hy = [num2str(n/b/unit) '*cos(' num2str(m*pi/a) '*x)*sin('  num2str(n*pi/b) '*y)'];
-
-disp([' Cutoff frequencies for this mode and wavguide is: ' num2str(fc/1e6) ' MHz']);
-
-if (f_start<fc)
-    warning('openEMS:example','f_start is smaller than the cutoff-frequency, this may result in a long simulation... ');
-end
-
 %% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-FDTD = InitFDTD(numTS,1e-5);
+FDTD = InitFDTD();
 FDTD = SetGaussExcite(FDTD,0.5*(f_start+f_stop),0.5*(f_stop-f_start));
 
 % boundary conditions
@@ -71,35 +47,14 @@ mesh.y = SmoothMeshLines([0 b], mesh_res(2));
 mesh.z = SmoothMeshLines([0 length], mesh_res(3));
 CSX = DefineRectGrid(CSX, unit,mesh);
 
-%% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% apply the waveguide port %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% xy-mode profile excitation located directly on top of pml (first 8 z-lines)
+start=[mesh.x(1)   mesh.y(1)   mesh.z(8)];
-CSX = AddExcitation(CSX,'excite',0,[1 1 0]);
+stop =[mesh.x(end) mesh.y(end) mesh.z(15)];
-weight{1} = func_Ex;
+[CSX, port{1}] = AddRectWaveGuidePort( CSX, 0, 1, start, stop, 2, a*unit, b*unit, TE_mode, 1);
-weight{2} = func_Ey;
-weight{3} = 0;
-CSX = SetExcitationWeight(CSX,'excite',weight);
-start=[mesh.x(1)   mesh.y(1)   mesh.z(8) ];
-stop =[mesh.x(end) mesh.y(end) mesh.z(8) ];
-CSX = AddBox(CSX,'excite',0 ,start,stop);
 
-%% voltage and current definitions using the mode matching probes %%%%%%%%%
+start=[mesh.x(1)   mesh.y(1)   mesh.z(end-13)];
-%port 1
+stop =[mesh.x(end) mesh.y(end) mesh.z(end-14)];
-start = [mesh.x(1)   mesh.y(1)   mesh.z(15)];
+[CSX, port{2}] = AddRectWaveGuidePort( CSX, 0, 2, start, stop, 2, a*unit, b*unit, TE_mode);
-stop  = [mesh.x(end) mesh.y(end) mesh.z(15)];
-CSX = AddProbe(CSX, 'ut1', 10, 'ModeFunction',{func_Ex,func_Ey,0});
-CSX = AddBox(CSX,  'ut1',  0 ,start,stop);
-CSX = AddProbe(CSX,'it1', 11, 'ModeFunction',{func_Hx,func_Hy,0});
-CSX = AddBox(CSX,'it1', 0 ,start,stop);
-
-%port 2
-start = [mesh.x(1)   mesh.y(1)   mesh.z(end-14)];
-stop  = [mesh.x(end) mesh.y(end) mesh.z(end-14)];
-CSX = AddProbe(CSX, 'ut2', 10, 'ModeFunction',{func_Ex,func_Ey,0});
-CSX = AddBox(CSX,  'ut2',  0 ,start,stop);
-CSX = AddProbe(CSX,'it2', 11, 'ModeFunction',{func_Hx,func_Hy,0});
-CSX = AddBox(CSX,'it2', 0 ,start,stop);
-
-port_dist = mesh.z(end-14) - mesh.z(15);
 
 %% define dump box... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 CSX = AddDump(CSX,'Et','FileType',1,'SubSampling','4,4,4');
@@ -108,7 +63,7 @@ stop  = [mesh.x(end) mesh.y(end) mesh.z(end)];
 CSX = AddBox(CSX,'Et',0 , start,stop);
 
 %% Write openEMS compatoble xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Sim_Path = 'tmp';
+Sim_Path = 'tmp_mod';
 Sim_CSX = 'rect_wg.xml';
 
 [status, message, messageid] = rmdir(Sim_Path,'s');
@@ -119,29 +74,16 @@ WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
 RunOpenEMS(Sim_Path, Sim_CSX)
 
 %% postproc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-U = ReadUI({'ut1','ut2'},[Sim_Path '/'],freq);
+freq = linspace(f_start,f_stop,201);
-I = ReadUI({'it1','it2'},[Sim_Path '/'],freq);
+port = calcPort( port, Sim_Path, freq);
-Exc = ReadUI('et',Sim_Path,freq);
 
-uf1 = U.FD{1}.val./Exc.FD{1}.val;
+s11 = port{1}.uf.ref./ port{1}.uf.inc;
-uf2 = U.FD{2}.val./Exc.FD{1}.val;
+s21 = port{2}.uf.ref./ port{1}.uf.inc;
-if1 = I.FD{1}.val./Exc.FD{1}.val;
+ZL = port{1}.uf.tot./port{1}.if.tot;
-if2 = I.FD{2}.val./Exc.FD{1}.val;
+ZL_a = port{1}.ZL; % analytic waveguide impedance
-
-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_inc - if1; 
-uf2_ref = uf2 - uf2_inc;
-if2_ref = if2_inc - if2;
 
 %% plot s-parameter %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 figure
-s11 = uf1_ref./uf1_inc;
-s21 = uf2_inc./uf1_inc;
 plot(freq*1e-6,20*log10(abs(s11)),'k-','Linewidth',2);
 xlim([freq(1) freq(end)]*1e-6);
 grid on;
@@ -154,7 +96,6 @@ xlabel('frequency (MHz) \rightarrow','FontSize',12);
 
 %% compare analytic and numerical wave-impedance %%%%%%%%%%%%%%%%%%%%%%%%%%
 figure
-ZL = uf1./if1;
 plot(freq*1e-6,real(ZL),'Linewidth',2);
 hold on;
 grid on;