From f290fd6db491237ebaa965d63698686600108862 Mon Sep 17 00:00:00 2001 From: Sebastian Held Date: Fri, 21 May 2010 11:41:33 +0200 Subject: [PATCH] matlab: added a microstrip port (with example) --- matlab/AddMSLPort.m | 161 ++++++++++++++++++++++++++++++ matlab/calcMSLPort.m | 90 +++++++++++++++++ matlab/examples/MSL2.m | 221 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 472 insertions(+) create mode 100644 matlab/AddMSLPort.m create mode 100644 matlab/calcMSLPort.m create mode 100644 matlab/examples/MSL2.m diff --git a/matlab/AddMSLPort.m b/matlab/AddMSLPort.m new file mode 100644 index 0000000..ca09ec4 --- /dev/null +++ b/matlab/AddMSLPort.m @@ -0,0 +1,161 @@ +function [CSX,port] = AddMSLPort( CSX, portnr, materialname, start, stop, dir, evec, excitename ) +% [CSX,port] = AddMSLPort( CSX, portnr, materialname, start, stop, dir, evec, excitename ) +% +% CSX: CSX-object created by InitCSX() +% portnr: (integer) number of the port +% materialname: property for the MSL (created by AddMetal() or AddMaterial()) +% 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]) +% evec: excitation vector, which defines the direction of the e-field (must be the same as used in AddExcitation()) +% 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) +% 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 +% +% See also InitCSX AddMetal AddMaterial AddExcitation calcMSLPort + +% check dir +if ~(dir(1) == dir(2) == 0) && ~(dir(1) == dir(3) == 0) && ~(dir(2) == dir(3) == 0) || (sum(dir) == 0) + 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) + error 'evec must have exactly one component ~= 0' +end +evec0 = evec ./ abs(sum(evec)); % evec0 is a unit vector + +% 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]; + +% determine index (1, 2 or 3) of width of MSL +idx_width = abs(cross(dir,evec0)) * [1;2;3]; + +% determine index (1, 2 or 3) of height +idx_height = abs(evec0) * [1;2;3]; + +% direction of propagation +if stop(idx_prop)-start(idx_prop) > 0 + direction = +1; +else + direction = -1; +end + +% create the metal/material for the MSL +MSL_start = start; +MSL_stop = stop; +MSL_stop(idx_height) = MSL_start(idx_height); +CSX = AddBox( CSX, materialname, 999, MSL_start, MSL_stop ); + +% FIXME +% openEMS v0.0.7 does not snap PEC + +% 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 = mesh{idx_prop}(meshlines-1:meshlines+1); % get three lines (approx. at center) +if direction == -1 + meshlines = fliplr(meshlines); +end +MSL_w2 = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), (nstart(idx_width)+nstop(idx_width))/2, 'nearest' ); +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_stop = v1_start; +v1_stop(idx_height) = nstart(idx_height); +v2_start = v1_start; +v2_stop = v1_stop; +v2_start(idx_prop) = meshlines(2); +v2_stop(idx_prop) = meshlines(2); +v3_start = v2_start; +v3_stop = v2_stop; +v3_start(idx_prop) = meshlines(3); +v3_stop(idx_prop) = meshlines(3); + +% calculate position of the current probes +idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstart(idx_width), 'nearest' ); +i1_start(idx_width) = mesh{idx_width}(idx) - diff(mesh{idx_width}(idx-1:idx))/2; +idx = interp1( mesh{idx_height}, 1:numel(mesh{idx_height}), start(idx_height), 'nearest' ); +i1_start(idx_height) = mesh{idx_height}(idx) - diff(mesh{idx_height}(idx-1:idx))/2; +i1_stop(idx_height) = mesh{idx_height}(idx) + diff(mesh{idx_height}(idx:idx+1))/2; +i1_start(idx_prop) = sum(meshlines(1:2))/2; +i1_stop(idx_prop) = i1_start(idx_prop); +idx = interp1( mesh{idx_width}, 1:numel(mesh{idx_width}), nstop(idx_width), 'nearest' ); +i1_stop(idx_width) = mesh{idx_width}(idx) + diff(mesh{idx_width}(idx:idx+1))/2; +i2_start = i1_start; +i2_stop = i1_stop; +i2_start(idx_prop) = sum(meshlines(2:3))/2; +i2_stop(idx_prop) = i2_start(idx_prop); + +% create the probes +name = ['port_ut' num2str(portnr) 'A']; +CSX = AddProbe( CSX, name, 0 ); +CSX = AddBox( CSX, name, 999, v1_start, v1_stop ); +name = ['port_ut' num2str(portnr) 'B']; +CSX = AddProbe( CSX, name, 0 ); +CSX = AddBox( CSX, name, 999, v2_start, v2_stop ); +name = ['port_ut' num2str(portnr) 'C']; +CSX = AddProbe( CSX, name, 0 ); +CSX = AddBox( CSX, name, 999, v3_start, v3_stop ); +name = ['port_it' num2str(portnr) 'A']; +CSX = AddProbe( CSX, name, 1 ); +CSX = AddBox( CSX, name, 999, i1_start, i1_stop ); +name = ['port_it' num2str(portnr) 'B']; +CSX = AddProbe( CSX, name, 1 ); +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.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.dir = dir; +port.evec = evec; +port.idx_prop = idx_prop; +port.idx_width = idx_width; +port.idx_height = idx_height; +port.excite = 0; + +% create excitation +if nargin >= 8 + % 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); + 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 = AddBox( CSX, excitename, 999, ex_start, ex_stop ); +end diff --git a/matlab/calcMSLPort.m b/matlab/calcMSLPort.m new file mode 100644 index 0000000..f78d84f --- /dev/null +++ b/matlab/calcMSLPort.m @@ -0,0 +1,90 @@ +function [S11,beta,ZL] = calcMSLPort( portstruct, SimDir, f, ref_shift ) +% [S11,beta,ZL] = calcMSLPort( portstruct, SimDir, [f], [ref_shift] ) +% +% portstruct: 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) +% +% 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 +% +% See also AddMSLPort + +% check +if portstruct.v_delta(1) ~= portstruct.v_delta(2) + warning( 'mesh is not equidistant; expect degraded accuracy' ); +end + +% read time domain data +filename = ['/port_ut' num2str(portstruct.nr)]; +U = ReadUI( {[filename 'A'],[filename 'B'],[filename 'C']}, SimDir ); +filename = ['/port_it' num2str(portstruct.nr)]; +I = ReadUI( {[filename 'A'],[filename 'B']}, SimDir ); + +if (nargin > 2) && ~isempty(f) + % freq vector given: use DFT + for n=1:numel(U.FD) + U.FD{n}.f = f; + U.FD{n}.val = DFT_time2freq( U.TD{n}.t, U.TD{n}.val, f ); + end + for n=1:numel(I.FD) + I.FD{n}.f = f; + I.FD{n}.val = DFT_time2freq( I.TD{n}.t, I.TD{n}.val, f ); + end +end + +delta_t = I.TD{1}.t(1) - U.TD{1}.t(1); +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); +Ht = (I.FD{1}.val + I.FD{2}.val)/2; % space averaging: Ht is now defined at the same pos as Et +Ht = Ht .* exp( -1i*2*pi*f * delta_t/2 ); % compensate time shift of Ht with respect to Et +dHt = (I.FD{2}.val - I.FD{1}.val) / (abs(portstruct.i_delta(1)) * portstruct.drawingunit); +dHt = dHt .* exp( -1i*2*pi*f * delta_t/2 ); % compensate time shift + +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 +if (nargin > 3) + % renormalize the shift to the measurement plane + if (portstruct.stop(portstruct.idx_prop) - portstruct.start(portstruct.idx_prop) > 0) + dir = +1; + else + dir = -1; + end + ref_shift = ref_shift - dir*(portstruct.v2_start(portstruct.idx_prop) - portstruct.start(portstruct.idx_prop)); + ref_shift = ref_shift * portstruct.drawingunit; + S11 = S11 .* exp(2i*real(beta)*ref_shift); + S11_corrected = S11_corrected .* exp(2i*real(beta)*ref_shift); +end diff --git a/matlab/examples/MSL2.m b/matlab/examples/MSL2.m new file mode 100644 index 0000000..fe4daaa --- /dev/null +++ b/matlab/examples/MSL2.m @@ -0,0 +1,221 @@ +% +% microstrip line example +% +% this example shows how to use a MSL port +% +% The MSL is excited at the center of the computational volume. The +% boundary at xmin is an absorbing boundary (Mur) and at xmax an electric +% wall. The reflection coefficient at this wall is S11 = -1. +% + + +close all +clear +clc + +physical_constants + + +postprocessing_only = 0; + + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +drawingunits = 1e-6; % specify everything in um +MSL_length = 10000; +MSL_width = 1000; +substrate_thickness = 254; + +mesh_res = [200 0 0]; +max_timesteps = 20000; +min_decrement = 1e-6; +f_max = 8e9; + +%% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitFDTD( max_timesteps, min_decrement, 'OverSampling', 10 ); +FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2 ); +BC = [2 0 0 0 0 1]; +FDTD = SetBoundaryCond( FDTD, BC ); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX(); +mesh.x = -MSL_length : mesh_res(1) : MSL_length; +mesh.y = linspace(-MSL_width/2,MSL_width/2,10); % discretize the width of the MSL with 10 cells +temp1 = linspace(-4*MSL_width,mesh.y(1),20); +temp2 = linspace(mesh.y(end),4*MSL_width,20); +mesh.y = [temp1(1:end-1), mesh.y, temp2(2:end)]; % add coarser discretization +mesh.z = linspace(0,substrate_thickness,5); % discretize the substrate with 5 cells +temp1 = linspace(substrate_thickness,2*substrate_thickness,5); +mesh.z = [mesh.z temp1(2:end)]; % add same space above the strip +temp1 = linspace(2*substrate_thickness,5*substrate_thickness,10); +mesh.z = [mesh.z temp1(2:end)]; % coarser discretization +CSX = DefineRectGrid( CSX, drawingunits, mesh ); + +%% Material definitions +CSX = AddMetal( CSX, 'PEC' ); +CSX = AddMaterial( CSX, 'RO4350B' ); + +%% substrate +CSX = SetMaterialProperty( CSX, 'RO4350B', 'Epsilon', 3.66 ); +start = [mesh.x(1), mesh.y(1), 0]; +stop = [mesh.x(end), mesh.y(end), substrate_thickness]; +CSX = AddBox( CSX, 'RO4350B', 0, start, stop ); + +%% MSL port +CSX = AddExcitation( CSX, 'excite', 0, [0 0 1]); +portstart = [ 0, -MSL_width/2, substrate_thickness]; +portstop = [ MSL_length, MSL_width/2, 0]; +[CSX,portstruct] = AddMSLPort( CSX, 1, 'PEC', portstart, portstop, [1 0 0], [0 0 1], 'excite' ); + +%% MSL +start = [-MSL_length, -MSL_width/2, substrate_thickness]; +stop = [ 0, MSL_width/2, substrate_thickness]; +CSX = AddBox( CSX, 'PEC', 0, start, stop ); + +%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = AddDump(CSX,'Et_','DumpType',0,'DumpMode',0); +start = [mesh.x(1) , mesh.y(1), substrate_thickness/2]; +stop = [mesh.x(end), mesh.y(end), substrate_thickness/2]; +CSX = AddBox(CSX,'Et_',0 , start,stop); + +CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',0); +CSX = AddBox(CSX,'Ht_',0,start,stop); + + +%% define openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +openEMS_opts = ''; +openEMS_opts = [openEMS_opts ' --disable-dumps']; +% openEMS_opts = [openEMS_opts ' --debug-material']; +% openEMS_opts = [openEMS_opts ' --debug-operator']; +% openEMS_opts = [openEMS_opts ' --debug-boxes']; +% openEMS_opts = [openEMS_opts ' --engine=sse-compressed']; +% openEMS_opts = [openEMS_opts ' --engine=multithreaded']; +openEMS_opts = [openEMS_opts ' --engine=fastest']; + +Sim_Path = 'tmp'; +Sim_CSX = 'MSL2.xml'; + +if ~postprocessing_only + rmdir(Sim_Path,'s'); +end +mkdir(Sim_Path); + +%% Write openEMS compatible xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); + +%% cd to working dir and run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +savePath = pwd; +cd(Sim_Path); %cd to working dir +args = [Sim_CSX ' ' openEMS_opts]; +if ~postprocessing_only + invoke_openEMS(args); +end +cd(savePath); + + + +%% postproc & do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +U = ReadUI({'port_ut1A','port_ut1B','et'},'tmp/'); +I = ReadUI({'port_it1A','port_it1B'},'tmp/'); +delta_t_2 = I.TD{1}.t(1) - U.TD{1}.t(1); % half time-step (s) + +% create finer frequency resolution +f = linspace( 0, f_max, 1601 ); +for n=1:numel(U.FD) + U.FD{n}.f = f; + U.FD{n}.val = DFT_time2freq( U.TD{n}.t, U.TD{n}.val, f ); +end +for n=1:numel(I.FD) + I.FD{n}.f = f; + I.FD{n}.val = DFT_time2freq( I.TD{n}.t, I.TD{n}.val, f ); + I.FD{n}.val = I.FD{n}.val .* exp(-1i*2*pi*I.FD{n}.f*delta_t_2); % compensate half time-step advance of H-field +end + +% interpolate et to the time spacing of the voltage probes +et = interp1( U.TD{3}.t, U.TD{3}.val, U.TD{1}.t ); + +f = U.FD{1}.f; + +% Z = (U.FD{1}.val+U.FD{2}.val)/2 ./ I.FD{1}.val; +% plot( f*1e-9, [real(Z);imag(Z)],'Linewidth',2); +% xlabel('frequency (GHz)'); +% ylabel('impedance (Ohm)'); +% grid on; +% legend( {'real','imaginary'}, 'location', 'northwest' ) +% title( 'line impedance (will fail in case of reflections!)' ); + +% figure +% plotyy(U.TD{1}.t/1e-6,[U.TD{1}.val;U.TD{2}.val],U.TD{1}.t/1e-6,et); +% xlabel('time (us)'); +% ylabel('amplitude (V)'); +% grid on; +% title( 'Time domain voltage probes and excitation signal' ); +% +% figure +% plot(I.TD{1}.t/1e-6,[I.TD{1}.val;I.TD{2}.val]); +% xlabel('time (us)'); +% ylabel('amplitude (A)'); +% grid on; +% title( 'Time domain current probes' ); + + +%% port analysis +[S11,beta,ZL] = calcMSLPort( portstruct, Sim_Path, f ); + +figure +plot( sin(0:0.01:2*pi), cos(0:0.01:2*pi), 'Color', [.7 .7 .7] ); +hold on +plot( 0.5+0.5*sin(0:0.01:2*pi), 0.5*cos(0:0.01:2*pi), 'Color', [.7 .7 .7] ); +plot( [-1 1], [0 0], 'Color', [.7 .7 .7] ); +plot( S11, 'k' ); +plot( real(S11(1)), imag(S11(1)), '*r' ); +axis equal +title( 'Reflection coefficient S11 at the measurement plane' ); + +figure +plot( f/1e9, [real(S11);imag(S11)], 'Linewidth',2 ); +legend( {'Re(S11)', 'Im(S11)'} ); +ylabel( 'amplitude' ); +xlabel( 'frequency (GHz)' ); +title( 'Reflection coefficient S11 at the measurement plane' ); + +figure +plotyy( f/1e9, 20*log10(abs(S11)), f/1e9, angle(S11)/pi*180 ); +legend( {'abs(S11)', 'angle(S11)'} ); +xlabel( 'frequency (GHz)' ); +title( 'Reflection coefficient S11 at the measurement plane' ); + +figure +plot( f/1e9, [real(beta);imag(beta)], 'Linewidth',2 ); +legend( 'Re(beta)', 'Im(beta)' ); +ylabel( 'propagation constant beta (1/m)' ); +xlabel( 'frequency (GHz)' ); +title( 'Propagation constant of the MSL' ); + +figure +plot( f/1e9, [real(ZL);imag(ZL)], 'Linewidth',2); +xlabel('frequency (GHz)'); +ylabel('impedance (Ohm)'); +grid on; +legend( {'real','imaginary'}, 'location', 'northeast' ) +title( 'Characteristic line impedance ZL' ); +ylim( [-2*mean(real(ZL)) 2*mean(real(ZL))] ); + +% reference plane shift (to the end of the port) +ref_shift = abs(portstop(1) - portstart(1)); +[S11,beta,ZL] = calcMSLPort( portstruct, Sim_Path, f, ref_shift ); + +figure +plotyy( f/1e9, 20*log10(abs(S11)), f/1e9, angle(S11)/pi*180 ); +legend( {'abs(S11)', 'angle(S11)'} ); +xlabel( 'frequency (GHz)' ); +title( 'Reflection coefficient S11 at the reference plane (at the electric wall)' ); + +figure +plot( sin(0:0.01:2*pi), cos(0:0.01:2*pi), 'Color', [.7 .7 .7] ); +hold on +plot( 0.5+0.5*sin(0:0.01:2*pi), 0.5*cos(0:0.01:2*pi), 'Color', [.7 .7 .7] ); +plot( [-1 1], [0 0], 'Color', [.7 .7 .7] ); +plot( S11, 'k' ); +plot( real(S11(1)), imag(S11(1)), '*r' ); +axis equal +title( 'Reflection coefficient S11 at the reference plane (at the electric wall)' );