matlab: revision of AddMSLPort + calcPort

all examples using this functions need to be revised!

matlab/AddMSLPort.m

@@ -1,38 +1,58 @@
-function [CSX,port] = AddMSLPort( CSX, prio, portnr, materialname, start, stop, dir, evec, refplaneshift, excitename )
+function [CSX,port] = AddMSLPort( CSX, prio, portnr, materialname, start, stop, dir, evec, varargin )
-% [CSX,port] = AddMSLPort( CSX, prio, portnr, materialname, start, stop, dir, evec, refplaneshift, excitename )
+% [CSX,port] = AddMSLPort( CSX, prio, portnr, materialname, start, stop, dir, evec, varargin )
 %
 % CSX:          CSX-object created by InitCSX()
 % prio:         priority for excitation and probe boxes
 % portnr:       (integer) number of the port
-% materialname: property for the MSL (created by AddMetal() or AddMaterial())
+% materialname: property for the MSL (created by AddMetal())
 % start:        3D start rowvector for port definition
 % stop:         3D end rowvector for port definition
-% dir:          direction of wave propagation (choices: [1 0 0], [0 1 0] or [0 0 1])
+% dir:          direction of wave propagation (choices: 0 1 2)
 % evec:         excitation vector, which defines the direction of the e-field (must be the same as used in AddExcitation())
-% refplaneshift (optional): if not specified or empty, the measurement
+% 
-%    plane is used; if specified, reference plane is shifted by
+% variable input:
-%    <refplaneshift> starting from <start> (thus refplaneshift is normally
+% 'ExcitePort'  necessary excitation name to make the port an active feeding port
-%    positive)
+% 'FeedShift'   shift to port from start by a given distance in drawing
-% excitename (optional): if specified, the port will be switched on (see AddExcitation())
+%               units. Default is 0. Only active if 'ExcitePort' is set!
+% 'Feed_R'      Specifiy a lumped port resistance. Default is no lumped
+%               port resistance --> port has to end in an ABC.
+%               Only active if 'ExcitePort' is set!
+% 'MeasPlaneShift'  Shift the measurement plane from start t a given distance 
+%               in drawing units. Default is the middle of start/stop.
+%               Only active if 'ExcitePort' is set!
 % 
 % the mesh must be already initialized
 %
 % example:
-%   start = [0 0 height]; stop = [length width 0]; dir = [1 0 0]; evec = [0 0 -1]
+%   start = [0 0 height]; 
-%   this defines a MSL in x-direction (dir) with an e-field excitation in -z-direction (evec)
+%   stop = [length width 0]; 
-%   the excitation is placed at x=start(1); the wave travels towards x=stop(1)
+%   CSX = AddMetal( CSX, 'metal' ); %create a PEC called 'metal'
-%   the MSL-metal is created in xy-plane at z=start(3)
+%   [CSX,port] = AddMSLPort( CSX, 0, 1, 'metal', start, stop,  ...
+%                0, [0 0 -1] , 'ExcitePort', 'excite', 'Feed_R', 50 )
 % 
-% Sebastian Held <sebastian.held@gmx.de>
+%   this defines a MSL in x-direction (dir=0) with an e-field excitation 
-% May 13 2010
+%   in -z-direction (evec=[0 0 -1]) the excitation is placed at x=start(1); 
+%   the wave travels towards x=stop(1) the MSL-metal is created 
+%   in xy-plane at z=start(3)
 %
-% See also InitCSX AddMetal AddMaterial AddExcitation calcMSLPort
+% Sebastian Held <sebastian.held@gmx.de> May 13 2010
+% Thorsten Liebig <thorsten.liebig@gmx.de> Sept 16 2011
+%
+% See also InitCSX AddMetal AddMaterial AddExcitation calcPort
+
+%% validate arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%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
 
 % check dir
-if ~(dir(1) == dir(2) == 0) && ~(dir(1) == dir(3) == 0) && ~(dir(2) == dir(3) == 0) || (sum(dir) == 0)
+if ~( (dir >= 0) && (dir <= 2) )
 	error 'dir must have exactly one component ~= 0'
 end
-dir = dir ./ sum(dir); % dir is now a unit vector
 
 % check evec
 if ~(evec(1) == evec(2) == 0) && ~(evec(1) == evec(3) == 0) && ~(evec(2) == evec(3) == 0) || (sum(evec) == 0)
@@ -40,14 +60,56 @@ if ~(evec(1) == evec(2) == 0) && ~(evec(1) == evec(3) == 0) && ~(evec(2) == evec
 end
 evec0 = evec ./ sum(evec); % evec0 is a unit vector
 
+%% read optional arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+n_conv_arg = 8; % number of conventional arguments
+
+%set defaults
+feed_shift = 0;
+feed_R = 0;
+excitename = '';
+measplanepos = nan;
+
+if (nargin>n_conv_arg)
+    for n=1:2:(nargin-n_conv_arg)
+        if (strcmp(varargin{n},'FeedShift')==1);
+            feed_shift = varargin{n+1};
+            if (numel(feed_shift)>1)
+                error 'FeedShift must be a scalar value'
+            end
+        end
+        
+        if (strcmp(varargin{n},'Feed_R')==1);
+            feed_R = varargin{n+1};
+            if (numel(feed_shift)>1)
+                error 'Feed_R must be a scalar value'
+            end
+        end
+        
+        if (strcmp(varargin{n},'MeasPlaneShift')==1);
+            measplanepos = varargin{n+1};
+            if (numel(feed_shift)>1)
+                error 'MeasPlaneShift must be a scalar value'
+            end
+        end
+                
+        if (strcmp(varargin{n},'ExcitePort')==1);
+            excitename = varargin{n+1};
+        end
+    end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 % normalize start and stop
 nstart = min( [start;stop] );
 nstop  = max( [start;stop] );
 
 % determine index (1, 2 or 3) of propagation (length of MSL)
-idx_prop = dir * [1;2;3];
+idx_prop = dir + 1;
 
 % determine index (1, 2 or 3) of width of MSL
+dir = [0 0 0];
+dir(idx_prop) = 1;
 idx_width = abs(cross(dir,evec0)) * [1;2;3];
 
 % determine index (1, 2 or 3) of height
@@ -66,14 +128,17 @@ MSL_stop = stop;
 MSL_stop(idx_height) = MSL_start(idx_height);
 CSX = AddBox( CSX, materialname, prio, MSL_start, MSL_stop );
 
-% FIXME
+if isnan(measplanepos)
-% openEMS v0.0.7 does not snap PEC
+    measplanepos = (nstart(idx_prop)+nstop(idx_prop))/2;
+else
+    measplanepos = start(idx_prop)+direction*measplanepos;
+end
 
 % calculate position of the voltage probes
 mesh{1} = sort(CSX.RectilinearGrid.XLines);
 mesh{2} = sort(CSX.RectilinearGrid.YLines);
 mesh{3} = sort(CSX.RectilinearGrid.ZLines);
-meshlines = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), (nstart(idx_prop)+nstop(idx_prop))/2, 'nearest' );
+meshlines = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), measplanepos, 'nearest' );
 meshlines = mesh{idx_prop}(meshlines-1:meshlines+1); % get three lines (approx. at center)
 if direction == -1
     meshlines = fliplr(meshlines);
@@ -82,9 +147,9 @@ MSL_w2 = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), (nstart(idx_width)+
 MSL_w2 = mesh{idx_width}(MSL_w2); % get e-line at center of MSL (MSL_width/2)
 v1_start(idx_prop)   = meshlines(1);
 v1_start(idx_width)  = MSL_w2;
-v1_start(idx_height) = nstop(idx_height);
+v1_start(idx_height) = start(idx_height);
 v1_stop  = v1_start;
-v1_stop(idx_height)  = nstart(idx_height);
+v1_stop(idx_height)  = stop(idx_height);
 v2_start = v1_start;
 v2_stop  = v1_stop;
 v2_start(idx_prop)   = meshlines(2);
@@ -111,7 +176,8 @@ i2_stop(idx_prop)    = i2_start(idx_prop);
 
 % create the probes
 name = ['port_ut' num2str(portnr) 'A'];
-weight = sum(evec);
+% weight = sign(stop(idx_height)-start(idx_height))
+weight = -1;
 CSX = AddProbe( CSX, name, 0, weight );
 CSX = AddBox( CSX, name, prio, v1_start, v1_stop );
 name = ['port_ut' num2str(portnr) 'B'];
@@ -121,7 +187,8 @@ name = ['port_ut' num2str(portnr) 'C'];
 CSX = AddProbe( CSX, name, 0, weight );
 CSX = AddBox( CSX, name, prio, v3_start, v3_stop );
 name = ['port_it' num2str(portnr) 'A'];
-weight = direction;
+
+weight = direction
 CSX = AddProbe( CSX, name, 1, weight );
 CSX = AddBox( CSX, name, prio, i1_start, i1_stop );
 name = ['port_it' num2str(portnr) 'B'];
@@ -131,41 +198,18 @@ CSX = AddBox( CSX, name, prio, i2_start, i2_stop );
 % create port structure
 port.nr = portnr;
 port.drawingunit = CSX.RectilinearGrid.ATTRIBUTE.DeltaUnit;
-% port.start = start;
-% port.stop = stop;
-% port.v1_start = v1_start;
-% port.v1_stop = v1_stop;
-% port.v2_start = v2_start;
-% port.v2_stop = v2_stop;
-% port.v3_start = v3_start;
-% port.v3_stop = v3_stop;
 port.v_delta = diff(meshlines);
-% port.i1_start = i1_start;
-% port.i1_stop = i1_stop;
-% port.i2_start = i2_start;
-% port.i2_stop = i2_stop;
 port.i_delta = diff( meshlines(1:end-1) + diff(meshlines)/2 );
 port.direction = direction;
-% port.dir = dir;
-% port.evec = evec;
-% port.idx_prop = idx_prop;
-% port.idx_width = idx_width;
-% port.idx_height = idx_height;
 port.excite = 0;
-port.refplaneshift = 0;
 port.measplanepos = abs(v2_start(idx_prop) - start(idx_prop));
 
-if (nargin >= 9) && (~isempty(refplaneshift))
-    % refplaneshift counts from start of port
-    port.refplaneshift = refplaneshift - direction*(v2_start(idx_prop) - start(idx_prop));
-end
-
 % create excitation
-if nargin >= 10
+if ~isempty(excitename)
 	% excitation of this port is enabled
 	port.excite = 1;
-    meshline = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), start(idx_prop), 'nearest' );
+    meshline = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), start(idx_prop) + feed_shift*direction, 'nearest' );
-    ex_start(idx_prop)   = mesh{idx_prop}(meshline+direction);
+    ex_start(idx_prop)   = mesh{idx_prop}(meshline) ;
 	ex_start(idx_width)  = nstart(idx_width);
 	ex_start(idx_height) = nstart(idx_height);
 	ex_stop(idx_prop)    = ex_start(idx_prop);
@@ -173,4 +217,10 @@ if nargin >= 10
 	ex_stop(idx_height)  = nstop(idx_height);
     CSX = AddExcitation( CSX, excitename, 0, evec );
 	CSX = AddBox( CSX, excitename, prio, ex_start, ex_stop );
+    
+    if feed_R > 0
+        CSX = AddLumpedElement( CSX, [excitename '_R'], idx_height-1, 'R', feed_R );
+        CSX = AddBox( CSX, [excitename '_R'], prio, ex_start, ex_stop );
+        port.Feed_R = port_R;
+    end
 end

matlab/calcPort.m

@@ -1,23 +1,33 @@
-function [S11,beta,ZL,vi] = calcPort( portstruct, SimDir, f, ref_shift )
+function [port] = calcPort( port, SimDir, f, varargin)
-%[S11,beta,ZL,vi] = calcPort( portstruct, SimDir, [f], [ref_shift] )
+% [port] = calcPort( port, SimDir, f, [ref_ZL], [ref_shift])
 %
-% Calculate the reflection coefficient S11, the propagation constant beta
+% Calculate voltages and currents, the propagation constant beta
-% of the MSL-port and the characteristic impedance ZL of the MSL-port.
+% and the characteristic impedance ZL of the given port.
-% The port is to be created by AddMSLPort().
+% The port has to be created by e.g. AddMSLPort().
 %
 % input:
-%   portstruct: return value of AddMSLPort()
+%   port:       return value of AddMSLPort()
 %   SimDir:     directory, where the simulation files are
-%   f:          (optional) frequency vector for DFT
+%   f:          frequency vector for DFT
-%   ref_shift:  (optional) reference plane shift measured from start of port (in drawing units)
+% 
+% 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 
 %
 % output: 
-%   S11:  reflection coefficient (normalized to ZL)
+%   port.f                  the given frequency fector
-%   beta: propagation constant
+%   port.uf.tot/inc/ref     total, incoming and reflected voltage
-%   ZL:   characteristic line impedance
+%   port.if.tot/inc/ref     total, incoming and reflected current
-%   vi:   structure of voltages and currents
+%   port.beta:              propagation constant
-%         vi.TD.v.{val,t}; vi.TD.i.{val,t};
+%   port.ZL:                characteristic line impedance
-%         vi.FD.v.{val,val_shifted,f}; vi.FD.i.{val,val_shifted,f}; 
+%
+% example:
+%   port{1} = calcPort( port{1}, Sim_Path, f, 'RefImpedance', 50);
+%  or 
 % 
 % reference: W. K. Gwarek, "A Differential Method of Reflection Coefficient Extraction From FDTD Simulations",
 %            IEEE Microwave and Guided Wave Letters, Vol. 6, No. 5, May 1996
@@ -28,83 +38,96 @@ function [S11,beta,ZL,vi] = calcPort( portstruct, SimDir, f, ref_shift )
 % See also AddMSLPort
 
 %DEBUG
-% save('/tmp/test.mat', 'portstruct', 'SimDir', 'f', 'nargin' )
+% save('/tmp/test.mat', 'port', 'SimDir', 'f', 'nargin' )
 % load('/tmp/test.mat')
 
 % check
-if portstruct.v_delta(1) ~= portstruct.v_delta(2)
+if abs((port.v_delta(1) - port.v_delta(2)) / port.v_delta(1))>1e-6
 	warning( 'openEMS:calcPort:mesh', 'mesh is not equidistant; expect degraded accuracy' );
 end
 
-if nargin < 3
+
-    f = [];
+%% read optional arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+n_conv_arg = 3; % number of conventional arguments
+
+%set defaults
+ref_ZL = 0;
+ref_shift = nan;
+
+if (nargin>n_conv_arg)
+    for n=1:2:(nargin-n_conv_arg)
+        if (strcmp(varargin{n},'RefPlaneShift')==1);
+            ref_shift = varargin{n+1};
+        end
+        
+        if (strcmp(varargin{n},'RefImpedance')==1);
+            ref_ZL = varargin{n+1};
+        end
+    end
 end
 
 % read time domain data
-filename = ['port_ut' num2str(portstruct.nr)];
+filename = ['port_ut' num2str(port.nr)];
 U = ReadUI( {[filename 'A'],[filename 'B'],[filename 'C']}, SimDir, f );
-filename = ['port_it' num2str(portstruct.nr)];
+filename = ['port_it' num2str(port.nr)];
 I = ReadUI( {[filename 'A'],[filename 'B']}, SimDir, f );
 
-% store the original time domain waveforms
-vi.TD.v = U.TD{2};
-vi.TD.i.t = I.TD{1}.t;
-vi.TD.i.val = (I.TD{1}.val + I.TD{2}.val) / 2; % shift to same position as v
-
 % store the original frequency domain waveforms
-vi.FD.v = U.FD{2};
+u_f = U.FD{2}.val;
-vi.FD.i = I.FD{1};
+i_f = (I.FD{1}.val + I.FD{2}.val) / 2; % shift to same position as v
-vi.FD.i.val = (I.FD{1}.val + I.FD{2}.val) / 2; % shift to same position as v
 
 f = U.FD{2}.f;
 Et = U.FD{2}.val;
-dEt = (U.FD{3}.val - U.FD{1}.val) / (sum(abs(portstruct.v_delta(1:2))) * portstruct.drawingunit);
+dEt = (U.FD{3}.val - U.FD{1}.val) / (sum(abs(port.v_delta(1:2))) * port.drawingunit);
 Ht = (I.FD{1}.val + I.FD{2}.val)/2; % space averaging: Ht is now defined at the same pos as Et
-dHt = (I.FD{2}.val - I.FD{1}.val) / (abs(portstruct.i_delta(1)) * portstruct.drawingunit);
+dHt = (I.FD{2}.val - I.FD{1}.val) / (abs(port.i_delta(1)) * port.drawingunit);
 
 beta = sqrt( - dEt .* dHt ./ (Ht .* Et) );
 beta(real(beta) < 0) = -beta(real(beta) < 0); % determine correct sign (unlike the paper)
 
-% determine S11
-A = sqrt( Et .* dHt ./ (Ht .* dEt) );
-A(imag(A) > 0) = -A(imag(A) > 0); % determine correct sign (unlike the paper)
-S11 = (A - 1) ./ (A + 1);
-
-% determine S11_corrected
-delta_e = sum(portstruct.v_delta(1:2))/2 * portstruct.drawingunit;
-delta_h = portstruct.i_delta(1) * portstruct.drawingunit;
-S11_corrected = sqrt( Et .* (dHt ./ (sin(beta.*delta_h*.5)/(beta*delta_h*.5))) ./ ((Ht ./ cos(beta*delta_h*.5)) .* (dEt ./ (sin(beta*delta_e)./(beta*delta_e)))));
-S11_corrected(imag(S11_corrected) > 0) = -S11_corrected(imag(S11_corrected) > 0); % determine correct sign (unlike the paper)
-S11_corrected = (S11_corrected-1) ./ (S11_corrected+1);
-
-% my own solution...
-temp = sqrt(-dHt .* dEt ./ (Ht .* Et));
-S11 = (-1i * dEt + Et .* temp) ./ (Et .* temp + 1i * dEt); % solution 1
-% S11 = (-1i * dEt - Et .* temp) ./ (-Et .* temp + 1i * dEt); % solution 2
-
-% % determine ZL
-% Et_forward = Et ./ (1 + S11);
-% Ht_forward = Ht ./ (1 - S11);
-% ZL = Et_forward ./ Ht_forward;
-% 
-% % determine ZL_corrected
-% Et_forward_corrected = Et ./ (1 + S11_corrected);
-% Ht_forward_corrected = Ht ./ (1 - S11_corrected);
-% ZL_corrected = Et_forward_corrected ./ Ht_forward_corrected;
-
 % determine ZL
 ZL = sqrt(Et .* dEt ./ (Ht .* dHt));
 
 % reference plane shift (lossless)
-if (nargin > 3)
+if ~isnan(ref_shift)
     % renormalize the shift to the measurement plane
-    ref_shift = ref_shift - portstruct.measplanepos;
+    ref_shift = ref_shift - port.measplanepos * port.drawingunit;
-    ref_shift = ref_shift * portstruct.drawingunit;
+%     ref_shift = ref_shift * port.drawingunit;
-    S11 = S11 .* exp(2i*real(beta)*ref_shift);
-    S11_corrected = S11_corrected .* exp(2i*real(beta)*ref_shift);
     
     % store the shifted frequency domain waveforms
     phase = real(beta)*ref_shift;
-    vi.FD.v.val_shifted = vi.FD.v.val .* cos(-phase) + 1i * vi.FD.i.val.*ZL .* sin(-phase);
+    U.FD{1}.val = u_f .* cos(-phase) + 1i * i_f.*ZL .* sin(-phase);
-    vi.FD.i.val_shifted = vi.FD.i.val .* cos(-phase) + 1i * vi.FD.v.val./ZL .* sin(-phase);
+    I.FD{1}.val = i_f .* cos(-phase) + 1i * u_f./ZL .* sin(-phase);
+    
+    u_f = U.FD{1}.val;
+    i_f = I.FD{1}.val;
 end
+
+if (ref_ZL == 0)
+    if isfield(port,'Feed_R')
+        ref_ZL = port.Feed_R;
+    else
+        ref_ZL = ZL;
+    end
+end
+
+port.ZL =  ZL;
+port.beta = beta;
+
+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;

matlab/examples/__deprecated__/MSL2.m

@@ -53,7 +53,7 @@ min_decrement = 1e-6;
 f_max         = 7e9;
 FDTD = InitFDTD( max_timesteps, min_decrement, 'OverSampling', 10 );
 FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2 );
-BC   = {'MUR' 'PEC' 'PEC' 'PEC' 'PEC' 'PMC'};
+BC   = {'MUR' 'MUR' 'PEC' 'PEC' 'PEC' 'PMC'};
 FDTD = SetBoundaryCond( FDTD, BC );
 
 %% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -161,9 +161,9 @@ grid on
 title( 'Frequency domain current probes' );
 legend( {'abs(if1A)','abs(if1B)','angle(if1A)','angle(if1B)'} );
 
-% port analysis
+%% port analysis
-[S11,beta,ZL,vi] = calcPort( portstruct, Sim_Path, f );
+[U,I,beta,ZL] = calcPort( portstruct, Sim_Path, f );
-% attention! the reflection coefficient S11 is normalized to ZL!
+%% attention! the reflection coefficient S11 is normalized to ZL!
 
 figure
 plot( sin(0:0.01:2*pi), cos(0:0.01:2*pi), 'Color', [.7 .7 .7] );
@@ -216,9 +216,10 @@ grid on;
 legend( {'real','imaginary'}, 'location', 'northeast' )
 title( 'Characteristic line impedance ZL' );
 
-% reference plane shift (to the end of the port)
+%% reference plane shift (to the end of the port)
 ref_shift = abs(portstop(1) - portstart(1));
-[S11,beta,ZL,vi] = calcPort( portstruct, Sim_Path, f, ref_shift );
+[U, I,beta,ZL] = calcPort( portstruct, Sim_Path, f );
+%%
 
 figure
 plotyy( f/1e9, 20*log10(abs(S11)), f/1e9, angle(S11)/pi*180 );