updated matlab port definition functions; Y-parameter calculation
This commit is contained in:
Sebastian Held
parent
3fd58b7e7c
commit
936983c331
@ -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
|
||||
|
||||
@ -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
|
||||
% <refplaneshift> starting from <start> (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
|
||||
|
||||
69
matlab/calc_ypar.m
Normal file
69
matlab/calc_ypar.m
Normal file
@ -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 <sebastian.held@uni-due.de>
|
||||
% 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
|
||||
@ -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)' );
|
||||
|
||||
Loading…
Reference in New Issue
Block a user