matlab: new coaxial port and example

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>
This commit is contained in:
Thorsten Liebig 2013-07-22 23:00:26 +02:00
parent 68e20f06ea
commit 2b3408b0b8
4 changed files with 262 additions and 95 deletions

232
matlab/AddCoaxialPort.m Normal file
View File

@ -0,0 +1,232 @@
function [CSX,port] = AddCoaxialPort( CSX, prio, portnr, pec_name, materialname, start, stop, dir, r_i, r_o, r_os, varargin )
% function [CSX,port] = AddCoaxialPort( CSX, prio, portnr, pec_name, materialname, start, stop, dir, r_i, r_o, r_os, varargin )
%
% CSX: CSX-object created by InitCSX()
% prio: priority for excitation and probe boxes
% portnr: (integer) number of the port
% pec_name: metal property for coaxial inner/outer conductor (created by AddMetal())
% materialname: substrate property for coaxial line (created by AddMaterial())
% Note: this may be empty for an "air filled" coaxial line
% start: 3D start rowvector for coaxial cable axis
% stop: 3D end rowvector for coaxial cable axis
% dir: direction of wave propagation (choices: 0, 1, 2 or 'x','y','z')
% r_i: inner coaxial radius (in drawing unit)
% r_o: outer coaxial radius (in drawing unit)
% r_os: outer shell coaxial radius (in drawing unit)
%
% variable input:
% varargin: optional additional excitations options, see also AddExcitation
% 'ExciteAmp' excitation amplitude of transversal electric field profile,
% set to 0 (default) for a passive port
% 'FeedShift' shift to port from start by a given distance in drawing
% units. Default is 0. Only active if 'ExciteAmp' is set!
% 'Feed_R' Specifiy a lumped port resistance. Default is no lumped
% port resistance --> port has to end in an ABC.
% 'MeasPlaneShift' Shift the measurement plane from start t a given distance
% in drawing units. Default is the middle of start/stop.
% 'PortNamePrefix' a prefix to the port name
%
% the mesh must be already initialized
%
% example:
%
% openEMS matlab interface
% -----------------------
% Thorsten Liebig <thorsten.liebig@gmx.de> (c) 2013
%
% See also InitCSX AddMetal AddMaterial AddExcitation calcPort
%% validate arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%check mesh
if ~isfield(CSX,'RectilinearGrid')
error 'mesh needs to be defined! Use DefineRectGrid() first!';
end
if (~isfield(CSX.RectilinearGrid,'XLines') || ~isfield(CSX.RectilinearGrid,'YLines') || ~isfield(CSX.RectilinearGrid,'ZLines'))
error 'mesh needs to be defined! Use DefineRectGrid() first!';
end
% check dir
dir = DirChar2Int(dir);
%set defaults
feed_shift = 0;
feed_R = inf; %(default is open, no resitance)
excite_amp = 0;
measplanepos = nan;
PortNamePrefix = '';
excite_args = {};
%% read optional arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n=1:2:numel(varargin)
if (strcmp(varargin{n},'FeedShift')==1);
feed_shift = varargin{n+1};
if (numel(feed_shift)>1)
error 'FeedShift must be a scalar value'
end
elseif (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
elseif (strcmp(varargin{n},'MeasPlaneShift')==1);
measplanepos = varargin{n+1};
if (numel(feed_shift)>1)
error 'MeasPlaneShift must be a scalar value'
end
elseif (strcmp(varargin{n},'ExciteAmp')==1);
excite_amp = varargin{n+1};
elseif (strcmpi(varargin{n},'PortNamePrefix'))
PortNamePrefix = varargin{n+1};
else
excite_args{end+1} = varargin{n};
excite_args{end+1} = varargin{n+1};
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine index (1, 2 or 3) of propagation (length of MSL)
idx_prop_n = dir + 1;
idx_prop_nP = mod((dir+1),3)+1;
idx_prop_nPP = mod((dir+2),3)+1;
% direction of propagation
if stop(idx_prop_n)-start(idx_prop_n) > 0
direction = +1;
else
direction = -1;
end
% create the metal for the coaxial line
CSX = AddCylinder( CSX, pec_name, prio, start, stop, r_i );
CSX = AddCylindricalShell( CSX, pec_name, prio, start, stop, 0.5*(r_o+r_os), r_os-r_o );
% create the material filling for the coaxial line
if (~isempty(materialname))
CSX = AddCylindricalShell( CSX, materialname, prio-1, start, stop, 0.5*(r_o+r_i), r_o-r_i );
end
if isnan(measplanepos)
measplanepos = (start(idx_prop_n)+stop(idx_prop_n))/2;
else
measplanepos = start(idx_prop_n)+direction*measplanepos;
end
% calculate position of the voltage probes
mesh{1} = sort(unique(CSX.RectilinearGrid.XLines));
mesh{2} = sort(unique(CSX.RectilinearGrid.YLines));
mesh{3} = sort(unique(CSX.RectilinearGrid.ZLines));
meshlines = interp1( mesh{idx_prop_n}, 1:numel(mesh{idx_prop_n}), measplanepos, 'nearest' );
meshlines = mesh{idx_prop_n}(meshlines-1:meshlines+1); % get three lines (approx. at center)
if direction == -1
meshlines = fliplr(meshlines);
end
v1_start(idx_prop_n) = meshlines(1);
v1_start(idx_prop_nP) = start(idx_prop_nP)+r_i;
v1_start(idx_prop_nPP) = start(idx_prop_nPP);
v1_stop = v1_start;
v1_stop(idx_prop_nP) = start(idx_prop_nP)+r_o;
v2_start = v1_start;
v2_stop = v1_stop;
v2_start(idx_prop_n) = meshlines(2);
v2_stop(idx_prop_n) = meshlines(2);
v3_start = v2_start;
v3_stop = v2_stop;
v3_start(idx_prop_n) = meshlines(3);
v3_stop(idx_prop_n) = meshlines(3);
% calculate position of the current probes
i1_start(idx_prop_n) = 0.5*(meshlines(1)+meshlines(2));
i1_start(idx_prop_nP) = start(idx_prop_nP)-r_i-0.1*(r_o-r_i);
i1_start(idx_prop_nPP) = start(idx_prop_nPP)-r_i-0.1*(r_o-r_i);
i1_stop = i1_start;
i1_stop(idx_prop_nP) = start(idx_prop_nP)+r_i+0.1*(r_o-r_i);
i1_stop(idx_prop_nPP) = start(idx_prop_nPP)+r_i+0.1*(r_o-r_i);
i2_start = i1_start;
i2_stop = i1_stop;
i2_start(idx_prop_n) = 0.5*(meshlines(2)+meshlines(3));
i2_stop(idx_prop_n) = 0.5*(meshlines(2)+meshlines(3));
% create the probes
port.U_filename{1} = [PortNamePrefix 'port_ut' num2str(portnr) 'A'];
weight = 1;
CSX = AddProbe( CSX, port.U_filename{1}, 0, 'weight', weight );
CSX = AddBox( CSX, port.U_filename{1}, prio, v1_start, v1_stop );
port.U_filename{2} = [PortNamePrefix 'port_ut' num2str(portnr) 'B'];
CSX = AddProbe( CSX, port.U_filename{2}, 0, 'weight', weight );
CSX = AddBox( CSX, port.U_filename{2}, prio, v2_start, v2_stop );
port.U_filename{3} = [PortNamePrefix 'port_ut' num2str(portnr) 'C'];
CSX = AddProbe( CSX, port.U_filename{3}, 0, 'weight', weight );
CSX = AddBox( CSX, port.U_filename{3}, prio, v3_start, v3_stop );
weight = direction;
port.I_filename{1} = [PortNamePrefix 'port_it' num2str(portnr) 'A'];
CSX = AddProbe( CSX, port.I_filename{1}, 1, 'weight', weight );
CSX = AddBox( CSX, port.I_filename{1}, prio, i1_start, i1_stop );
port.I_filename{2} = [PortNamePrefix 'port_it' num2str(portnr) 'B'];
CSX = AddProbe( CSX, port.I_filename{2}, 1,'weight', weight );
CSX = AddBox( CSX, port.I_filename{2}, prio, i2_start, i2_stop );
% create port structure
port.LengthScale = 1;
port.nr = portnr;
port.type = 'Coaxial';
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;
port.direction = direction;
port.excite = 0;
port.measplanepos = abs(v2_start(idx_prop_n) - start(idx_prop_n))*port.LengthScale;
port.r_i = r_i;
port.r_o = r_o;
% create excitation (if enabled) and port resistance
meshline = interp1( mesh{idx_prop_n}, 1:numel(mesh{idx_prop_n}), start(idx_prop_n) + feed_shift*direction, 'nearest' );
min_cell_prop = min(diff(mesh{idx_prop_n}));
ex_start = start;
ex_start(idx_prop_n) = mesh{idx_prop_n}(meshline) - 0.01*min_cell_prop;
ex_stop = ex_start;
ex_stop(idx_prop_n) = mesh{idx_prop_n}(meshline) + 0.01*min_cell_prop;
port.excite = 0;
if (excite_amp~=0)
dir_names={'x','y','z'};
nameX = ['(' dir_names{idx_prop_nP} '-' num2str(start(idx_prop_nP)) ')'];
nameY = ['(' dir_names{idx_prop_nPP} '-' num2str(start(idx_prop_nPP)) ')'];
func_Ex = [ nameX '/(' nameX '*' nameX '+' nameY '*' nameY ') * (sqrt(' nameX '*' nameX '+' nameY '*' nameY ')<' num2str(r_o) ') * (sqrt(' nameX '*' nameX '+' nameY '*' nameY ')>' num2str(r_i) ')'];
func_Ey = [ nameY '/(' nameX '*' nameX '+' nameY '*' nameY ') * (sqrt(' nameX '*' nameX '+' nameY '*' nameY ')<' num2str(r_o) ') * (sqrt(' nameX '*' nameX '+' nameY '*' nameY ')>' num2str(r_i) ')'];
func_E{idx_prop_n} = 0;
func_E{idx_prop_nP} = func_Ex;
func_E{idx_prop_nPP} = func_Ey;
port.excite = 1;
evec = [1 1 1];
evec(idx_prop_n) = 0;
CSX = AddExcitation( CSX, [PortNamePrefix 'port_excite_' num2str(portnr)], 0, evec, excite_args{:} );
CSX = SetExcitationWeight(CSX, [PortNamePrefix 'port_excite_' num2str(portnr)], func_E );
CSX = AddCylindricalShell(CSX,[PortNamePrefix 'port_excite_' num2str(portnr)],0 ,ex_start,ex_stop,0.5*(r_i+r_o),(r_o-r_i));
end
%% resitance at start of coaxial line
ex_start = start;
ex_stop = stop;
ex_stop(idx_prop_n) = ex_start(idx_prop_n);
if (feed_R > 0) && ~isinf(feed_R)
error 'feed_R not yet implemented'
elseif isinf(feed_R)
% do nothing --> open port
elseif feed_R == 0
%port "resistance" as metal
CSX = AddBox( CSX, pec_name, prio, ex_start, ex_stop );
CSX = AddCylindricalShell(CSX, pec_name, prio ,ex_start, ex_stop, 0.5*(r_i+r_o),(r_o-r_i));
else
error('openEMS:AddMSLPort','MSL port with resitance <= 0 it not possible');
end
end

View File

@ -46,7 +46,7 @@ if (iscell(port))
return;
end
if strcmpi(port.type,'MSL')
if (strcmpi(port.type,'MSL') || strcmpi(port.type,'Coaxial'))
port = calcTLPort( port, SimDir, f, varargin{:});
return
elseif strcmpi(port.type,'WaveGuide')

View File

@ -48,7 +48,7 @@ if (iscell(port))
return;
end
if (strcmpi(port.type,'MSL')~=1)
if ((strcmpi(port.type,'MSL')~=1) && (strcmpi(port.type,'Coaxial')~=1))
error('openEMS:calcTLPort','error, type is not a transmission line port');
end
@ -102,6 +102,10 @@ beta(real(beta) < 0) = -beta(real(beta) < 0); % determine correct sign (unlike t
% determine ZL
ZL = sqrt(Et .* dEt ./ (Ht .* dHt));
% if (strcmpi(port.type,'Coaxial'))
% port.ZL = Z0/2/pi/ref_index*log(port.r_o/port.r_i);
% end
% reference plane shift (lossless)
if ~isnan(ref_shift)
ref_shift = ref_shift * port.LengthScale;

View File

@ -50,35 +50,28 @@ end
%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FDTD = InitFDTD(numTS,1e-5);
FDTD = SetGaussExcite(FDTD,f0,f0);
BC = {'PEC','PEC','PEC','PEC','PEC','MUR'};
BC = {'PEC','PEC','PEC','PEC','MUR','MUR'};
if (use_pml>0)
BC = {'PEC','PEC','PEC','PEC','PEC','PML_8'};
BC = {'PEC','PEC','PEC','PEC','PML_8','PML_8'};
end
FDTD = SetBoundaryCond(FDTD,BC);
%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CSX = InitCSX();
mesh.x = -2.5*mesh_res(1)-coax_rad_aa : mesh_res(1) : coax_rad_aa+2.5*mesh_res(1);
mesh.x = -coax_rad_aa : mesh_res(1) : coax_rad_aa;
mesh.y = mesh.x;
mesh.z = 0 : mesh_res(3) : length;
mesh.z = SmoothMeshLines([0 length], mesh_res(3));
CSX = DefineRectGrid(CSX, unit, mesh);
%%% coax
CSX = AddMaterial(CSX,'copper');
CSX = SetMaterialProperty(CSX,'copper','Kappa',56e6);
start = [0, 0 , 0];stop = [0, 0 , length];
CSX = AddCylinder(CSX,'copper',0 ,start,stop,coax_rad_i);
CSX = AddCylindricalShell(CSX,'copper',0 ,start,stop,0.5*(coax_rad_aa+coax_rad_ai),(coax_rad_aa-coax_rad_ai));
CSX = AddMetal(CSX,'copper');
start = [0,0,0];
stop = [0,0,length/2];
[CSX,port{1}] = AddCoaxialPort( CSX, 10, 1, 'copper', '', start, stop, 'z', coax_rad_i, coax_rad_ai, coax_rad_aa, 'ExciteAmp', 1,'FeedShift', 10*mesh_res(1) );
%% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
start = [0,0,-0.1];
stop = [0,0, 0.1];
CSX = AddExcitation(CSX,'excite',0,[1 1 0]);
weight{1} = '(x)/(x*x+y*y)';
weight{2} = 'y/pow(rho,2)';
weight{3} = 0;
CSX = SetExcitationWeight(CSX, 'excite', weight );
CSX = AddCylindricalShell(CSX,'excite',0 ,start,stop,0.5*(coax_rad_i+coax_rad_ai),(coax_rad_ai-coax_rad_i));
start = [0,0,length/2];
stop = [0,0,length];
[CSX,port{2}] = AddCoaxialPort( CSX, 10, 2, 'copper', '', start, stop, 'z', coax_rad_i, coax_rad_ai, coax_rad_aa );
%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CSX = AddDump(CSX,'Et_','DumpMode',2);
@ -89,101 +82,39 @@ CSX = AddBox(CSX,'Et_',0 , start,stop);
CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',2);
CSX = AddBox(CSX,'Ht_',0,start,stop);
%% define voltage calc boxes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%voltage calc
CSX = AddProbe(CSX,'ut1',0);
start = [ coax_rad_i 0 mesh.z(10) ];
stop = [ coax_rad_ai 0 mesh.z(10) ];
CSX = AddBox(CSX,'ut1', 0 ,start,stop);
CSX = AddProbe(CSX,'ut2',0);
start = [ coax_rad_i 0 mesh.z(end-10)];
stop = [ coax_rad_ai 0 mesh.z(end-10)];
CSX = AddBox(CSX,'ut2', 0 ,start,stop);
%current calc, for each position there are two currents, which will get
%averaged to match the voltage position in between (!Yee grid!)
CSX = AddProbe(CSX,'it1a',1);
mid = 0.5*(coax_rad_i+coax_rad_ai);
start = [ -mid -mid mesh.z(9) ];
stop = [ mid mid mesh.z(9) ];
CSX = AddBox(CSX,'it1a', 0 ,start,stop);
CSX = AddProbe(CSX,'it1b',1);
start = [ -mid -mid mesh.z(10) ];
stop = [ mid mid mesh.z(10) ];
CSX = AddBox(CSX,'it1b', 0 ,start,stop);
CSX = AddProbe(CSX,'it2a',1);
start = [ -mid -mid mesh.z(end-11) ];
stop = [ mid mid mesh.z(end-11) ];
CSX = AddBox(CSX,'it2a', 0 ,start,stop);
CSX = AddProbe(CSX,'it2b',1);
start = [ -mid -mid mesh.z(end-10) ];
stop = [ mid mid mesh.z(end-10) ];
CSX = AddBox(CSX,'it2b', 0 ,start,stop);
%% Write openEMS
if (postprocessing_only==0)
WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
CSXGeomPlot([Sim_Path '/' Sim_CSX]);
RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts);
end
%%
freq = linspace(0,2*f0,201);
U = ReadUI({'ut1','ut2'},[Sim_Path '/'],freq);
I = ReadUI({'it1a','it1b','it2a','it2b'},[Sim_Path '/'],freq);
Exc = ReadUI('et',Sim_Path,freq);
port = calcPort(port, Sim_Path, freq);
%% plot voltages
%% plot s-parameter
figure
plot(U.TD{1}.t, U.TD{1}.val,'Linewidth',2);
hold on;
grid on;
plot(U.TD{2}.t, U.TD{2}.val,'r--','Linewidth',2);
xlabel('time (s)')
ylabel('voltage (V)')
legend('u_1(t)','u_2(t)')
%% calculate incoming and reflected voltages & currents
ZL_a = ones(size(freq))*Z0/2/pi/sqrt(epsR)*log(coax_rad_ai/coax_rad_i); %analytic line-impedance of a coax
uf1 = U.FD{1}.val./Exc.FD{1}.val;
uf2 = U.FD{2}.val./Exc.FD{1}.val;
if1 = 0.5*(I.FD{1}.val+I.FD{2}.val)./Exc.FD{1}.val;
if2 = 0.5*(I.FD{3}.val+I.FD{4}.val)./Exc.FD{1}.val;
uf1_inc = 0.5 * ( uf1 + if1 .* ZL_a );
if1_inc = 0.5 * ( if1 + uf1 ./ ZL_a );
uf2_inc = 0.5 * ( uf2 + if2 .* ZL_a );
if2_inc = 0.5 * ( if2 + uf2 ./ ZL_a );
uf1_ref = uf1 - uf1_inc;
if1_ref = if1 - if1_inc;
uf2_ref = uf2 - uf2_inc;
if2_ref = if2 - if2_inc;
% plot s-parameter
figure
s11 = uf1_ref./uf1_inc;
s21 = uf2_inc./uf1_inc;
s11 = port{1}.uf.ref./port{1}.uf.inc;
s21 = port{2}.uf.inc./port{1}.uf.inc;
plot(freq,20*log10(abs(s11)),'Linewidth',2);
hold on
grid on
plot(freq,20*log10(abs(s21)),'r--','Linewidth',2);
xlim([freq(1) freq(end)]);
xlabel('frequency (Hz)')
ylabel('s-para (dB)');
% ylim([-40 5]);
grid on;
hold on;
plot(freq,20*log10(abs(s21)),'r','Linewidth',2);
legend('s11','s21','Location','SouthEast');
% plot line-impedance comparison
%% plot line-impedance comparison
figure()
ZL = uf1./if1;
plot(freq,real(ZL),'Linewidth',2);
ZL_a = ones(size(freq))*Z0/2/pi/sqrt(epsR)*log(coax_rad_ai/coax_rad_i); %analytic line-impedance of a coax
ZL = port{2}.uf.tot./port{2}.if.tot;
plot(freq,real(port{1}.ZL),'Linewidth',2);
xlim([freq(1) freq(end)]);
xlabel('frequency (Hz)')
ylabel('line-impedance (\Omega)');
grid on;
hold on;
plot(freq,imag(ZL),'r--','Linewidth',2);
plot(freq,imag(port{1}.ZL),'r--','Linewidth',2);
plot(freq,ZL_a,'g-.','Linewidth',2);
legend('\Re\{ZL\}','\Im\{ZL\}','ZL-analytic','Location','Best');