From dd7269d40a4ba3d6fa3160cf474168133c0c8c27 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Mon, 19 Sep 2011 10:14:27 +0200 Subject: [PATCH] matlab: revision of AddMSLPort + calcPort all examples using this functions need to be revised! --- matlab/AddMSLPort.m | 156 ++++++++++++------ matlab/calcPort.m | 153 +++++++++-------- .../{microstrip => __deprecated__}/MSL2.m | 13 +- 3 files changed, 198 insertions(+), 124 deletions(-) rename matlab/examples/{microstrip => __deprecated__}/MSL2.m (96%) diff --git a/matlab/AddMSLPort.m b/matlab/AddMSLPort.m index c7dfa75..e0d0bc0 100644 --- a/matlab/AddMSLPort.m +++ b/matlab/AddMSLPort.m @@ -1,38 +1,58 @@ -function [CSX,port] = AddMSLPort( CSX, prio, portnr, materialname, start, stop, dir, evec, refplaneshift, excitename ) -% [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, 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 -% starting from (thus refplaneshift is normally -% positive) -% excitename (optional): if specified, the port will be switched on (see AddExcitation()) -% +% +% variable input: +% 'ExcitePort' necessary excitation name to make the port an active feeding port +% 'FeedShift' shift to port from start by a given distance in drawing +% 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] -% this defines a MSL in x-direction (dir) with an e-field excitation in -z-direction (evec) -% 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) +% start = [0 0 height]; +% stop = [length width 0]; +% CSX = AddMetal( CSX, 'metal' ); %create a PEC called 'metal' +% [CSX,port] = AddMSLPort( CSX, 0, 1, 'metal', start, stop, ... +% 0, [0 0 -1] , 'ExcitePort', 'excite', 'Feed_R', 50 ) +% +% this defines a MSL in x-direction (dir=0) with an e-field excitation +% 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) % -% Sebastian Held -% May 13 2010 +% Sebastian Held May 13 2010 +% Thorsten Liebig Sept 16 2011 % -% See also InitCSX AddMetal AddMaterial AddExcitation calcMSLPort +% 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 -% openEMS v0.0.7 does not snap PEC +if isnan(measplanepos) + 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' ); - ex_start(idx_prop) = mesh{idx_prop}(meshline+direction); + 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); 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 diff --git a/matlab/calcPort.m b/matlab/calcPort.m index 782e2e6..f6a5b9b 100644 --- a/matlab/calcPort.m +++ b/matlab/calcPort.m @@ -1,24 +1,34 @@ -function [S11,beta,ZL,vi] = calcPort( portstruct, SimDir, f, ref_shift ) -%[S11,beta,ZL,vi] = calcPort( portstruct, SimDir, [f], [ref_shift] ) +function [port] = calcPort( port, SimDir, f, varargin) +% [port] = calcPort( port, SimDir, f, [ref_ZL], [ref_shift]) % -% Calculate the reflection coefficient S11, the propagation constant beta -% of the MSL-port and the characteristic impedance ZL of the MSL-port. -% The port is to be created by AddMSLPort(). +% Calculate voltages and currents, the propagation constant beta +% and the characteristic impedance ZL of the given port. +% 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 -% ref_shift: (optional) reference plane shift measured from start of port (in drawing units) +% 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 % -% output: -% S11: reflection coefficient (normalized to ZL) -% beta: propagation constant -% ZL: characteristic line impedance -% vi: structure of voltages and currents -% vi.TD.v.{val,t}; vi.TD.i.{val,t}; -% vi.FD.v.{val,val_shifted,f}; vi.FD.i.{val,val_shifted,f}; +% 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} = 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}; -vi.FD.i = I.FD{1}; -vi.FD.i.val = (I.FD{1}.val + I.FD{2}.val) / 2; % shift to same position as v +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(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 * portstruct.drawingunit; - S11 = S11 .* exp(2i*real(beta)*ref_shift); - S11_corrected = S11_corrected .* exp(2i*real(beta)*ref_shift); + ref_shift = ref_shift - port.measplanepos * port.drawingunit; +% ref_shift = ref_shift * port.drawingunit; % 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); - vi.FD.i.val_shifted = vi.FD.i.val .* cos(-phase) + 1i * vi.FD.v.val./ZL .* sin(-phase); + 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; diff --git a/matlab/examples/microstrip/MSL2.m b/matlab/examples/__deprecated__/MSL2.m similarity index 96% rename from matlab/examples/microstrip/MSL2.m rename to matlab/examples/__deprecated__/MSL2.m index d12f4f1..31a2600 100644 --- a/matlab/examples/microstrip/MSL2.m +++ b/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 -[S11,beta,ZL,vi] = calcPort( portstruct, Sim_Path, f ); -% attention! the reflection coefficient S11 is normalized to ZL! +%% port analysis +[U,I,beta,ZL] = calcPort( portstruct, Sim_Path, f ); +%% 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 );