diff --git a/matlab/AddCoaxialPort.m b/matlab/AddCoaxialPort.m new file mode 100644 index 0000000..45c0d2b --- /dev/null +++ b/matlab/AddCoaxialPort.m @@ -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 (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 diff --git a/matlab/calcPort.m b/matlab/calcPort.m index 2033f69..14593c5 100644 --- a/matlab/calcPort.m +++ b/matlab/calcPort.m @@ -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') diff --git a/matlab/calcTLPort.m b/matlab/calcTLPort.m index 6df2d1c..8c10174 100644 --- a/matlab/calcTLPort.m +++ b/matlab/calcTLPort.m @@ -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; diff --git a/matlab/examples/waveguide/Coax.m b/matlab/examples/waveguide/Coax.m index d7ce8a7..3db4d6a 100644 --- a/matlab/examples/waveguide/Coax.m +++ b/matlab/examples/waveguide/Coax.m @@ -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');