From 936983c3310c674f9907fe5a9f37f7084f30b7f8 Mon Sep 17 00:00:00 2001 From: Sebastian Held Date: Fri, 18 Jun 2010 14:26:05 +0200 Subject: [PATCH] updated matlab port definition functions; Y-parameter calculation --- matlab/AddLumpedPort.m | 58 ++++++++++++-------- matlab/AddMSLPort.m | 73 +++++++++++++++---------- matlab/calc_ypar.m | 69 +++++++++++++++++++++++ matlab/examples/gauss_excitation_test.m | 29 ++++++++-- 4 files changed, 171 insertions(+), 58 deletions(-) create mode 100644 matlab/calc_ypar.m diff --git a/matlab/AddLumpedPort.m b/matlab/AddLumpedPort.m index 83cd0f8..ceb5e2c 100644 --- a/matlab/AddLumpedPort.m +++ b/matlab/AddLumpedPort.m @@ -6,7 +6,7 @@ function [CSX,port] = AddLumpedPort( CSX, portnr, R, start, stop, dir, excitenam % R: internal resistance of the port % 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 of port (choices: [1 0 0], [0 1 0] or [0 0 1]) % excitename (optional): if specified, the port will be switched on (see AddExcitation()) % % the mesh must be already initialized @@ -74,21 +74,25 @@ delta2_n = mesh{idx_plane}(idx) - mesh{idx_plane}(idx-1); delta2_p = mesh{idx_plane}(idx+1) - mesh{idx_plane}(idx); m_start = nstart; m_stop = nstop; -m_start = m_start - delta2_n/2; -m_stop = m_stop + delta2_p/2; +m_start(idx_plane) = m_start(idx_plane) - delta2_n/2; +m_stop(idx_plane) = m_stop(idx_plane) + delta2_p/2; % calculate kappa l = (m_stop(idx_cal) - m_start(idx_cal)) * drawingunit; % length of the "sheet" A = (m_stop(idx1) - m_start(idx1)) * (m_stop(idx_plane) - m_start(idx_plane)) * drawingunit^2; % area of the "sheet" kappa = l/A / R; % [kappa] = S/m -CSX = AddMaterial( CSX, ['port' num2str(portnr) '_sheet_resistance'] ); -CSX = SetMaterialProperty( CSX, ['port' num2str(portnr) '_sheet_resistance'], 'Kappa', kappa ); +CSX = AddMaterial( CSX, ['port' num2str(portnr) '_sheet_resistance'], 'Isotropy', 0 ); +kappa_cell = {}; +kappa_cell{1} = kappa*dir(1); +kappa_cell{2} = kappa*dir(2); +kappa_cell{3} = kappa*dir(3); +CSX = SetMaterialProperty( CSX, ['port' num2str(portnr) '_sheet_resistance'], 'Kappa', kappa_cell ); CSX = AddBox( CSX, ['port' num2str(portnr) '_sheet_resistance'], 0, m_start, m_stop ); % calculate position of the voltage probe -center1 = interp1( mesh{idx1}, 1:numel(mesh{idx1}), (nstart(idx1)+nstop(idx1))/2, 'nearest' ); v_start(idx_plane) = start(idx_plane); -v_start(idx1) = center1; +center1 = interp1( mesh{idx1}, 1:numel(mesh{idx1}), (nstart(idx1)+nstop(idx1))/2, 'nearest' ); +v_start(idx1) = mesh{idx1}(center1); v_stop = v_start; v_start(idx_cal) = nstart(idx_cal); v_stop(idx_cal) = nstop(idx_cal); @@ -98,9 +102,10 @@ idx = interp1( mesh{idx1}, 1:numel(mesh{idx1}), nstart(idx1), 'nearest' ); delta1_n = mesh{idx1}(idx) - mesh{idx1}(idx-1); idx = interp1( mesh{idx1}, 1:numel(mesh{idx1}), nstop(idx1), 'nearest' ); delta1_p = mesh{idx1}(idx+1) - mesh{idx1}(idx); -idx = interp1( mesh{idx_cal}, 1:numel(mesh{idx_cal}), (nstart(idx_cal)+nstop(idx_cal))/2, 'nearest' ); -i_start(idx_cal) = mesh{idx_cal}(idx-1); -i_stop(idx_cal) = mesh{idx_cal}(idx-1); +h_offset = diff(mesh{idx_cal}); +idx = interp1( mesh{idx_cal} + [h_offset h_offset(end)]/2, 1:numel(mesh{idx_cal}), (nstart(idx_cal)+nstop(idx_cal))/2, 'nearest' ); +i_start(idx_cal) = mesh{idx_cal}(idx) + h_offset(idx)/2; +i_stop(idx_cal) = i_start(idx_cal); i_start(idx1) = nstart(idx1) - delta1_n/2; i_start(idx_plane) = nstart(idx_plane) - delta2_n/2; i_stop(idx1) = nstop(idx1) + delta1_p/2; @@ -108,31 +113,38 @@ i_stop(idx_plane) = nstop(idx_plane) + delta2_p/2; % create the probes name = ['port_ut' num2str(portnr)]; -CSX = AddProbe( CSX, name, 0 ); +weight = -direction; +CSX = AddProbe( CSX, name, 0, weight ); CSX = AddBox( CSX, name, 999, v_start, v_stop ); name = ['port_it' num2str(portnr)]; -CSX = AddProbe( CSX, name, 1 ); +weight = direction; +CSX = AddProbe( CSX, name, 1, weight ); CSX = AddBox( CSX, name, 999, i_start, i_stop ); % create port structure port.nr = portnr; port.drawingunit = CSX.RectilinearGrid.ATTRIBUTE.DeltaUnit; -port.start = start; -port.stop = stop; -port.v_start = v_start; -port.v_stop = v_stop; -port.i_start = i_start; -port.i_stop = i_stop; -port.dir = dir; +% port.start = start; +% port.stop = stop; +% port.v_start = v_start; +% port.v_stop = v_stop; +% port.i_start = i_start; +% port.i_stop = i_stop; +% port.dir = dir; port.direction = direction; -port.idx_cal = idx_cal; -port.idx1 = idx1; -port.idx1 = idx1; +% port.idx_cal = idx_cal; +% port.idx1 = idx1; +% port.idx1 = idx1; port.excite = 0; % create excitation if (nargin >= 7) && ~isempty(excitename) % excitation of this port is enabled port.excite = 1; - CSX = AddBox( CSX, excitename, 999, v_start, v_stop ); + e_start = nstart; + e_stop = nstop; + e_start(idx_plane) = start(idx_plane); % excitation-plane is determined by start vector + e_stop(idx_plane) = start(idx_plane); + CSX = AddExcitation( CSX, excitename, 0, -dir*direction); + CSX = AddBox( CSX, excitename, 999, e_start, e_stop ); end diff --git a/matlab/AddMSLPort.m b/matlab/AddMSLPort.m index ca09ec4..591cf70 100644 --- a/matlab/AddMSLPort.m +++ b/matlab/AddMSLPort.m @@ -1,4 +1,4 @@ -function [CSX,port] = AddMSLPort( CSX, portnr, materialname, start, stop, dir, evec, excitename ) +function [CSX,port] = AddMSLPort( CSX, portnr, materialname, start, stop, dir, evec, refplaneshift, excitename ) % [CSX,port] = AddMSLPort( CSX, portnr, materialname, start, stop, dir, evec, excitename ) % % CSX: CSX-object created by InitCSX() @@ -8,13 +8,17 @@ function [CSX,port] = AddMSLPort( CSX, portnr, materialname, start, stop, dir, e % stop: 3D end rowvector for port definition % dir: direction of wave propagation (choices: [1 0 0], [0 1 0] or [0 0 1]) % 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()) % % 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) +% 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) % @@ -33,7 +37,7 @@ dir = dir ./ sum(dir); % dir is now a unit vector if ~(evec(1) == evec(2) == 0) && ~(evec(1) == evec(3) == 0) && ~(evec(2) == evec(3) == 0) || (sum(evec) == 0) error 'evec must have exactly one component ~= 0' end -evec0 = evec ./ abs(sum(evec)); % evec0 is a unit vector +evec0 = evec ./ sum(evec); % evec0 is a unit vector % normalize start and stop nstart = min( [start;stop] ); @@ -106,56 +110,65 @@ i2_stop(idx_prop) = i2_start(idx_prop); % create the probes name = ['port_ut' num2str(portnr) 'A']; -CSX = AddProbe( CSX, name, 0 ); +weight = sum(evec); +CSX = AddProbe( CSX, name, 0, weight ); CSX = AddBox( CSX, name, 999, v1_start, v1_stop ); name = ['port_ut' num2str(portnr) 'B']; -CSX = AddProbe( CSX, name, 0 ); +CSX = AddProbe( CSX, name, 0, weight ); CSX = AddBox( CSX, name, 999, v2_start, v2_stop ); name = ['port_ut' num2str(portnr) 'C']; -CSX = AddProbe( CSX, name, 0 ); +CSX = AddProbe( CSX, name, 0, weight ); CSX = AddBox( CSX, name, 999, v3_start, v3_stop ); name = ['port_it' num2str(portnr) 'A']; -CSX = AddProbe( CSX, name, 1 ); +weight = direction; +CSX = AddProbe( CSX, name, 1, weight ); CSX = AddBox( CSX, name, 999, i1_start, i1_stop ); name = ['port_it' num2str(portnr) 'B']; -CSX = AddProbe( CSX, name, 1 ); +CSX = AddProbe( CSX, name, 1, weight ); CSX = AddBox( CSX, name, 999, 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.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.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.dir = dir; -port.evec = evec; -port.idx_prop = idx_prop; -port.idx_width = idx_width; -port.idx_height = idx_height; +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; + +if (nargin >= 8) && (~isempty(refplaneshift)) + % refplaneshift counts from start of port + port.refplaneshift = refplaneshift - direction*(v2_start(idx_prop) - start(idx_prop)); +end % create excitation -if nargin >= 8 +if nargin >= 9 % 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*2); % excitation is placed two cells away from the start of the port (to be able to use the MUR_ABC) - ex_start(idx_prop) = start(idx_prop); + meshline = interp1( mesh{idx_prop}, 1:numel(mesh{idx_prop}), start(idx_prop), 'nearest' ); + ex_start(idx_prop) = mesh{idx_prop}(meshline+direction); ex_start(idx_width) = nstart(idx_width); ex_start(idx_height) = nstart(idx_height); ex_stop(idx_prop) = ex_start(idx_prop); ex_stop(idx_width) = nstop(idx_width); ex_stop(idx_height) = nstop(idx_height); + CSX = AddExcitation( CSX, excitename, 0, evec ); CSX = AddBox( CSX, excitename, 999, ex_start, ex_stop ); end diff --git a/matlab/calc_ypar.m b/matlab/calc_ypar.m new file mode 100644 index 0000000..f6ae1bd --- /dev/null +++ b/matlab/calc_ypar.m @@ -0,0 +1,69 @@ +function Y = calc_ypar( f, ports, Sim_Path_Prefix ) +% Y = calc_ypar( f, ports, Sim_Path_Prefix ) +% +% f: frequency vector (Hz) +% ports: cell array of ports (see AddMSLPort() and AddLumpedPort()) +% Sim_Path_Prefix: prefix of the simulation dirs (will be postfixed by +% excitation port number) +% +% This function calculates the Y-matrix representation of the ports +% +% It is assumed that each port (inside ports) is excited and the +% corresponding simulation was carried out at Sim_Path + portnr (e.g. for +% port 2: '/tmp/sim2') +% +% Sebastian Held +% Jun 9 2010 +% +% See also AddMSLPort AddLumpedPort + +% sanitize input arguments +f = reshape(f,1,[]); % make it a row vector + +% prepare result matrix +maxportnr = max( cellfun(@(x) x.nr, ports) ); +Y = ones(maxportnr,maxportnr,numel(f)) * NaN; +U = ones(maxportnr,maxportnr,numel(f)) * NaN; +I = ones(maxportnr,maxportnr,numel(f)) * NaN; + +% read time domain simulation results +for runnr = 1:numel(ports) + Sim_Path = [Sim_Path_Prefix num2str(ports{runnr}.nr)]; + for pnr = 1:numel(ports) + if isfield( ports{pnr}, 'v_delta' ) + % this is an MSLPort + temp_U = ReadUI( ['port_ut' num2str(ports{pnr}.nr) 'B'], Sim_Path ); + temp = ReadUI( {['port_it' num2str(ports{pnr}.nr) 'A'],['port_it' num2str(ports{pnr}.nr) 'B']}, Sim_Path ); + temp_I.TD{1}.t = temp.TD{1}.t; + temp_I.TD{1}.val = (temp.TD{1}.val + temp.TD{2}.val) / 2; % space averaging + else + % this is a lumped port + temp_U = ReadUI( ['port_ut' num2str(ports{pnr}.nr)], Sim_Path ); + temp_I = ReadUI( ['port_it' num2str(ports{pnr}.nr)], Sim_Path ); + +% % correct the orientation of the probes (FIXME to be done inside +% % openEMS) +% temp_U.TD{1}.val = temp_U.TD{1}.val * (-ports{pnr}.direction); + end + +% % correct the orientation of the probes (FIXME to be done inside +% % openEMS) +% temp_I.TD{1}.val = temp_I.TD{1}.val * ports{pnr}.direction; +% if runnr == 5 % DEBUG +% temp_I.TD{1}.val = temp_I.TD{1}.val * -1; +% end + + % time domain -> frequency domain + U(ports{pnr}.nr,ports{runnr}.nr,:) = DFT_time2freq( temp_U.TD{1}.t, temp_U.TD{1}.val, f ); + I(ports{pnr}.nr,ports{runnr}.nr,:) = DFT_time2freq( temp_I.TD{1}.t, temp_I.TD{1}.val, f ); + + % compensate H-field time advance + delta_t_2 = temp_I.TD{1}.t(1) - temp_U.TD{1}.t(1); % half time-step (s) + I(ports{pnr}.nr,ports{runnr}.nr,:) = squeeze(I(ports{pnr}.nr,ports{runnr}.nr,:)).' .* exp(-1i*2*pi*f*delta_t_2); + end +end + +% calc Y-parameters +for a=1:numel(f) + Y(:,:,a) = I(:,:,a) / U(:,:,a); +end diff --git a/matlab/examples/gauss_excitation_test.m b/matlab/examples/gauss_excitation_test.m index 87bc25f..aea2b89 100644 --- a/matlab/examples/gauss_excitation_test.m +++ b/matlab/examples/gauss_excitation_test.m @@ -1,10 +1,19 @@ +% +% this script evaluates the same gaussian excitation function, as openEMS does +% + clear close all -f0 = 1e9/2; -fc = 1e9/2; +f0 = 0; +fc = 10e9; dT = 8e-12; % sample time-step + + + + + len = 2 * 9/(2*pi*fc) / dT; % gauss length for n=1:len @@ -22,16 +31,26 @@ val = DFT_time2freq( t_, ex, [f0-fc f0 f0+fc] ); disp( ['Amplitude at f=f0-fc: ' num2str(20*log10(abs(val(1))/abs(val(2)))) ' dB'] ); disp( ['Amplitude at f=f0+fc: ' num2str(20*log10(abs(val(3))/abs(val(2)))) ' dB'] ); +% calculate frequency domain via slow DFT freq = linspace(f0-fc,f0+fc,1000); val = DFT_time2freq( t_, ex, freq ); figure plot( freq/1e9, abs(val) ) -[f,val] = FFT_time2freq( t_, ex ); -val = val((f0-fc<=f) & (f<=f0+fc)); +% overlay the FFT result +[f,val_fft] = FFT_time2freq( t_, ex ); +val_fft = val_fft((f0-fc<=f) & (f<=f0+fc)); f = f((f0-fc<=f) & (f<=f0+fc)); hold on -plot( f/1e9, abs(val), 'r' ) +plot( f/1e9, abs(val_fft), 'r' ) xlabel( 'frequency (GHz)' ); ylabel( 'amplitude' ); + +% dB +figure +val = val(freq>=0); +freq = freq(freq>=0); +plot( freq/1e9, 20*log10(abs(val)/max(abs(val))), 'r' ) +xlabel( 'frequency (GHz)' ); +ylabel( 'amplitude (dB)' );