matlab: update handling ports

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

matlab/AddCurvePort.m

@@ -27,6 +27,9 @@ function [CSX,port] = AddCurvePort( CSX, prio, portnr, R, start, stop, excite, v
 % (C) 2010 Sebastian Held <sebastian.held@uni-due.de>
 % See also InitCSX AddExcitation
 
+port.type='Lumped';
+port.nr=portnr;
+
 % make row vector
 start = reshape( start, 1, [] );
 stop  = reshape( stop , 1, [] );
@@ -98,12 +101,7 @@ m_stop(dir1)  = m_stop(dir1)  + delta1_p/2;
 m_start(dir2) = m_start(dir2) - delta2_n/2;
 m_stop(dir2)  = m_stop(dir2)  + delta2_p/2;
 
-% calculate kappa
+% calculate position of the voltage probe & excitation
-materialname = ['port' num2str(portnr) '_sheet_resistance'];
-CSX = AddLumpedElement( CSX, materialname, dir-1, 'R', R);
-CSX = AddBox( CSX, materialname, prio, m_start, m_stop );
-
-% calculate position of the voltage probe
 v_start = [mesh{1}(port_start_idx(1)), mesh{2}(port_start_idx(2)), mesh{3}(port_start_idx(3))];
 v_stop  = [mesh{1}(port_stop_idx(1)),  mesh{2}(port_stop_idx(2)),  mesh{3}(port_stop_idx(3))];
 
@@ -124,7 +122,6 @@ CSX = AddProbe( CSX, name, 1, weight );
 CSX = AddBox( CSX, name, prio, i_start, i_stop );
 
 % create port structure
-port.nr = portnr;
 port.drawingunit = unit;
 % port.start = start;
 % port.stop = stop;
@@ -137,7 +134,6 @@ port.drawingunit = unit;
 % port.idx_cal = idx_cal;
 % port.idx1 = idx1;
 % port.idx1 = idx1;
-port.excite = 0;
 
 if (nargin < 7)
     excite = false;
@@ -153,6 +149,8 @@ if ischar(excite)
     end
 end
 
+port.excite = excite;
+
 % create excitation
 if (excite)
 	% excitation of this port is enabled
@@ -162,3 +160,11 @@ if (excite)
     CSX = AddExcitation( CSX, ['port_excite_' num2str(portnr)], 0, start_idx ~= stop_idx, varargin{:});
 	CSX = AddBox( CSX, ['port_excite_' num2str(portnr)], prio, e_start, e_stop );
 end
+
+port.Feed_R = R;
+if (R>0 && (~isinf(R))) 
+    CSX = AddLumpedElement( CSX, ['port_resist_' int2str(portnr)], dir-1, 'R', R);
+    CSX = AddBox( CSX, ['port_resist_' int2str(portnr)], prio, v_start, v_stop );
+elseif (R==0)
+    CSX = AddBox(CSX,metalname, prio, v_start, v_stop);
+end

matlab/AddLumpedPort.m

@@ -1,5 +1,5 @@
-function [CSX] = AddLumpedPort( CSX, prio, portnr, R, start, stop, dir, excite, varargin )
+function [CSX, port] = AddLumpedPort( CSX, prio, portnr, R, start, stop, dir, excite, varargin )
-% [CSX] = AddLumpedPort( CSX, prio, portnr, R, start, stop, dir, excite, varargin )
+% [CSX, port] = AddLumpedPort( CSX, prio, portnr, R, start, stop, dir, excite, varargin )
 %
 % Add a 3D lumped port as an excitation.
 %
@@ -31,6 +31,9 @@ function [CSX] = AddLumpedPort( CSX, prio, portnr, R, start, stop, dir, excite,
 
 % check dir
 
+port.type='Lumped';
+port.nr=portnr;
+
 if (dir(1)~=0) && (dir(2) == 0) && (dir(3)==0)
     n_dir = 1;
 elseif (dir(1)==0) && (dir(2) ~= 0) && (dir(3)==0)
@@ -50,7 +53,9 @@ if (stop(n_dir)-start(n_dir)) > 0
 else
     direction = -1;
 end
+port.direction = direction;
 
+port.Feed_R = R;
 if (R>0 && (~isinf(R)))
     CSX = AddLumpedElement(CSX,['port_resist_' int2str(portnr)],  n_dir-1, 'Caps', 1, 'R', R);
     CSX = AddBox(CSX,['port_resist_' int2str(portnr)], prio, start, stop);
@@ -73,7 +78,7 @@ if ischar(excite)
     end
 end
 
-
+port.excite = excite;
 % create excitation
 if (excite)
     CSX = AddExcitation( CSX, ['port_excite_' num2str(portnr)], 0, -dir*direction, varargin{:});

matlab/AddMSLPort.m

@@ -66,7 +66,7 @@ n_conv_arg = 8; % number of conventional arguments
 
 %set defaults
 feed_shift = 0;
-feed_R = 0;
+feed_R = inf; %(default is open, no resitance)
 excite = false;
 measplanepos = nan;
 
@@ -215,6 +215,7 @@ if ((CSX.ATTRIBUTE.CoordSystem==1) && (idx_prop==2))
     port.LengthScale = MSL_stop(idx_height);
 end
 port.nr = portnr;
+port.type = 'MSL';
 port.drawingunit = CSX.RectilinearGrid.ATTRIBUTE.DeltaUnit;
 port.v_delta = diff(meshlines)*port.LengthScale;
 port.i_delta = diff( meshlines(1:end-1) + diff(meshlines)/2 )*port.LengthScale;
@@ -223,9 +224,7 @@ port.excite = 0;
 port.measplanepos = abs(v2_start(idx_prop) - start(idx_prop))*port.LengthScale;
 % port
 
-% create excitation
+% create excitation (if enabled) and port resistance
-% excitation of this port is enabled
-port.excite = 1;
 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) ;
 ex_start(idx_width)  = nstart(idx_width);
@@ -234,13 +233,26 @@ ex_stop(idx_prop)    = ex_start(idx_prop);
 ex_stop(idx_width)   = nstop(idx_width);
 ex_stop(idx_height)  = nstop(idx_height);
 
+port.excite = 0;
 if excite
+    port.excite = 1;
     CSX = AddExcitation( CSX, ['port_excite_' num2str(portnr)], 0, evec, excite_args{:} );
     CSX = AddBox( CSX, ['port_excite_' num2str(portnr)], prio, ex_start, ex_stop );
 end
-if feed_R > 0
+
+%% MSL resitance at start of MSL line
+ex_start(idx_prop) = start(idx_prop);
+ex_stop(idx_prop) = ex_start(idx_prop);
+
+if (feed_R > 0) && ~isinf(feed_R)
     CSX = AddLumpedElement( CSX, 'port_R', idx_height-1, 'R', feed_R );
     CSX = AddBox( CSX, 'port_R', prio, ex_start, ex_stop );
-    port.Feed_R = feed_R;
+elseif isinf(feed_R)
+    % do nothing --> open port
+elseif feed_R == 0
+    %port "resistance" as metal
+    CSX = AddBox( CSX, materialname, prio, ex_start, ex_stop );
+else
+    error('openEMS:AddMSLPort','MSL port with resitance <= 0 it not possible');
 end
 end

matlab/calcLumpedPort.m

@@ -0,0 +1,84 @@
+function [port] = calcLumpedPort( port, SimDir, f, varargin)
+% [port] = calcLumpedPort( port, SimDir, f, varargin)
+%
+% Calculate voltages and currents of given lumped port.
+%
+% The port has to be created by e.g. AddLumpedPort().
+%
+% input:
+%   port:       return value of e.g. AddMSLPort()
+%   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
+%
+% 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
+%
+% example:
+%   port{1} = calcLumpedPort( port{1}, Sim_Path, f, 'RefImpedance', 50);
+%
+% openEMS matlab interface
+% -----------------------
+% (C) 2012 Thorsten Liebig <thorsten.liebig@gmx.de>
+%
+% See also AddLumpedPort, calcPort
+
+if (iscell(port))
+    for n=1:numel(port)
+        port{n}=calcLumpedPort(port{n}, SimDir, f, varargin{:});
+    end
+    return;
+end
+
+if (strcmpi(port.type,'Lumped')~=1 && strcmpi(port.type,'Curve')~=1)
+    error('openEMS:calcLumpedPort','error, type is not a lumped port');
+end
+
+
+%% read optional arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+n_conv_arg = 3; % number of conventional arguments
+
+%set defaults
+ref_ZL = port.Feed_R;
+
+if (nargin>n_conv_arg)
+    for n=1:2:(nargin-n_conv_arg)
+        if (strcmp(varargin{n},'RefImpedance')==1);
+            ref_ZL = varargin{n+1};
+        end
+    end
+end
+
+% read time domain data
+filename = ['port_ut' num2str(port.nr)];
+U = ReadUI(filename, SimDir, f );
+filename = ['port_it' num2str(port.nr)];
+I = ReadUI(filename, SimDir, f );
+
+% store the original frequency domain waveforms
+u_f = U.FD{1}.val;
+i_f = I.FD{1}.val; % shift to same position as v
+
+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/calcPort.m

@@ -1,9 +1,11 @@
 function [port] = calcPort( port, SimDir, f, varargin)
 % [port] = calcPort( port, SimDir, f, varargin)
 %
-% Calculate voltages and currents, the propagation constant beta
+% Calculate:
-% and the characteristic impedance ZL of the given port.
+%   - voltages and currents
-% The port has to be created by e.g. AddMSLPort().
+%   - the propagation constant and the characteristic impedance (if applicable)
+%
+% The port has to be created by e.g. AddMSLPort(), AddLumpedPort() or AddCurvePort
 %
 % input:
 %   port:       return value of AddMSLPort()
@@ -14,122 +16,39 @@ function [port] = calcPort( port, SimDir, f, varargin)
 %   '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
+%   'RefPlaneShift': for transmission lines only, See also calcTLPort for
-%                      for a desired phase correction
+%                    more details
-%                    - default is the measurement plane
-%                    - 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
+%
+%   if port is a transmission line port:
 %   port.beta:              propagation constant
 %   port.ZL:                characteristic line impedance
 %
 % example:
-%   port{1} = calcPort( port{1}, Sim_Path, f, 'RefImpedance', 50);
+%   port = calcPort(port, 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
 %
 % openEMS matlab interface
 % -----------------------
-% (C) 2010 Sebastian Held <sebastian.held@uni-due.de>
+% (C) 2012 Thorsten Liebig <thorsten.liebig@gmx.de>
-% See also AddMSLPort
+%
+% See also AddMSLPort, AddLumpedPort, AddCurvePort, calcTLPort, calcLumpedPort
 
-%DEBUG
+if (iscell(port))
-% save('/tmp/test.mat', 'port', 'SimDir', 'f', 'nargin' )
+    for n=1:numel(port)
-% load('/tmp/test.mat')
+        port{n}=calcPort(port{n}, SimDir, f, varargin{:});
-
+    end
-% check
+    return;
-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 strcmpi(port.type,'MSL')
-%% read optional arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    port = calcTLPort( port, SimDir, f, varargin{:});
-n_conv_arg = 3; % number of conventional arguments
+    return
-
+elseif (strcmpi(port.type,'Lumped') || strcmpi(port.type,'Curve'))
-%set defaults
+    port = calcLumpedPort( port, SimDir, f, varargin{:});
-ref_ZL = 0;
+    return
-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(port.nr)];
-U = ReadUI( {[filename 'A'],[filename 'B'],[filename 'C']}, SimDir, f );
-filename = ['port_it' num2str(port.nr)];
-I = ReadUI( {[filename 'A'],[filename 'B']}, SimDir, f );
-
-% store the original frequency domain waveforms
-u_f = U.FD{2}.val;
-i_f = (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(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(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 ZL
-ZL = sqrt(Et .* dEt ./ (Ht .* dHt));
-
-% reference plane shift (lossless)
-if ~isnan(ref_shift)
-    ref_shift = ref_shift * port.LengthScale;
-    % shift to the beginning of MSL
-    ref_shift = ref_shift - port.measplanepos;
-    ref_shift = ref_shift * port.drawingunit;
-    
-    % store the shifted frequency domain waveforms
-    phase = real(beta)*ref_shift;
-    U.FD{1}.val = u_f .* cos(-phase) + 1i * i_f.*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/calcTLPort.m

@@ -0,0 +1,137 @@
+function [port] = calcTLPort( 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 transmission line port.
+%
+% The port has to be created by e.g. AddMSLPort().
+%
+% input:
+%   port:       return value of e.g. AddMSLPort()
+%   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
+%                    - 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
+%
+% example:
+%   port{1} = calcTLPort( port{1}, Sim_Path, f, 'RefImpedance', 50);
+%
+% 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
+%
+% openEMS matlab interface
+% -----------------------
+% (C) 2010 Sebastian Held <sebastian.held@uni-due.de>
+%
+% See also AddMSLPort, calcPort
+
+if (iscell(port))
+    for n=1:numel(port)
+        port{n}=calcTLPort(port{n}, SimDir, f, varargin{:});
+    end
+    return;
+end
+
+if (strcmpi(port.type,'MSL')~=1)
+    error('openEMS:calcTLPort','error, type is not a transmission line port');
+end
+
+% check
+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
+
+
+%% read optional arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+n_conv_arg = 3; % number of conventional arguments
+
+%set defaults
+ref_ZL = -1;
+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(port.nr)];
+U = ReadUI( {[filename 'A'],[filename 'B'],[filename 'C']}, SimDir, f );
+filename = ['port_it' num2str(port.nr)];
+I = ReadUI( {[filename 'A'],[filename 'B']}, SimDir, f );
+
+% store the original frequency domain waveforms
+u_f = U.FD{2}.val;
+i_f = (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(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(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 ZL
+ZL = sqrt(Et .* dEt ./ (Ht .* dHt));
+
+% reference plane shift (lossless)
+if ~isnan(ref_shift)
+    ref_shift = ref_shift * port.LengthScale;
+    % shift to the beginning of MSL
+    ref_shift = ref_shift - port.measplanepos;
+    ref_shift = ref_shift * port.drawingunit;
+
+    % store the shifted frequency domain waveforms
+    phase = real(beta)*ref_shift;
+    U.FD{1}.val = u_f .* cos(-phase) + 1i * i_f.*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)
+    ref_ZL = ZL;
+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;