diff --git a/matlab/AddCircWaveGuidePort.m b/matlab/AddCircWaveGuidePort.m new file mode 100644 index 0000000..4993489 --- /dev/null +++ b/matlab/AddCircWaveGuidePort.m @@ -0,0 +1,109 @@ +function [CSX,port] = AddCircWaveGuidePort( CSX, prio, portnr, start, stop, radius, mode_name, pol_ang, exc_amp, varargin ) +% function [CSX,port] = AddCircWaveGuidePort( CSX, prio, portnr, start, stop, radius, mode_name, pol_ang, exc_amp, varargin ) +% +% Create a circular waveguide port, including an optional excitation and probes +% +% Note: - The excitation will be located at the start position in the given direction +% - The voltage and current probes at the stop position in the given direction +% +% input: +% CSX: complete CSX structure (must contain a mesh) +% prio: priority of primitives +% start: start coordinates of waveguide port box +% stop: stop coordinates of waveguide port box +% radius: circular waveguide radius (in meter) +% mode_name: mode name, e.g. 'TE11' or 'TM21' +% pol_ang: polarization angle (e.g. 0 = horizontal, pi/2 = vertical) +% exc_amp: excitation amplitude (set 0 to be passive) +% varargin: optional additional excitations options, see also AddExcitation +% +% output: +% CSX: modified CSX structure +% port: port structure to use with calcPort +% +% example: +% % create a TE11 circular waveguide mode, using cylindircal coordinates +% start=[mesh.r(1) mesh.a(1) 0 ]; +% stop =[mesh.r(end) mesh.a(end) 100]; +% [CSX,port] = AddCircWaveGuidePort( CSX, 99, 1, start, stop, 320e-3, 'TE11', 0, 1); +% +% openEMS matlab interface +% ----------------------- +% (c) 2013 Thorsten Liebig (thorsten.liebig@gmx.de) +% +% See also InitCSX, AddExcitation, calcWGPort, calcPort + +if (~strcmpi(mode_name(1:2),'TE')) + error 'currently only TE type modes are supported' +end + +if (nargin<9) + exc_amp = 0; +end +if (nargin<8) + pol_ang = 0; +end + +pnm = 0; +n = str2double(mode_name(3)); +m = str2double(mode_name(4)); + +% values by David M. Pozar, Microwave Engineering, third edition +if ((n==0) && (m==1)) + pnm = 3.832; +elseif ((n==1) && (m==1)) + pnm = 1.841; +elseif ((n==2) && (m==1)) + pnm = 3.054; +elseif ((n==0) && (m==2)) + pnm = 7.016; +elseif ((n==1) && (m==2)) + pnm = 5.331; +elseif ((n==2) && (m==2)) + pnm = 6.706; +elseif ((n==0) && (m==3)) + pnm = 10.174; +elseif ((n==1) && (m==3)) + pnm = 8.536; +elseif ((n==2) && (m==3)) + pnm = 9.970; +else + error 'invalid TE_nm mode' +end + +if ~isfield(CSX,'RectilinearGrid') + error 'mesh needs to be defined! Use DefineRectGrid() first!'; + if (~isfield(CSX.RectilinearGrid,'XLines') || ~isfield(CSX.RectilinearGrid,'YLines') || ~isfield(CSX.RectilinearGrid,'ZLines')) + error 'mesh needs to be defined! Use DefineRectGrid() first!'; + end +end + +unit = CSX.RectilinearGrid.ATTRIBUTE.DeltaUnit; +kc = pnm/radius; +kc_draw = kc*unit; + +angle = ['a-' num2str(pol_ang)]; +% functions by David M. Pozar, Microwave Engineering, third edition +% electric field mode profile +func_Er = [ num2str(-1/kc_draw^2,15) '/rho*cos(' angle ')*j1(' num2str(kc_draw,15) '*rho)']; +func_Ea = [ num2str(1/kc_draw,15) '*sin(' angle ')*0.5*(j0(' num2str(kc_draw,15) '*rho)-jn(2,' num2str(kc_draw,15) '*rho))']; + +% magnetic field mode profile +func_Hr = [ num2str(-1/kc_draw,15) '*sin(' angle ')*0.5*(j0(' num2str(kc_draw,15) '*rho)-jn(2,' num2str(kc_draw,15) '*rho))']; +func_Ha = [ num2str(-1/kc_draw^2,15) '/rho*cos(' angle ')*j1(' num2str(kc_draw,15) '*rho)']; + +if (CSX.ATTRIBUTE.CoordSystem==1) + func_E = {func_Er, func_Ea, 0}; + func_H = {func_Hr, func_Ha, 0}; +else + func_Ex = ['(' func_Er '*cos(a) - ' func_Ea '*sin(a) ) * (rho<' num2str(radius/unit) ')']; + func_Ey = ['(' func_Er '*sin(a) + ' func_Ea '*cos(a) ) * (rho<' num2str(radius/unit) ')']; + func_E = {func_Ex, func_Ey, 0}; + + func_Hx = ['(' func_Hr '*cos(a) - ' func_Ha '*sin(a) ) * (rho<' num2str(radius/unit) ')']; + func_Hy = ['(' func_Hr '*sin(a) + ' func_Ha '*cos(a) ) * (rho<' num2str(radius/unit) ')']; + func_H = {func_Hx, func_Hy, 0}; +end + +[CSX,port] = AddWaveGuidePort( CSX, prio, portnr, start, stop, 2, func_E, func_H, kc, exc_amp, varargin{:} ); + diff --git a/matlab/AddRectWaveGuidePort.m b/matlab/AddRectWaveGuidePort.m new file mode 100644 index 0000000..83f87bb --- /dev/null +++ b/matlab/AddRectWaveGuidePort.m @@ -0,0 +1,89 @@ +function [CSX,port] = AddRectWaveGuidePort( CSX, prio, portnr, start, stop, dir, a, b, mode_name, exc_amp, varargin ) +% function [CSX,port] = AddRectWaveGuidePort( CSX, prio, portnr, start, stop, dir, a, b, mode_name, exc_amp, varargin ) +% +% Create a rectangular waveguide port, including an optional excitation and probes +% +% Note: - The excitation will be located at the start position in the given direction +% - The voltage and current probes at the stop position in the given direction +% +% input: +% CSX: complete CSX structure (must contain a mesh) +% prio: priority of primitives +% start: start coordinates of waveguide port box +% stop: stop coordinates of waveguide port box +% dir: direction of port (0/1/2 for x/y/z-direction) +% a,b: rectangular waveguide width and height (in meter) +% mode_name: mode name, e.g. 'TE11' or 'TM21' +% exc_amp: excitation amplitude (set 0 to be passive) +% varargin: optional additional excitations options, see also AddExcitation +% +% output: +% CSX: modified CSX structure +% port: port structure to use with calcPort +% +% example: +% % create a TE10 circular waveguide mode, using cylindircal coordinates +% start=[mesh.r(1) mesh.a(1) 0 ]; +% stop =[mesh.r(end) mesh.a(end) 100]; +% [CSX,port] = AddCircWaveGuidePort( CSX, 99, 1, start, stop, 320e-3, 'TE11', 0, 1); +% +% openEMS matlab interface +% ----------------------- +% (c) 2013 Thorsten Liebig (thorsten.liebig@gmx.de) +% +% See also InitCSX, AddExcitation, calcWGPort, calcPort + +if (~strcmpi(mode_name(1:2),'TE')) + error 'currently only TE type modes are supported' +end + +if (nargin<10) + exc_amp = 0; +end + +m = str2double(mode_name(3)); +n = str2double(mode_name(4)); + +% values by David M. Pozar, Microwave Engineering, third edition +kc = sqrt((m*pi/a)^2 + (n*pi/b)^2); + +if ~isfield(CSX,'RectilinearGrid') + error 'mesh needs to be defined! Use DefineRectGrid() first!'; + if (~isfield(CSX.RectilinearGrid,'XLines') || ~isfield(CSX.RectilinearGrid,'YLines') || ~isfield(CSX.RectilinearGrid,'ZLines')) + error 'mesh needs to be defined! Use DefineRectGrid() first!'; + end +end + +unit = CSX.RectilinearGrid.ATTRIBUTE.DeltaUnit; +kc_draw = kc*unit; + +dir_names={'x','y','z'}; + +dirP = mod((dir+1),3)+1; +dirPP = mod((dir+2),3)+1; +nameX = ['(' dir_names{dirP} '-' num2str(start(dirP)) ')']; +nameY = ['(' dir_names{dirPP} '-' num2str(start(dirPP)) ')']; + +%convert a&b to drawing units +a = a/unit; +b = b/unit; +% functions by David M. Pozar, Microwave Engineering, third edition +% electric field mode profile +func_Ex = [num2str( n/b) '*cos(' num2str(m*pi/a) '*' nameX ')*sin(' num2str(n*pi/b) '*' nameY ')']; +func_Ey = [num2str(-m/a) '*sin(' num2str(m*pi/a) '*' nameX ')*cos(' num2str(n*pi/b) '*' nameY ')']; + +% magnetic field mode profile +func_Hx = [num2str(m/a) '*sin(' num2str(m*pi/a) '*' nameX ')*cos(' num2str(n*pi/b) '*' nameY ')']; +func_Hy = [num2str(n/b) '*cos(' num2str(m*pi/a) '*' nameX ')*sin(' num2str(n*pi/b) '*' nameY ')']; + + +func_E{dir+1} = 0; +func_E{dirP} = func_Ex; +func_E{dirPP} = func_Ey; + +func_H{dir+1} = 0; +func_H{dirP} = func_Hx; +func_H{dirPP} = func_Hy; + +[CSX,port] = AddWaveGuidePort( CSX, prio, portnr, start, stop, dir, func_E, func_H, kc, exc_amp, varargin{:} ); + diff --git a/matlab/AddWaveGuidePort.m b/matlab/AddWaveGuidePort.m new file mode 100644 index 0000000..cc012c2 --- /dev/null +++ b/matlab/AddWaveGuidePort.m @@ -0,0 +1,111 @@ +function [CSX,port] = AddWaveGuidePort( CSX, prio, portnr, start, stop, dir, E_WG_func, H_WG_func, kc, exc_amp, varargin ) +% function [CSX,port] = AddWaveGuidePort( CSX, prio, portnr, start, stop, dir, E_WG_func, H_WG_func, kc, exc_amp, varargin ) +% +% Create a waveguide port, including an optional excitation and probes +% +% Note: - The excitation will be located at the start position in the given direction +% - The voltage and current probes at the stop position in the given direction +% +% input: +% CSX: complete CSX structure (must contain a mesh) +% prio: priority of primitives +% start: start coordinates of waveguide port box +% stop: stop coordinates of waveguide port box +% dir: direction of port (0/1/2 for x/y/z-direction) +% E_WG_func: electric field mode profile function as a string +% H_WG_func: magnetic field mode profile function as a string +% kc: cutoff wavenumber (defined by the waveguide dimensions) +% exc_amp: excitation amplitude (set 0 to be passive) +% varargin: optional additional excitations options, see also AddExcitation +% +% output: +% CSX: modified CSX structure +% port: port structure to use with calcPort +% +% example: +% % create a TE11 circular waveguide mode, using cylindircal coordinates +% p11 = 1.841; +% kc = p11 / radius; % cutoff wavenumber with radius in meter +% kc_draw = kc*unit; % cutoff wavenumber in drawing units +% +% % electric field mode profile +% func_E{1} = [ num2str(-1/kc_draw^2,15) '/rho*cos(a)*j1(' num2str(kc_draw,15) '*rho)']; +% func_E{2} = [ num2str(1/kc_draw,15) '*sin(a)*0.5*(j0(' num2str(kc_draw,15) '*rho)-jn(2,' num2str(kc_draw,15) '*rho))']; +% func_E{3} = 0; +% +% % magnetic field mode profile +% func_H{1} = [ '-1*' num2str(1/kc_draw,15) '*sin(a)*0.5*(j0(' num2str(kc_draw,15) '*rho)-jn(2,' num2str(kc_draw,15) '*rho))']; +% func_H{2} = [ num2str(-1/kc_draw^2,15) '/rho*cos(a)*j1(' num2str(kc_draw,15) '*rho)']; +% func_H{3} = 0; +% +% start=[mesh.r(1) mesh.a(1) 0 ]; +% stop =[mesh.r(end) mesh.a(end) 100]; +% [CSX, port{1}] = AddWaveGuidePort(CSX, 0, 1, start, stop, 2, func_E, func_H, kc, 1); +% +% openEMS matlab interface +% ----------------------- +% (c) 2013 Thorsten Liebig (thorsten.liebig@gmx.de) +% +% See also InitCSX, AddExcitation, calcWGPort, calcPort + +%check mesh +if ~isfield(CSX,'RectilinearGrid') + error 'mesh needs to be defined! Use DefineRectGrid() first!'; + if (~isfield(CSX.RectilinearGrid,'XLines') || ~isfield(CSX.RectilinearGrid,'YLines') || ~isfield(CSX.RectilinearGrid,'ZLines')) + error 'mesh needs to be defined! Use DefineRectGrid() first!'; + end +end + +port.type='WaveGuide'; +port.nr=portnr; +port.kc = kc; +port.dir = dir; +port.drawingunit = CSX.RectilinearGrid.ATTRIBUTE.DeltaUnit; + +if ~( (dir==0) || (dir==1) || (dir==2) ) + error 'dir must be 0, 1 or 2' +end + +% matlab adressing +dir = dir + 1; +dir_sign = sign(stop(dir) - start(dir)); +if (dir_sign==0) + dir_sign = 1; +end + +port.direction = dir_sign; + +E_WG_func{dir} = 0; +H_WG_func{dir} = 0; + +port.excite = 0; +if (exc_amp~=0) + if (start(dir)==stop(dir)) + error 'if waveguide port is to be excited, the length in propagation direction must not be zero' + end + e_start = start; + e_stop = stop; + e_stop(dir) = e_start(dir); + port.excite = 1; + port.excitepos = e_start(dir); + e_vec = [1 1 1]*exc_amp; + e_vec(dir) = 0; + exc_name = ['port_excite_' num2str(portnr)]; + CSX = AddExcitation( CSX, exc_name, 0, e_vec, varargin{:}); + CSX = SetExcitationWeight(CSX, exc_name, E_WG_func ); + CSX = AddBox( CSX, exc_name, prio, e_start, e_stop); +end + +% voltage/current planes +m_start = start; +m_stop = stop; +m_start(dir) = stop(dir); + +port.measplanepos = m_start(dir); +port.U_filename = ['port_ut' int2str(portnr)]; +CSX = AddProbe(CSX, port.U_filename, 10, 'ModeFunction', E_WG_func); +CSX = AddBox(CSX, port.U_filename, 0 ,m_start, m_stop); + +port.I_filename = ['port_it' int2str(portnr)]; +CSX = AddProbe(CSX, port.I_filename, 11, 'ModeFunction', H_WG_func, 'weight', dir_sign); +CSX = AddBox(CSX, port.I_filename, 0 ,m_start, m_stop); diff --git a/matlab/calcPort.m b/matlab/calcPort.m index 98d90f4..418e095 100644 --- a/matlab/calcPort.m +++ b/matlab/calcPort.m @@ -48,6 +48,9 @@ end if strcmpi(port.type,'MSL') port = calcTLPort( port, SimDir, f, varargin{:}); return +elseif strcmpi(port.type,'WaveGuide') + port = calcWGPort( port, SimDir, f, varargin{:}); + return elseif (strcmpi(port.type,'Lumped') || strcmpi(port.type,'Curve')) port = calcLumpedPort( port, SimDir, f, varargin{:}); return diff --git a/matlab/calcWGPort.m b/matlab/calcWGPort.m new file mode 100644 index 0000000..05fca16 --- /dev/null +++ b/matlab/calcWGPort.m @@ -0,0 +1,121 @@ +function [port] = calcWGPort( port, SimDir, f, varargin) +% [port] = calcTLPort( port, SimDir, f, varargin) +% +% Calculate voltages and currents, the propagation constant beta +% and the characteristic impedance ZL of the given waveguide port. +% +% The port has to be created by e.g. AddWaveGuidePort(). +% +% input: +% port: return value of e.g. AddWaveGuidePort() +% SimDir: directory, where the simulation files are +% f: frequency vector for DFT +% +% variable input: +% 'RefImpedance': - use a given reference impedance to calculate inc and +% ref voltages and currents +% - default is given port or calculated line impedance +% 'RefPlaneShift': - use a given reference plane shift from port beginning +% for a desired phase correction +% - default is the measurement plane at the end of the +% port +% - the plane shift has to be given in drawing units! +% +% output: +% port.f the given frequency fector +% port.uf.tot/inc/ref total, incoming and reflected voltage +% port.if.tot/inc/ref total, incoming and reflected current +% port.beta: propagation constant +% port.ZL: characteristic line impedance +% port.ZL_ref used reference impedance +% +% example: +% port{1} = calcWGPort( port{1}, Sim_Path, f, 'RefImpedance', 50); +% +% openEMS matlab interface +% ----------------------- +% (C) 2013 Thorsten Liebig (thorsten.liebig@gmx.de) +% +% See also AddWaveGuidePort, calcPort + +if (iscell(port)) + for n=1:numel(port) + port{n}=calcTLPort(port{n}, SimDir, f, varargin{:}); + end + return; +end + +if (strcmpi(port.type,'WaveGuide')~=1) + error('openEMS:calcWGPort','error, type is not a waveguide port'); +end + +%set defaults +ref_ZL = -1; +ref_shift = nan; + +UI_args = {}; + +for n=1:2:numel(varargin) + if (strcmp(varargin{n},'RefPlaneShift')==1); + ref_shift = varargin{n+1}; + elseif (strcmp(varargin{n},'RefImpedance')==1); + ref_ZL = varargin{n+1}; + else + UI_args(end+1) = varargin(n); + UI_args(end+1) = varargin(n+1); + end +end + +% read time domain data +U = ReadUI( port.U_filename, SimDir, f, UI_args{:} ); +I = ReadUI( port.I_filename, SimDir, f, UI_args{:} ); + +% store the original frequency domain waveforms +u_f = U.FD{1}.val; +i_f = I.FD{1}.val; + +physical_constants +k = 2*pi*f/C0; +fc = C0*port.kc/2/pi; +port.beta = sqrt(k.^2 - port.kc^2); +port.ZL = k * Z0 ./ port.beta; %analytic waveguide impedance + + +% reference plane shift (lossless) +if ~isnan(ref_shift) + % shift relative to the beginning of the waveguide + ref_shift = ref_shift - port.measplanepos; + ref_shift = ref_shift * port.drawingunit; + + % store the shifted frequency domain waveforms + phase = real(beta)*ref_shift; + u_f_shift = u_f .* cos(-phase) + 1i * i_f.*port.ZL .* sin(-phase); + i_f_shift = i_f .* cos(-phase) + 1i * u_f./port.ZL .* sin(-phase); + + u_f = u_f_shift; + i_f = i_f_shift; +end + +if (ref_ZL < 0) + ref_ZL = port.ZL; +end + +port.ZL_ref = ref_ZL; + +port.f = f; +uf_inc = 0.5 * ( u_f + i_f .* ref_ZL ); +if_inc = 0.5 * ( i_f + u_f ./ ref_ZL ); + +uf_ref = u_f - uf_inc; +if_ref = if_inc - i_f; + +port.uf.tot = u_f; +port.uf.inc = uf_inc; +port.uf.ref = uf_ref; + +port.if.tot = i_f; +port.if.inc = if_inc; +port.if.ref = if_ref; + +port.raw.U = U; +port.raw.I = I;