diff --git a/matlab/examples/resistance_sheet.m b/matlab/examples/resistance_sheet.m new file mode 100644 index 0000000..b63c00e --- /dev/null +++ b/matlab/examples/resistance_sheet.m @@ -0,0 +1,207 @@ +% +% resistance "sheet" example +% +% this example calculates the reflection coefficient of a sheet resistance +% at the end of a parallel plate wave guide +% +% play around with the R and epr values +% + +close all +clear +clc + +physical_constants + + +postprocessing_only = 0; + + + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +epr = 1; % relative permittivity of the material inside the parallel plate waveguide + +% define the resistance +R = sqrt(MUE0/(EPS0*epr)); % matched load (no reflections) (vacuum: approx. 377 Ohm) +% R = 1e-10; % short circuit (reflection coefficient = -1) +% R = 1e10; % open circuit (reflection coefficient = 1) + + +drawingunit = 1e-6; % specify everything in um +length = 10000; +mesh_res = [200 200 200]; +max_timesteps = 100000; +min_decrement = 1e-6; +f_max = 1e9; + +%% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitFDTD( max_timesteps, min_decrement ); +FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2 ); +BC = [1 2 1 1 0 0]; % 0:PEC 1:PMC 2:MUR-ABC +FDTD = SetBoundaryCond( FDTD, BC ); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX(); +mesh.x = 0 : mesh_res(1) : length; +mesh.y = -2*mesh_res(2) : mesh_res(2) : 2*mesh_res(2); +mesh.z = 0 : mesh_res(3) : 4*mesh_res(3); +CSX = DefineRectGrid( CSX, drawingunit, mesh ); + +%% measurement plane & reference plane +meas_plane_xidx = interp1( mesh.x, 1:numel(mesh.x), length*1/3, 'nearest' ); +ref_plane_xidx = 3; + +%% fill the parallel plate waveguide with material +CSX = AddMaterial( CSX, 'm1' ); +CSX = SetMaterialProperty( CSX, 'm1', 'Epsilon', epr ); +start = [mesh.x(1), mesh.y(1), mesh.z(1)]; +stop = [mesh.x(end), mesh.y(end), mesh.z(end)]; +CSX = AddBox( CSX, 'm1', -1, start, stop ); + +%% excitation +CSX = AddExcitation( CSX, 'excitation1', 0, [0 0 1]); +idx = interp1( mesh.x, 1:numel(mesh.x), length*2/3, 'nearest' ); +start = [mesh.x(idx), mesh.y(1), mesh.z(1)]; +stop = [mesh.x(idx), mesh.y(end), mesh.z(end)]; +CSX = AddBox( CSX, 'excitation1', 0, start, stop ); + +%% define the sheet resistance +start = [mesh.x(ref_plane_xidx-1), mesh.y(1), mesh.z(1)]; +stop = [mesh.x(ref_plane_xidx), mesh.y(end), mesh.z(end)]; +l = abs(mesh.z(end) - mesh.z(1)) * drawingunit; % length of the "sheet" +A = abs(start(1) - stop(1)) * abs(mesh.y(end) - mesh.y(1)) * drawingunit^2; % area of the "sheet" +kappa = l/A / R; % [kappa] = S/m +CSX = AddMaterial( CSX, 'sheet_resistance' ); +CSX = SetMaterialProperty( CSX, 'sheet_resistance', 'Kappa', kappa ); +CSX = AddBox( CSX, 'sheet_resistance', 0, start, stop ); + +%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = AddDump( CSX, 'Et_', 'DumpMode', 2 ); +start = [mesh.x(1), mesh.y(1), mesh.z(3)]; +stop = [mesh.x(end), mesh.y(end), mesh.z(3)]; +CSX = AddBox( CSX, 'Et_', 0, start, stop ); + +CSX = AddDump( CSX, 'Ht_', 'DumpType', 1, 'DumpMode', 2 ); +CSX = AddBox( CSX, 'Ht_', 0, start, stop ); + +% hdf5 file +CSX = AddDump( CSX, 'E', 'DumpType', 0, 'DumpMode', 2, 'FileType', 1 ); +start = [mesh.x(meas_plane_xidx), mesh.y(3), mesh.z(1)]; +stop = [mesh.x(meas_plane_xidx), mesh.y(3), mesh.z(end)]; +CSX = AddBox( CSX, 'E', 0, start, stop ); + +% hdf5 file +CSX = AddDump( CSX, 'H', 'DumpType', 1, 'DumpMode', 2, 'FileType', 1 ); +start = [mesh.x(meas_plane_xidx), mesh.y(1), mesh.z(3)]; +stop = [mesh.x(meas_plane_xidx), mesh.y(end), mesh.z(3)]; +CSX = AddBox( CSX, 'H', 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 ' --showProbeDiscretization']; +openEMS_opts = [openEMS_opts ' --engine=fastest']; + +Sim_Path = 'tmp'; +Sim_CSX = 'tmp.xml'; + +if ~postprocessing_only + [~,~,~] = rmdir(Sim_Path,'s'); + [~,~,~] = mkdir(Sim_Path); +end + +%% Write openEMS compatible xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + +%% cd to working dir and run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if ~postprocessing_only + savePath = pwd; + cd(Sim_Path); %cd to working dir + args = [Sim_CSX ' ' openEMS_opts]; + invoke_openEMS(args); + cd(savePath) +end + + +%% postproc & do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% E_coords = ReadHDF5Mesh( [Sim_Path '/E.h5'] ); +% H_coords = ReadHDF5Mesh( [Sim_Path '/H.h5'] ); +E = ReadHDF5FieldData( [Sim_Path '/E.h5'] ); +H = ReadHDF5FieldData( [Sim_Path '/H.h5'] ); +E_val = cellfun( @(x) squeeze(x(1,1,:,3)), E.values, 'UniformOutput', false ); +H_val = cellfun( @(x) squeeze(x(1,:,1,2)), H.values, 'UniformOutput', false ); +E_val = cell2mat(E_val); +H_val = cell2mat(H_val.'); + +% pick center point +Et = E_val(3,:); +Ht = H_val(:,3).'; + +delta_t_2 = H.time(1) - E.time(1); % half time-step (s) + +% create finer frequency resolution +f = linspace( 0, f_max, 201 ); +Ef = DFT_time2freq( E.time, Et, f ); +Hf = DFT_time2freq( H.time, Ht, f ); +Hf = Hf .* exp(-1i*2*pi*f*delta_t_2); % compensate half time-step advance of H-field + +% H is now time interpolated, but the position is not corrected with +% respect to E + +% figure +% plot( E.time/1e-6, Et ); +% xlabel('time (us)'); +% ylabel('amplitude (V)'); +% grid on; +% title( 'Time domain voltage probe' ); +% +% figure +% plot( H.time/1e-6, Ht ); +% xlabel('time (us)'); +% ylabel('amplitude (A)'); +% grid on; +% title( 'Time domain current probe' ); + + +Z0 = sqrt(MUE0/(EPS0*epr)); % line impedance +Z = Ef ./ Hf; % impedance at measurement plane +gamma = (Z - Z0) ./ (Z + Z0); + +% reference plane shift +beta = 2*pi*f * sqrt(MUE0*(EPS0*epr)); % TEM wave +meas_plane_x = mesh.x(meas_plane_xidx); +ref_plane_x = mesh.x(ref_plane_xidx); +gamma_refplane = gamma .* exp(2i*beta* (meas_plane_x-ref_plane_x)*drawingunit); +Z_refplane = Z0 * (1+gamma_refplane)./(1-gamma_refplane); + +% smith chart +figure +if exist( 'smith', 'file' ) + % smith chart + % www.ece.rutgers.edu/~orfanidi/ewa + % or cmt toolbox from git.ate.uni-duisburg.de + smith +else + % poor man smith chart + plot( sin(0:0.01:2*pi), cos(0:0.01:2*pi), 'Color', [.7 .7 .7] ); + hold on +% plot( 0.25+0.75*sin(0:0.01:2*pi), 0.75*cos(0:0.01:2*pi), 'Color', [.7 .7 .7] ); + plot( 0.5+0.5*sin(0:0.01:2*pi), 0.5*cos(0:0.01:2*pi), 'Color', [.7 .7 .7] ); +% plot( 0.75+0.25*sin(0:0.01:2*pi), 0.25*cos(0:0.01:2*pi), 'Color', [.7 .7 .7] ); + plot( [-1 1], [0 0], 'Color', [.7 .7 .7] ); + axis equal +end +plot( real(gamma_refplane), imag(gamma_refplane), 'r*' ); +% plot( real(gamma), imag(gamma), 'k*' ); +title( 'reflection coefficient S11 at reference plane' ) + +figure +plot( f/1e9, [real(Z_refplane);imag(Z_refplane)],'Linewidth',2); +xlabel('frequency (GHz)'); +ylabel('impedance (Ohm)'); +grid on; +title( 'Impedance at reference plane' ); +legend( {'real','imag'} );