diff options
Diffstat (limited to 'openEMS/matlab/examples')
27 files changed, 4324 insertions, 0 deletions
diff --git a/openEMS/matlab/examples/__deprecated__/MSL2.m b/openEMS/matlab/examples/__deprecated__/MSL2.m new file mode 100644 index 0000000..31a2600 --- /dev/null +++ b/openEMS/matlab/examples/__deprecated__/MSL2.m @@ -0,0 +1,254 @@ +% +% EXAMPLE / microstrip / MSL2 +% +% This example shows how to use the 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. +% Direction of propagation is x. +% +% This example demonstrates: +% - simple microstrip geometry (made of PEC) +% - MSL port +% - MSL analysis +% +% You may modify the PEC boundary condition at xmax to become a MUR +% boundary. This resembles a matched microstrip line. +% +% Tested with +% - Matlab 2009b +% - Octave 3.3.52 +% - openEMS v0.0.14 +% +% (C) 2010 Sebastian Held <sebastian.held@uni-due.de> + +close all +clear +clc + +%% switches +postproc_only = 0; + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +physical_constants; +unit = 1e-6; % specify everything in um +MSL_length = 10000; +MSL_width = 1000; +substrate_thickness = 254; +substrate_epr = 3.66; + +% mesh_res = [200 0 0]; + +%% prepare simulation folder +Sim_Path = 'tmp'; +Sim_CSX = 'msl2.xml'; +if ~postproc_only + [status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory + [status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder +end + +%% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%% +max_timesteps = 20000; +min_decrement = 1e-6; +f_max = 7e9; +FDTD = InitFDTD( max_timesteps, min_decrement, 'OverSampling', 10 ); +FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2 ); +BC = {'MUR' 'MUR' 'PEC' 'PEC' 'PEC' 'PMC'}; +FDTD = SetBoundaryCond( FDTD, BC ); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX(); +resolution = c0/(f_max*sqrt(substrate_epr))/unit /50; % resolution of lambda/50 +mesh.x = SmoothMeshLines( [-MSL_length MSL_length], resolution ); +mesh.y = SmoothMeshLines( [-4*MSL_width -MSL_width/2 MSL_width/2 4*MSL_width], resolution ); +mesh.z = SmoothMeshLines( [linspace(0,substrate_thickness,5) 10*substrate_thickness], resolution ); +CSX = DefineRectGrid( CSX, unit, mesh ); + +%% substrate +CSX = AddMaterial( CSX, 'RO4350B' ); +CSX = SetMaterialProperty( CSX, 'RO4350B', 'Epsilon', substrate_epr ); +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 = AddMetal( CSX, 'PEC' ); +portstart = [ 0, -MSL_width/2, substrate_thickness]; +portstop = [ MSL_length, MSL_width/2, 0]; +[CSX,portstruct] = AddMSLPort( CSX, 999, 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', 999, start, stop ); % priority needs to be higher than + +%% define dump boxes +start = [mesh.x(1), mesh.y(1), substrate_thickness/2]; +stop = [mesh.x(end), mesh.y(end), substrate_thickness/2]; +CSX = AddDump( CSX, 'Et_', 'DumpType', 0,'DumpMode', 2 ); % cell interpolated +CSX = AddBox( CSX, 'Et_', 0, start, stop ); +CSX = AddDump( CSX, 'Ht_', 'DumpType', 1,'DumpMode', 2 ); % cell interpolated +CSX = AddBox( CSX, 'Ht_', 0, start, stop ); + +%% write openEMS compatible xml-file +WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); + +%% show the structure +if ~postproc_only + CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); +end + +%% run openEMS +openEMS_opts = ''; +openEMS_opts = [openEMS_opts ' --engine=fastest']; +% openEMS_opts = [openEMS_opts ' --debug-material']; +% openEMS_opts = [openEMS_opts ' --debug-boxes']; +% openEMS_opts = [openEMS_opts ' --debug-PEC']; +if ~postproc_only + RunOpenEMS( Sim_Path, Sim_CSX, openEMS_opts ); +end + + +%% postprocess +f = linspace( 1e6, f_max, 1601 ); +U = ReadUI( {'port_ut1A','port_ut1B','port_ut1C','et'}, 'tmp/', f ); +I = ReadUI( {'port_it1A','port_it1B'}, 'tmp/', 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 +ax = plotyy( U.TD{1}.t/1e-6, [U.TD{1}.val;U.TD{2}.val;U.TD{3}.val], U.TD{4}.t/1e-6, U.TD{4}.val ); +xlabel( 'time (us)' ); +ylabel( 'amplitude (V)' ); +grid on +title( 'Time domain voltage probes and excitation signal' ); +legend( {'ut1A','ut1B','ut1C','excitation'} ); +% now make the y-axis symmetric to y=0 (align zeros of y1 and y2) +y1 = ylim(ax(1)); +y2 = ylim(ax(2)); +ylim( ax(1), [-max(abs(y1)) max(abs(y1))] ); +ylim( ax(2), [-max(abs(y2)) max(abs(y2))] ); + +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' ); +legend( {'it1A','it1B'} ); + +figure +ax = plotyy( U.FD{1}.f/1e9, abs([U.FD{1}.val;U.FD{2}.val;U.FD{3}.val]), U.FD{1}.f/1e9, angle([U.FD{1}.val;U.FD{2}.val;U.FD{3}.val])/pi*180 ); +xlabel( 'frequency (GHz)' ); +ylabel( ax(1), 'amplitude (A)' ); +ylabel( ax(2), 'phase (deg)' ); +grid on +title( 'Frequency domain voltage probes' ); +legend( {'abs(uf1A)','abs(uf1B)','abs(uf1C)','angle(uf1A)','angle(uf1B)','angle(uf1C)'} ); + +figure +ax = plotyy( I.FD{1}.f/1e9, abs([I.FD{1}.val;I.FD{2}.val]), I.FD{1}.f/1e9, angle([I.FD{1}.val;I.FD{2}.val])/pi*180 ); +xlabel( 'frequency (GHz)' ); +ylabel( ax(1), 'amplitude (A)' ); +ylabel( ax(2), 'phase (deg)' ); +grid on +title( 'Frequency domain current probes' ); +legend( {'abs(if1A)','abs(if1B)','angle(if1A)','angle(if1B)'} ); + +%% port analysis +[U,I,beta,ZL] = calcPort( portstruct, Sim_Path, f ); +%% attention! the reflection coefficient S11 is normalized to ZL! + +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( 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] ); +Z = vi.FD.v.val ./ vi.FD.i.val; +S11_ = (Z-ZL) ./ (Z+ZL); +plot( S11_, 'k' ); +plot( real(S11_(1)), imag(S11_(1)), '*r' ); +axis equal +title( {'Reflection coefficient S11 at the measurement plane' 'calculated from voltages and currents'} ); + +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( {'|S11|', 'angle(S11)'} ); +xlabel( 'frequency (GHz)' ); +ylabel( '|S11| (dB)' ); +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' ); + +%% reference plane shift (to the end of the port) +ref_shift = abs(portstop(1) - portstart(1)); +[U, I,beta,ZL] = calcPort( portstruct, Sim_Path, f ); +%% + +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)' ); + +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] ); +Z = vi.FD.v.val_shifted ./ vi.FD.i.val_shifted; +S11_ = (Z-ZL) ./ (Z+ZL); +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)' 'calculated from shifted voltages and currents'} ); + +%% visualize electric and magnetic fields +% you will find vtk dump files in the simulation folder (tmp/) +% use paraview to visualize them diff --git a/openEMS/matlab/examples/antennas/Bi_Quad_Antenna.m b/openEMS/matlab/examples/antennas/Bi_Quad_Antenna.m new file mode 100644 index 0000000..80ae97f --- /dev/null +++ b/openEMS/matlab/examples/antennas/Bi_Quad_Antenna.m @@ -0,0 +1,139 @@ +% +% Tutorials / bi-quad antenna +% +% Tested with +% - Octave 3.8.1 +% - openEMS v0.0.32 +% +% (C) 2011-2014 Thorsten Liebig <thorsten.liebig@uni-due.de> + +close all +clear +clc + +%% setup the simulation +physical_constants; +unit = 1e-3; % all length in mm + +quad_size = 110; +port_length = 10; +quad_mesh = 5; + +Feed_R = 75; + +% size of the simulation box +SimBox = [800 800 400]; + +% frequency range of interest +f_start = 400e6; +f_stop = 1000e6; + +% frequency of interest +f0 = 700e6; +freq = linspace(f_start,f_stop,201); + +%% setup FDTD parameter & excitation function +FDTD = InitFDTD( 'endCriteria', 1e-4 ); +FDTD = SetGaussExcite(FDTD,0.5*(f_start+f_stop),0.5*(f_stop-f_start)); +BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'}; % boundary conditions +FDTD = SetBoundaryCond( FDTD, BC ); + +%% setup CSXCAD geometry & mesh +CSX = InitCSX(); + +%create fixed lines for the antenna outline and port +mesh.x = [-quad_size*sqrt(2) -quad_size/sqrt(2) 0 quad_size/sqrt(2) quad_size*sqrt(2)]; +mesh.y = [-quad_size/sqrt(2) -port_length/2 0 port_length/2 quad_size/sqrt(2)]; +mesh.z = [0]; + +mesh = SmoothMesh(mesh, quad_mesh, 1.3); + +% add air box +mesh.x = [mesh.x -SimBox(1)/2 SimBox(1)/2]; +mesh.y = [mesh.y -SimBox(2)/2 SimBox(2)/2]; +mesh.z = [-SimBox(3)/2 0 SimBox(3)/2]; + +max_res = c0 / (f_stop) / unit / 20; % cell size: lambda/20 +mesh = SmoothMesh(mesh, max_res, 1.4); + +CSX = DefineRectGrid( CSX, unit, mesh ); + +%% create bi-quad +points(1,1) = 0; +points(2,1) = port_length/2; +points(3,1) = 0; +points(1,end+1) = quad_size/sqrt(2); +points(2,end) = quad_size/sqrt(2); +points(1,end+1) = quad_size*sqrt(2); +points(2,end) = 0; +points(1,end+1) = quad_size/sqrt(2); +points(2,end) = -quad_size/sqrt(2); +points(1,end+1) = 0; +points(2,end) = -port_length/2; +points(1,end+1) = -quad_size/sqrt(2); +points(2,end) = -quad_size/sqrt(2); +points(1,end+1) = -quad_size*sqrt(2); +points(2,end) = 0; +points(1,end+1) = -quad_size/sqrt(2); +points(2,end) = quad_size/sqrt(2); +points(1,end+1) = 0; +points(2,end) = port_length/2; + +% create a thin metal wire... +CSX = AddMetal(CSX,'metal'); %create PEC with propName 'metal' +CSX = AddCurve(CSX,'metal',10, points); + +%% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +start = [0 -port_length/2 0]; +stop = [0 port_length/2 0]; +[CSX port] = AddLumpedPort(CSX,10,0,Feed_R,start,stop,[0 1 0], true); + +%% nf2ff calc +start = [mesh.x(9) mesh.y(9) mesh.z(9)]; +stop = [mesh.x(end-8) mesh.y(end-8) mesh.z(end-8)]; +[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop); + +%% prepare simulation folder +Sim_Path = 'tmp'; +Sim_CSX = 'bi_quad_ant.xml'; + +[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory +[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder + +%% write openEMS compatible xml-file +WriteOpenEMS([Sim_Path '/' Sim_CSX], FDTD, CSX); + +%% show the structure +CSXGeomPlot([Sim_Path '/' Sim_CSX]); + +%% run openEMS +RunOpenEMS(Sim_Path, Sim_CSX); + +%% postprocessing & do the plots +port = calcPort(port, Sim_Path, freq); +s11 = port.uf.ref ./ port.uf.inc; + +% plot reflection coefficient S11 +figure +plot( freq/1e9, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 ); +ylim([-30 0]); +grid on +title( 'reflection coefficient S_{11}' ); +xlabel( 'frequency f / GHz' ); +ylabel( 'reflection coefficient |S_{11}|' ); + +%% calculate 3D far field pattern +phiRange = -180:2.5:180; +thetaRange = 0:2.5:180; + +nf2ff = CalcNF2FF(nf2ff, Sim_Path, f0, thetaRange*pi/180, phiRange*pi/180); + +disp( ['directivity: Dmax = ' num2str(10*log10(nf2ff.Dmax)) ' dBi'] ); + +% plot far-field pattern with Matlab +figure +plotFF3D(nf2ff, 'logscale', -20) + +%% +disp( 'Dumping far-field pattern to vtk (use Paraview to visualize)...' ); +DumpFF2VTK('Bi_Quad_Pattern.vtk', nf2ff.E_norm{1} / max(nf2ff.E_norm{1}(:)) * nf2ff.Dmax, thetaRange, phiRange, 'scale', 0.05); diff --git a/openEMS/matlab/examples/antennas/Patch_Antenna.m b/openEMS/matlab/examples/antennas/Patch_Antenna.m new file mode 100644 index 0000000..2011d6f --- /dev/null +++ b/openEMS/matlab/examples/antennas/Patch_Antenna.m @@ -0,0 +1,218 @@ +% +% EXAMPLE / antennas / patch antenna +% +% This example demonstrates how to: +% - calculate the reflection coefficient of a patch antenna +% +% +% Tested with +% - Matlab 2009b +% - Octave 3.3.52 +% - openEMS v0.0.23 +% +% (C) 2010,2011 Thorsten Liebig <thorsten.liebig@uni-due.de> + +close all +clear +clc + +%% switches & options... +postprocessing_only = 0; +draw_3d_pattern = 0; % this may take a while... +use_pml = 0; % use pml boundaries instead of mur +openEMS_opts = ''; + +%% setup the simulation +physical_constants; +unit = 1e-3; % all length in mm + +% width in x-direction +% length in y-direction +% main radiation in z-direction +patch.width = 32.86; % resonant length +patch.length = 41.37; + +substrate.epsR = 3.38; +substrate.kappa = 1e-3 * 2*pi*2.45e9 * EPS0*substrate.epsR; +substrate.width = 60; +substrate.length = 60; +substrate.thickness = 1.524; +substrate.cells = 4; + +feed.pos = -5.5; +feed.width = 2; +feed.R = 50; % feed resistance + +% size of the simulation box +SimBox = [100 100 25]; + +%% prepare simulation folder +Sim_Path = 'tmp'; +Sim_CSX = 'patch_ant.xml'; +if (postprocessing_only==0) + [status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory + [status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder +end + +%% setup FDTD parameter & excitation function +max_timesteps = 30000; +min_decrement = 1e-5; % equivalent to -50 dB +f0 = 0e9; % center frequency +fc = 3e9; % 20 dB corner frequency (in this case 0 Hz - 3e9 Hz) +FDTD = InitFDTD( 'NrTS', max_timesteps, 'EndCriteria', min_decrement ); +FDTD = SetGaussExcite( FDTD, f0, fc ); +BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % boundary conditions +if (use_pml>0) + BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'}; % use pml instead of mur +end +FDTD = SetBoundaryCond( FDTD, BC ); + +%% setup CSXCAD geometry & mesh +% currently, openEMS cannot automatically generate a mesh +max_res = c0 / (f0+fc) / unit / 20; % cell size: lambda/20 +CSX = InitCSX(); +mesh.x = [-SimBox(1)/2 SimBox(1)/2 -substrate.width/2 substrate.width/2 feed.pos]; +% add patch mesh with 2/3 - 1/3 rule +mesh.x = [mesh.x -patch.width/2-max_res/2*0.66 -patch.width/2+max_res/2*0.33 patch.width/2+max_res/2*0.66 patch.width/2-max_res/2*0.33]; +mesh.x = SmoothMeshLines( mesh.x, max_res, 1.4); % create a smooth mesh between specified mesh lines +mesh.y = [-SimBox(2)/2 SimBox(2)/2 -substrate.length/2 substrate.length/2 -feed.width/2 feed.width/2]; +% add patch mesh with 2/3 - 1/3 rule +mesh.y = [mesh.y -patch.length/2-max_res/2*0.66 -patch.length/2+max_res/2*0.33 patch.length/2+max_res/2*0.66 patch.length/2-max_res/2*0.33]; +mesh.y = SmoothMeshLines( mesh.y, max_res, 1.4 ); +mesh.z = [-SimBox(3)/2 linspace(0,substrate.thickness,substrate.cells) SimBox(3) ]; +mesh.z = SmoothMeshLines( mesh.z, max_res, 1.4 ); +mesh = AddPML( mesh, [8 8 8 8 8 8] ); % add equidistant cells (air around the structure) +CSX = DefineRectGrid( CSX, unit, mesh ); + +%% create patch +CSX = AddMetal( CSX, 'patch' ); % create a perfect electric conductor (PEC) +start = [-patch.width/2 -patch.length/2 substrate.thickness]; +stop = [ patch.width/2 patch.length/2 substrate.thickness]; +CSX = AddBox(CSX,'patch',10,start,stop); + +%% create substrate +CSX = AddMaterial( CSX, 'substrate' ); +CSX = SetMaterialProperty( CSX, 'substrate', 'Epsilon', substrate.epsR, 'Kappa', substrate.kappa ); +start = [-substrate.width/2 -substrate.length/2 0]; +stop = [ substrate.width/2 substrate.length/2 substrate.thickness]; +CSX = AddBox( CSX, 'substrate', 0, start, stop ); + +%% create ground (same size as substrate) +CSX = AddMetal( CSX, 'gnd' ); % create a perfect electric conductor (PEC) +start(3)=0; +stop(3) =0; +CSX = AddBox(CSX,'gnd',10,start,stop); + +%% apply the excitation & resist as a current source +start = [feed.pos-.1 -feed.width/2 0]; +stop = [feed.pos+.1 +feed.width/2 substrate.thickness]; +[CSX] = AddLumpedPort(CSX, 5 ,1 ,feed.R, start, stop, [0 0 1], true); + +%% dump magnetic field over the patch antenna +CSX = AddDump( CSX, 'Ht_', 'DumpType', 1, 'DumpMode', 2); % cell interpolated +start = [-patch.width -patch.length substrate.thickness+1]; +stop = [ patch.width patch.length substrate.thickness+1]; +CSX = AddBox( CSX, 'Ht_', 0, start, stop ); + +%%nf2ff calc +[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', -SimBox/2, SimBox/2); + +if (postprocessing_only==0) + %% write openEMS compatible xml-file + WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); + + %% show the structure + CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); + + %% run openEMS + RunOpenEMS( Sim_Path, Sim_CSX, openEMS_opts ); +end + +%% postprocessing & do the plots +freq = linspace( max([1e9,f0-fc]), f0+fc, 501 ); +U = ReadUI( {'port_ut1','et'}, 'tmp/', freq ); % time domain/freq domain voltage +I = ReadUI( 'port_it1', 'tmp/', freq ); % time domain/freq domain current (half time step is corrected) + +% plot time domain voltage +figure +[ax,h1,h2] = plotyy( U.TD{1}.t/1e-9, U.TD{1}.val, U.TD{2}.t/1e-9, U.TD{2}.val ); +set( h1, 'Linewidth', 2 ); +set( h1, 'Color', [1 0 0] ); +set( h2, 'Linewidth', 2 ); +set( h2, 'Color', [0 0 0] ); +grid on +title( 'time domain voltage' ); +xlabel( 'time t / ns' ); +ylabel( ax(1), 'voltage ut1 / V' ); +ylabel( ax(2), 'voltage et / V' ); +% now make the y-axis symmetric to y=0 (align zeros of y1 and y2) +y1 = ylim(ax(1)); +y2 = ylim(ax(2)); +ylim( ax(1), [-max(abs(y1)) max(abs(y1))] ); +ylim( ax(2), [-max(abs(y2)) max(abs(y2))] ); + +% plot feed point impedance +figure +Zin = U.FD{1}.val ./ I.FD{1}.val; +plot( freq/1e6, real(Zin), 'k-', 'Linewidth', 2 ); +hold on +grid on +plot( freq/1e6, imag(Zin), 'r--', 'Linewidth', 2 ); +title( 'feed point impedance' ); +xlabel( 'frequency f / MHz' ); +ylabel( 'impedance Z_{in} / Ohm' ); +legend( 'real', 'imag' ); + +% plot reflection coefficient S11 +figure +uf_inc = 0.5*(U.FD{1}.val + I.FD{1}.val * 50); +if_inc = 0.5*(I.FD{1}.val - U.FD{1}.val / 50); +uf_ref = U.FD{1}.val - uf_inc; +if_ref = I.FD{1}.val - if_inc; +s11 = uf_ref ./ uf_inc; +plot( freq/1e6, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 ); +grid on +title( 'reflection coefficient S_{11}' ); +xlabel( 'frequency f / MHz' ); +ylabel( 'reflection coefficient |S_{11}|' ); + +P_in = 0.5*U.FD{1}.val .* conj( I.FD{1}.val ); + +%% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +f_res_ind = find(s11==min(s11)); +f_res = freq(f_res_ind); + +% calculate the far field at phi=0 degrees and at phi=90 degrees +thetaRange = (0:2:359) - 180; +phiRange = [0 90]; +disp( 'calculating far field at phi=[0 90] deg...' ); +nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, thetaRange*pi/180, phiRange*pi/180); + +Dlog=10*log10(nf2ff.Dmax); + +% display power and directivity +disp( ['radiated power: Prad = ' num2str(nf2ff.Prad) ' Watt']); +disp( ['directivity: Dmax = ' num2str(Dlog) ' dBi'] ); +disp( ['efficiency: nu_rad = ' num2str(100*nf2ff.Prad./real(P_in(f_res_ind))) ' %']); + +% display phi +figure +plotFFdB(nf2ff,'xaxis','theta','param',[1 2]); +drawnow + +if (draw_3d_pattern==0) + return +end + +%% calculate 3D pattern +phiRange = 0:2:360; +thetaRange = 0:2:180; +disp( 'calculating 3D far field...' ); +nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, thetaRange*pi/180, phiRange*pi/180, 'Verbose',2,'Outfile','nf2ff_3D.h5'); +figure +plotFF3D(nf2ff); + + +%% visualize magnetic fields +% you will find vtk dump files in the simulation folder (tmp/) +% use paraview to visulaize them diff --git a/openEMS/matlab/examples/antennas/Patch_Antenna_Array.m b/openEMS/matlab/examples/antennas/Patch_Antenna_Array.m new file mode 100644 index 0000000..40b3b46 --- /dev/null +++ b/openEMS/matlab/examples/antennas/Patch_Antenna_Array.m @@ -0,0 +1,256 @@ +% +% EXAMPLE / antennas / patch antenna array +% +% This example demonstrates how to: +% - calculate the reflection coefficient of a patch antenna array +% +% +% Tested with +% - Matlab 2009b +% - Octave 3.3.52 +% - openEMS v0.0.23 +% +% (C) 2010 Thorsten Liebig <thorsten.liebig@uni-due.de> + +close all +clear +clc + +%% switches & options... +postprocessing_only = 0; +draw_3d_pattern = 0; % this may take a (very long) while... +use_pml = 0; % use pml boundaries instead of mur +openEMS_opts = ''; + +%% setup the simulation +physical_constants; +unit = 1e-3; % all length in mm + +% width in x-direction +% length in y-direction +% main radiation in z-direction +patch.width = 32.86; % resonant length +patch.length = 41.37; + +% define array size and dimensions +array.xn = 4; +array.yn = 4; +array.x_spacing = patch.width * 3; +array.y_spacing = patch.length * 3; + +substrate.epsR = 3.38; +substrate.kappa = 1e-3 * 2*pi*2.45e9 * EPS0*substrate.epsR; +substrate.width = 60 + (array.xn-1) * array.x_spacing; +substrate.length = 60 + (array.yn-1) * array.y_spacing; +substrate.thickness = 1.524; +substrate.cells = 4; + +feed.pos = -5.5; +feed.width = 2; +feed.R = 50; % feed resistance + +% size of the simulation box around the array +SimBox = [50+substrate.width 50+substrate.length 25]; + +%% prepare simulation folder +Sim_Path = 'tmp'; +Sim_CSX = 'patch_array.xml'; +if (postprocessing_only==0) + [status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory + [status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder +end + +%% setup FDTD parameter & excitation function +max_timesteps = 30000; +min_decrement = 1e-5; % equivalent to -50 dB +f0 = 0e9; % center frequency +fc = 3e9; % 10 dB corner frequency (in this case 0 Hz - 3e9 Hz) +FDTD = InitFDTD( 'NrTS', max_timesteps, 'EndCriteria', min_decrement ); +FDTD = SetGaussExcite( FDTD, f0, fc ); +BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % boundary conditions +if (use_pml>0) + BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'}; % use pml instead of mur +end +FDTD = SetBoundaryCond( FDTD, BC ); + +%% setup CSXCAD geometry & mesh +% currently, openEMS cannot automatically generate a mesh +max_res = c0 / (f0+fc) / unit / 20; % cell size: lambda/20 +CSX = InitCSX(); +mesh.x = [-SimBox(1)/2 SimBox(1)/2 -substrate.width/2 substrate.width/2]; +mesh.y = [-SimBox(2)/2 SimBox(2)/2 -substrate.length/2 substrate.length/2]; + +mesh.z = [-SimBox(3)/2 linspace(0,substrate.thickness,substrate.cells) SimBox(3) ]; +mesh.z = SmoothMeshLines( mesh.z, max_res, 1.4 ); + +for xn=1:array.xn + for yn=1:array.yn + midX = (array.xn/2 - xn + 1/2) * array.x_spacing; + midY = (array.yn/2 - yn + 1/2) * array.y_spacing; + + % feeding mesh + mesh.x = [mesh.x midX+feed.pos]; + mesh.y = [mesh.y midY-feed.width/2 midY+feed.width/2]; + + % add patch mesh with 2/3 - 1/3 rule + mesh.x = [mesh.x midX-patch.width/2-max_res/2*0.66 midX-patch.width/2+max_res/2*0.33 midX+patch.width/2+max_res/2*0.66 midX+patch.width/2-max_res/2*0.33]; + % add patch mesh with 2/3 - 1/3 rule + mesh.y = [mesh.y midY-patch.length/2-max_res/2*0.66 midY-patch.length/2+max_res/2*0.33 midY+patch.length/2+max_res/2*0.66 midY+patch.length/2-max_res/2*0.33]; + end +end +mesh.x = SmoothMeshLines( mesh.x, max_res, 1.4); % create a smooth mesh between specified mesh lines +mesh.y = SmoothMeshLines( mesh.y, max_res, 1.4 ); + +mesh = AddPML( mesh, [8 8 8 8 8 8] ); % add equidistant cells (air around the structure) +CSX = DefineRectGrid( CSX, unit, mesh ); + +%% create substrate +CSX = AddMaterial( CSX, 'substrate' ); +CSX = SetMaterialProperty( CSX, 'substrate', 'Epsilon', substrate.epsR, 'Kappa', substrate.kappa); +start = [-substrate.width/2 -substrate.length/2 0]; +stop = [ substrate.width/2 substrate.length/2 substrate.thickness]; +CSX = AddBox( CSX, 'substrate', 0, start, stop ); + +%% create ground (same size as substrate) +CSX = AddMetal( CSX, 'gnd' ); % create a perfect electric conductor (PEC) +start(3)=0; +stop(3) =0; +CSX = AddBox(CSX,'gnd',10,start,stop); +%% +CSX = AddMetal( CSX, 'patch' ); % create a perfect electric conductor (PEC) +number = 1; +for xn=1:array.xn + for yn=1:array.yn + + midX = (array.xn/2 - xn + 1/2) * array.x_spacing; + midY = (array.yn/2 - yn + 1/2) * array.y_spacing; + + % create patch + start = [midX-patch.width/2 midY-patch.length/2 substrate.thickness]; + stop = [midX+patch.width/2 midY+patch.length/2 substrate.thickness]; + CSX = AddBox(CSX,'patch',10,start,stop); + + % apply the excitation & resist as a current source + start = [midX+feed.pos-feed.width/2 midY-feed.width/2 0]; + stop = [midX+feed.pos+feed.width/2 midY+feed.width/2 substrate.thickness]; + [CSX] = AddLumpedPort(CSX, 5, number,feed.R, start, stop,[0 0 1],true); + number=number+1; + end +end + +%%nf2ff calc +[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', -SimBox/2, SimBox/2); + +if (postprocessing_only==0) + %% write openEMS compatible xml-file + WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); + + %% show the structure + CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); + + %% run openEMS + RunOpenEMS( Sim_Path, Sim_CSX, openEMS_opts ); +end + +%% postprocessing & do the plots +freq = linspace( max([1e9,f0-fc]), f0+fc, 501 ); +U = ReadUI( {'port_ut1','et'}, 'tmp/', freq ); % time domain/freq domain voltage +I = ReadUI( 'port_it1', 'tmp/', freq ); % time domain/freq domain current (half time step is corrected) + +% plot time domain voltage +figure +[ax,h1,h2] = plotyy( U.TD{1}.t/1e-9, U.TD{1}.val, U.TD{2}.t/1e-9, U.TD{2}.val ); +set( h1, 'Linewidth', 2 ); +set( h1, 'Color', [1 0 0] ); +set( h2, 'Linewidth', 2 ); +set( h2, 'Color', [0 0 0] ); +grid on +title( 'time domain voltage' ); +xlabel( 'time t / ns' ); +ylabel( ax(1), 'voltage ut1 / V' ); +ylabel( ax(2), 'voltage et / V' ); +% now make the y-axis symmetric to y=0 (align zeros of y1 and y2) +y1 = ylim(ax(1)); +y2 = ylim(ax(2)); +ylim( ax(1), [-max(abs(y1)) max(abs(y1))] ); +ylim( ax(2), [-max(abs(y2)) max(abs(y2))] ); + +% plot feed point impedance +figure +Zin = U.FD{1}.val ./ I.FD{1}.val; +plot( freq/1e6, real(Zin), 'k-', 'Linewidth', 2 ); +hold on +grid on +plot( freq/1e6, imag(Zin), 'r--', 'Linewidth', 2 ); +title( 'feed point impedance' ); +xlabel( 'frequency f / MHz' ); +ylabel( 'impedance Z_{in} / Ohm' ); +legend( 'real', 'imag' ); + +% plot reflection coefficient S11 +figure +uf_inc = 0.5*(U.FD{1}.val + I.FD{1}.val * 50); +if_inc = 0.5*(I.FD{1}.val - U.FD{1}.val / 50); +uf_ref = U.FD{1}.val - uf_inc; +if_ref = I.FD{1}.val - if_inc; +s11 = uf_ref ./ uf_inc; +plot( freq/1e6, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 ); +grid on +title( 'reflection coefficient S_{11}' ); +xlabel( 'frequency f / MHz' ); +ylabel( 'reflection coefficient |S_{11}|' ); + +%% +number = 1; +P_in = 0; +for xn=1:array.xn + for yn=1:array.yn + + U = ReadUI( ['port_ut' int2str(number)], 'tmp/', freq ); % time domain/freq domain voltage + I = ReadUI( ['port_it' int2str(number)], 'tmp/', freq ); % time domain/freq domain current (half time step is corrected) + + P_in = P_in + 0.5*U.FD{1}.val .* conj( I.FD{1}.val ); + number=number+1; + end +end + + +%% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +f_res_ind = find(s11==min(s11)); +f_res = freq(f_res_ind); + +% calculate the far field at phi=0 degrees and at phi=90 degrees +thetaRange = (0:2:359) - 180; +phiRange = [0 90]; +r = 1; % evaluate fields at radius r +disp( 'calculating far field at phi=[0 90] deg...' ); + +nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, thetaRange*pi/180, phiRange*pi/180); + +Dlog=10*log10(nf2ff.Dmax); + +% display power and directivity +disp( ['radiated power: Prad = ' num2str(nf2ff.Prad) ' Watt']); +disp( ['directivity: Dmax = ' num2str(Dlog) ' dBi'] ); +disp( ['efficiency: nu_rad = ' num2str(100*nf2ff.Prad./real(P_in(f_res_ind))) ' %']); + +% display phi +figure +plotFFdB(nf2ff,'xaxis','theta','param',[1 2]); +drawnow + +if (draw_3d_pattern==0) + return +end + +%% calculate 3D pattern +phiRange = 0:3:360; +thetaRange = unique([0:0.5:15 10:3:180]); +disp( 'calculating 3D far field...' ); +nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, thetaRange*pi/180, phiRange*pi/180, 'Verbose',2,'Outfile','nf2ff_3D.h5'); +figure +plotFF3D(nf2ff); + +%% visualize magnetic fields +% you will find vtk dump files in the simulation folder (tmp/) +% use paraview to visulaize them diff --git a/openEMS/matlab/examples/antennas/infDipol.m b/openEMS/matlab/examples/antennas/infDipol.m new file mode 100644 index 0000000..0a43bc8 --- /dev/null +++ b/openEMS/matlab/examples/antennas/infDipol.m @@ -0,0 +1,121 @@ +% +% infinitesimal dipole example +% + +close all +clear +clc + +postprocessing_only = 0; + +physical_constants + +% setup the simulation +drawingunit = 1e-6; % specify everything in um +Sim_Path = 'tmp'; +Sim_CSX = 'tmp.xml'; + +f_max = 1e9; +lambda = c0/f_max; + +% setup geometry values +dipole_length = lambda/50 /drawingunit; + + +dipole_orientation = 3; % 1,2,3: x,y,z + + +CSX = InitCSX(); + +% create an equidistant mesh +mesh.x = -dipole_length*10:dipole_length/2:dipole_length*10; +mesh.y = -dipole_length*10:dipole_length/2:dipole_length*10; +mesh.z = -dipole_length*10:dipole_length/2:dipole_length*10; + +% excitation +ex_vector = [0 0 0]; +ex_vector(dipole_orientation) = 1; +start = ex_vector * -dipole_length/2; +stop = ex_vector * dipole_length/2; +CSX = AddExcitation( CSX, 'infDipole', 1, ex_vector ); +% enlarge the box to be sure that one mesh line is covered by it +start = start - [0.1 0.1 0.1] * dipole_length/2; +stop = stop + [0.1 0.1 0.1] * dipole_length/2; +CSX = AddBox( CSX, 'infDipole', 1, start, stop ); + +% NFFF contour +start = [mesh.x(1) mesh.y(1) mesh.z(1) ]; +stop = [mesh.x(end) mesh.y(end) mesh.z(end) ]; +[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop); + +% add space for PML +mesh = AddPML( mesh, [8 8 8 8 8 8] ); +% define the mesh +CSX = DefineRectGrid( CSX, drawingunit, mesh ); + +if ~postprocessing_only + % setup FDTD parameters & excitation function + max_timesteps = 2000; + min_decrement = 1e-6; + FDTD = InitFDTD( 'NrTS', max_timesteps, 'EndCriteria', min_decrement, 'OverSampling',10 ); + FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2 ); + BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'}; + FDTD = SetBoundaryCond( FDTD, BC ); + + % Write openEMS compatible xml-file + [~,~,~] = rmdir(Sim_Path,'s'); + [~,~,~] = mkdir(Sim_Path); + WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + + % take a view at the "structure" + CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); + + % define openEMS options and start simulation + openEMS_opts = ''; + RunOpenEMS( Sim_Path, Sim_CSX, openEMS_opts ); +end + +%% post processing +disp( ' ' ); +disp( ' ********************************************************** ' ); +disp( ' ' ); + +% calculate the far field at phi=0 degrees and at phi=90 degrees +thetaRange = 0:0.5:359; +disp( 'calculating far field at phi=[0 90] deg..' ); +nf2ff = CalcNF2FF( nf2ff, Sim_Path, f_max, thetaRange/180*pi, [0 pi/2], 'Mode', 1 ); +Prad = nf2ff.Prad; +Dmax = nf2ff.Dmax; + +theta_HPBW = interp1(nf2ff.E_norm{1}(find(thetaRange<90),1)/max(nf2ff.E_norm{1}(find(thetaRange<90),1)),thetaRange(find(thetaRange<90)),1/sqrt(2))*2; + +% display power and directivity +disp( ['radiated power: Prad = ' num2str(Prad)] ); +disp( ['directivity: Dmax = ' num2str(Dmax)] ); +disp( ['theta_HPBW = ' num2str(theta_HPBW) ' °']); + +% display polar plot for the e-field magnitude for phi = 0 & 90 deg +figure +polarFF(nf2ff,'xaxis','theta','param',[1 2]); + +%% calculate the far field at theta=90 degrees +phiRange = 0:2:359; +disp( 'calculating far field at theta=90 deg..' ); +nf2ff = CalcNF2FF( nf2ff, Sim_Path, f_max, 90/180*pi, phiRange/180*pi, 'Mode', 1 ); + +% display polar plot +figure +polarFF(nf2ff,'xaxis','phi','param',1); + +%% calculate 3D pattern +phiRange = 0:5:360; +thetaRange = 0:5:180; +disp( 'calculating 3D far field...' ); +nf2ff = CalcNF2FF( nf2ff, Sim_Path, f_max, thetaRange/180*pi, phiRange/180*pi, 'Mode', 1 ); +figure +plotFF3D(nf2ff) + +%% +E_far_normalized = nf2ff.E_norm{1} / max(nf2ff.E_norm{1}(:)); +DumpFF2VTK([Sim_Path '/FF_pattern.vtk'],E_far_normalized, thetaRange, phiRange); +disp(['view the farfield pattern "' Sim_Path '/FF_pattern.vtk" using paraview' ]); diff --git a/openEMS/matlab/examples/antennas/inverted_f.m b/openEMS/matlab/examples/antennas/inverted_f.m new file mode 100644 index 0000000..175f94c --- /dev/null +++ b/openEMS/matlab/examples/antennas/inverted_f.m @@ -0,0 +1,205 @@ +% +% EXAMPLE / antennas / inverted-f antenna (ifa) 2.4GHz +% +% This example demonstrates how to: +% - calculate the reflection coefficient of an ifa +% - calculate farfield of an ifa +% +% Tested with +% - Octave 3.7.5 +% - openEMS v0.0.30+ (git 10.07.2013) +% +% (C) 2013 Stefan Mahr <dac922@gmx.de> + +close all +clear +clc + +%% setup the simulation +physical_constants; +unit = 1e-3; % all length in mm + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% substrate.width +% _______________________________________________ __ substrate. +% | A ifa.l |\ __ thickness +% | |ifa.e __________________________ | | +% | | | ___ _________________| w2 | | +% | | ifa.h | | || | | +% |_V_____________|___|___||______________________| | +% | .w1 .wf\ | | +% | |.fp| \ | | +% | | feed point | | +% | | | | substrate.length +% |<- substrate.width/2 ->| | | +% | | | +% |_______________________________________________| | +% \_______________________________________________\| +% +% Note: It's not checked whether your settings make sense, so check +% graphical output carefully. +% +substrate.width = 80; % width of substrate +substrate.length = 80; % length of substrate +substrate.thickness = 1.5; % thickness of substrate +substrate.cells = 4; % use 4 cells for meshing substrate + +ifa.h = 8; % height of short circuit stub +ifa.l = 22.5; % length of radiating element +ifa.w1 = 4; % width of short circuit stub +ifa.w2 = 2.5; % width of radiating element +ifa.wf = 1; % width of feed element +ifa.fp = 4; % position of feed element relative to short + % circuit stub +ifa.e = 10; % distance to edge + + +% substrate setup +substrate.epsR = 4.3; +substrate.kappa = 1e-3 * 2*pi*2.45e9 * EPS0*substrate.epsR; + +%setup feeding +feed.R = 50; %feed resistance + +%open AppCSXCAD and show ifa +show = 1; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% size of the simulation box +SimBox = [substrate.width*2 substrate.length*2 150]; + +%% setup FDTD parameter & excitation function +f0 = 2.5e9; % center frequency +fc = 1e9; % 20 dB corner frequency + +FDTD = InitFDTD('NrTS', 60000 ); +FDTD = SetGaussExcite( FDTD, f0, fc ); +BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % boundary conditions +FDTD = SetBoundaryCond( FDTD, BC ); + +%% setup CSXCAD geometry & mesh +CSX = InitCSX(); + +%initialize the mesh with the "air-box" dimensions +mesh.x = [-SimBox(1)/2 SimBox(1)/2]; +mesh.y = [-SimBox(2)/2 SimBox(2)/2]; +mesh.z = [-SimBox(3)/2 SimBox(3)/2]; + +%% create substrate +CSX = AddMaterial( CSX, 'substrate'); +CSX = SetMaterialProperty( CSX, 'substrate', 'Epsilon',substrate.epsR, 'Kappa', substrate.kappa); +start = [-substrate.width/2 -substrate.length/2 0]; +stop = [ substrate.width/2 substrate.length/2 substrate.thickness]; +CSX = AddBox( CSX, 'substrate', 1, start, stop ); +% add extra cells to discretize the substrate thickness +mesh.z = [linspace(0,substrate.thickness,substrate.cells+1) mesh.z]; + +%% create ground plane +CSX = AddMetal( CSX, 'groundplane' ); % create a perfect electric conductor (PEC) +start = [-substrate.width/2 -substrate.length/2 substrate.thickness]; +stop = [ substrate.width/2 substrate.length/2-ifa.e substrate.thickness]; +CSX = AddBox(CSX, 'groundplane', 10, start,stop); + +%% create ifa +CSX = AddMetal( CSX, 'ifa' ); % create a perfect electric conductor (PEC) +tl = [0,substrate.length/2-ifa.e,substrate.thickness]; % translate +start = [0 0.5 0] + tl; +stop = start + [ifa.wf ifa.h-0.5 0]; +CSX = AddBox( CSX, 'ifa', 10, start, stop); % feed element +start = [-ifa.fp 0 0] + tl; +stop = start + [-ifa.w1 ifa.h 0]; +CSX = AddBox( CSX, 'ifa', 10, start, stop); % short circuit stub +start = [(-ifa.fp-ifa.w1) ifa.h 0] + tl; +stop = start + [ifa.l -ifa.w2 0]; +CSX = AddBox( CSX, 'ifa', 10, start, stop); % radiating element + +ifa_mesh = DetectEdges(CSX, [], 'SetProperty','ifa'); +mesh.x = [mesh.x SmoothMeshLines(ifa_mesh.x, 0.5)]; +mesh.y = [mesh.y SmoothMeshLines(ifa_mesh.y, 0.5)]; + +%% apply the excitation & resist as a current source +start = [0 0 0] + tl; +stop = start + [ifa.wf 0.5 0]; +[CSX port] = AddLumpedPort(CSX, 5 ,1 ,feed.R, start, stop, [0 1 0], true); + +%% finalize the mesh +% generate a smooth mesh with max. cell size: lambda_min / 20 +mesh = DetectEdges(CSX, mesh); +mesh = SmoothMesh(mesh, c0 / (f0+fc) / unit / 20); +CSX = DefineRectGrid(CSX, unit, mesh); + +%% add a nf2ff calc box; size is 3 cells away from MUR boundary condition +start = [mesh.x(4) mesh.y(4) mesh.z(4)]; +stop = [mesh.x(end-3) mesh.y(end-3) mesh.z(end-3)]; +[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop); + +%% prepare simulation folder +Sim_Path = 'tmp_IFA'; +Sim_CSX = 'IFA.xml'; + +try confirm_recursive_rmdir(false,'local'); end + +[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory +[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder + +%% write openEMS compatible xml-file +WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); + +%% show the structure +if (show == 1) + CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); +end + + +%% run openEMS +RunOpenEMS( Sim_Path, Sim_CSX); %RunOpenEMS( Sim_Path, Sim_CSX, '--debug-PEC -v'); + +%% postprocessing & do the plots +freq = linspace( max([1e9,f0-fc]), f0+fc, 501 ); +port = calcPort(port, Sim_Path, freq); + +Zin = port.uf.tot ./ port.if.tot; +s11 = port.uf.ref ./ port.uf.inc; +P_in = real(0.5 * port.uf.tot .* conj( port.if.tot )); % antenna feed power + +% plot feed point impedance +figure +plot( freq/1e6, real(Zin), 'k-', 'Linewidth', 2 ); +hold on +grid on +plot( freq/1e6, imag(Zin), 'r--', 'Linewidth', 2 ); +title( 'feed point impedance' ); +xlabel( 'frequency f / MHz' ); +ylabel( 'impedance Z_{in} / Ohm' ); +legend( 'real', 'imag' ); + +% plot reflection coefficient S11 +figure +plot( freq/1e6, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 ); +grid on +title( 'reflection coefficient S_{11}' ); +xlabel( 'frequency f / MHz' ); +ylabel( 'reflection coefficient |S_{11}|' ); + +drawnow + +%% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%find resonance frequncy from s11 +f_res_ind = find(s11==min(s11)); +f_res = freq(f_res_ind); + +%% +disp( 'calculating 3D far field pattern and dumping to vtk (use Paraview to visualize)...' ); +thetaRange = (0:2:180); +phiRange = (0:2:360) - 180; +nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, thetaRange*pi/180, phiRange*pi/180,'Verbose',1,'Outfile','3D_Pattern.h5'); + +plotFF3D(nf2ff) + +% display power and directivity +disp( ['radiated power: Prad = ' num2str(nf2ff.Prad) ' Watt']); +disp( ['directivity: Dmax = ' num2str(nf2ff.Dmax) ' (' num2str(10*log10(nf2ff.Dmax)) ' dBi)'] ); +disp( ['efficiency: nu_rad = ' num2str(100*nf2ff.Prad./real(P_in(f_res_ind))) ' %']); + +E_far_normalized = nf2ff.E_norm{1} / max(nf2ff.E_norm{1}(:)) * nf2ff.Dmax; +DumpFF2VTK([Sim_Path '/3D_Pattern.vtk'],E_far_normalized,thetaRange,phiRange,1e-3); diff --git a/openEMS/matlab/examples/optimizer/optimizer_asco.m b/openEMS/matlab/examples/optimizer/optimizer_asco.m new file mode 100644 index 0000000..e405653 --- /dev/null +++ b/openEMS/matlab/examples/optimizer/optimizer_asco.m @@ -0,0 +1,36 @@ +% +% asco optimizer example -- optimize the turn number of a coil +% +% You need asco from http://asco.sf.net +% This is the main script. +% - optimizer_simfun.m starts the simulator with a parameter set from +% asco +% +% The goal is evaluated inside optimizer_simfun() to get as close to 2 uH. + +% clear +clear +close all +clc + +% setup the parameters +params = []; +params(end+1).name = 'turns'; +params(end).range = [1 30]; +params(end).value = 4; +params(end).step = 1; +params(end).active = 1; % this parameter is to be optimized + +% setup the simulation function +folder = fileparts( mfilename('fullpath') ); +options.simfun = [folder '/optimizer_simfun.m']; + +% additional options +% options.octave_exe = 'octave'; % must be newer than 3.2.4 (3.3.54 works) +options.clean = 1; + +% start the optimization +[params_opt,result] = optimize( 'opttmp', params, options, 'asco' ); + +% display best value +disp( ['ASCO found the optimum turn number: ' num2str(params_opt(1).value) ' result: ' num2str(result)] ); diff --git a/openEMS/matlab/examples/optimizer/optimizer_simfun.m b/openEMS/matlab/examples/optimizer/optimizer_simfun.m new file mode 100644 index 0000000..14152de --- /dev/null +++ b/openEMS/matlab/examples/optimizer/optimizer_simfun.m @@ -0,0 +1,134 @@ +function result = optimizer_simfun( folder, params ) +% +% simulation function +% +% the variable params contains the simulation parameters + +disp( [mfilename ': SIMULATING...'] ); + +if nargin == 0 + % visualize the structure if called without parameters + folder = 'tmp'; + params.turns = 10; +end + +oldpwd = pwd; +[a,a,a] = mkdir( folder ); +cd( folder ); + +% create the structure +f_max = 50e6; +structure( params, 'tmp.xml', f_max ); + +if nargin == 0 + % visualize the structure + CSXGeomPlot('tmp.xml'); + return; +end + +% start simulation +RunOpenEMS( '.', 'tmp.xml', '--engine=fastest' ); + +% postprocess the results +L = postproc( 'tmp.xml', f_max ); +disp( ['DONE. L = ' num2str(L(1)/1e-6) ' uH'] ); + +% calculate result +goal = 2e-6; % specify the goal: 2 uH +result = abs(goal - L(1)); % costs must not be negative + +% restore curent folder +cd( oldpwd ); + + + + +% ------------------------------------------------------------------------- +% ------------------------------------------------------------------------- + +function CSX = structure( params, filename, f_max ) +% CSX = structure( params, filename ) + +unit = 1e-3; % specify length in mm +lambda = 3e8/f_max; +resolution = lambda/15/unit; +mesh_size = 1.5; +radius = 10; +turns = params.turns; +height = 4 * turns; +feed_length = 10; + +CSX = InitCSX(); + +%% create coil +p1 = create_coil( radius, height, turns ); +p = p1(:,end) + [feed_length;0;0]; +p1 = [p1 p]; +p = p1(:,1) + [feed_length;0;0]; +p1 = [p p1]; +CSX = AddMetal(CSX,'PEC1'); +CSX = AddCurve(CSX, 'PEC1', 0, p1); + +%% create mesh +extraspace = 5*radius; +mesh.x = linspace(-radius,radius,ceil(2*radius/mesh_size)); +mesh.x = [mesh.x mesh.x(1)-extraspace mesh.x(end)+extraspace]; +mesh.x = [mesh.x p1(1,1) p1(1,1)-mesh_size p1(1,1)+mesh_size]; +mesh.x = SmoothMeshLines2( mesh.x, resolution ); +mesh.y = linspace(-radius,radius,ceil(2*radius/mesh_size)); +mesh.y = [mesh.y mesh.y(1)-extraspace mesh.y(end)+extraspace]; +mesh.y = SmoothMeshLines2( mesh.y, resolution ); +mesh.z = linspace(0,height,ceil(height/mesh_size)); +mesh.z = [mesh.z mesh.z(1)-extraspace mesh.z(end)+extraspace]; +% mesh.z = [mesh.z p1(3,1) p1(3,1)-mesh_size p1(3,1)+mesh_size]; +mesh.z = SmoothMeshLines2( mesh.z, resolution ); +CSX = DefineRectGrid(CSX, unit, mesh); + + +%% create port +[CSX,port] = AddCurvePort( CSX, 10, 1, 50, p1(:,1), p1(:,end), 'excite' ); + +if nargin > 1 + max_timesteps = 100000; + min_decrement = 1e-5; + FDTD = InitFDTD( max_timesteps, min_decrement ); + FDTD = SetGaussExcite( FDTD, 0, f_max ); + BC = {'PEC' 'PEC' 'PEC' 'PEC' 'PMC' 'PMC'}; + FDTD = SetBoundaryCond( FDTD, BC ); + + WriteOpenEMS( filename, FDTD, CSX ); +end + + +function p = create_coil(coil_rad,coil_length,coil_turns,coil_res,winding_direction,direction,offset,angle_offset) +if nargin < 8, angle_offset = 0; end +if nargin < 7, offset = [0; 0; 0]; end +if nargin < 6, direction = +1; end +if nargin < 5, winding_direction = +1; end +if nargin < 4, coil_res = 30; end +dt = 1/coil_res; +height = 0; + +p = []; +while abs(height) < coil_length + angle = height / (coil_length/coil_turns) * 2*pi; + p(1,end+1) = coil_rad * cos(angle*winding_direction+angle_offset); + p(2,end) = coil_rad * sin(angle*winding_direction+angle_offset); + p(3,end) = height * direction; + p(:,end) = p(:,end) + offset; + height = height + coil_length/coil_turns * dt; +end + + + +function L = postproc( filename, f_max ) +freq = linspace(0,f_max,201); +freq(1) = []; % delete DC component + +folder = fileparts(filename); +U = ReadUI( 'port_ut1', folder, freq ); +I = ReadUI( 'port_it1', folder, freq ); +Z = U.FD{1}.val ./ I.FD{1}.val; + +L = imag(Z) ./ (2*pi*freq); +L = reshape( L, 1, [] ); % row vector diff --git a/openEMS/matlab/examples/other/Helix.m b/openEMS/matlab/examples/other/Helix.m new file mode 100644 index 0000000..18e97c9 --- /dev/null +++ b/openEMS/matlab/examples/other/Helix.m @@ -0,0 +1,154 @@ +close all +clear +clc + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +feed_length=10; +wire_rad = sqrt(1.4/pi); +mesh_size = wire_rad; +coil_rad = 10; +coil_length = 50; +coil_turns = 8; +coil_res = 10; +port_length = mesh_size; %coil_length/2; +port_resist = 1000; + +f_max = 100e6; +f_excite = 300e6; + +%% define openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +openEMS_opts = ''; +% openEMS_opts = [openEMS_opts ' --debug-material']; +% openEMS_opts = [openEMS_opts ' --debug-boxes']; +% openEMS_opts = [openEMS_opts ' --debug-operator']; + +openEMS_opts = [openEMS_opts ' --disable-dumps --engine=fastest']; +% openEMS_opts = [openEMS_opts ' --engine=sse-compressed']; + +Sim_Path = 'tmp'; +Sim_CSX = 'helix.xml'; + +[status, message, messageid] = rmdir(Sim_Path,'s'); +[status,message,messageid] = mkdir(Sim_Path); + +%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitFDTD(30000,1e-6); +FDTD = SetGaussExcite(FDTD,f_excite/2,f_excite/2); +BC = [1 1 1 1 1 1]; +FDTD = SetBoundaryCond(FDTD,BC); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +add_Lines = mesh_size * 1.5.^(1:10); +add_Lines = add_Lines(find(add_Lines<(3e8/f_excite)/10*1e3)); + +CSX = InitCSX(); +mesh.x = -coil_rad-mesh_size : mesh_size : coil_rad+mesh_size+feed_length; +mesh.x = [mesh.x(1)-add_Lines mesh.x mesh.x(end)+add_Lines ]; +mesh.y = -coil_rad-mesh_size : mesh_size : coil_rad+mesh_size; +mesh.y = [mesh.y(1)-add_Lines mesh.y mesh.y(end)+add_Lines ]; +mesh.z = -mesh_size : mesh_size : coil_length+mesh_size; +mesh.z = [mesh.z(1)-add_Lines mesh.z mesh.z(end)+add_Lines ]; +CSX = DefineRectGrid(CSX, 1e-3,mesh); + +%% build/define helix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = AddMaterial(CSX,'copper'); +CSX = SetMaterialProperty(CSX,'copper','Kappa',56e6); + +dt = 1.0/coil_res; +height=0; +wire.Vertex = {}; +p(1,1) = coil_rad + feed_length; +p(2,1) = 0; +p(3,1) = 0.5*(coil_length-port_length); +p(1,2) = coil_rad + feed_length; +p(2,2) = 0; +p(3,2) = 0; +count=2; +for n=0:coil_turns-1 + for m=0:coil_res + count = count + 1; + p(1,count) = coil_rad * cos(2*pi*dt*m); + p(2,count) = coil_rad * sin(2*pi*dt*m); + p(3,count) = height + coil_length/coil_turns * dt*m; + end + height = height + coil_length/coil_turns; +end +p(1,count+1) = coil_rad + feed_length; +p(2,count+1) = 0; +p(3,count+1) = coil_length; +p(1,count+2) = coil_rad + feed_length; +p(2,count+2) = 0; +p(3,count+2) = 0.5*(coil_length+port_length); +CSX = AddWire(CSX, 'copper', 0, p, wire_rad); + +%% apply the excitation & resist as a current source%%%%%%%%%%%%%%%%%%%%%%% +CSX = AddMaterial(CSX,'resist'); +kappa = port_length/port_resist/wire_rad^2/pi/1e-3; +CSX = SetMaterialProperty(CSX,'resist','Kappa',kappa); + +start=[coil_rad+feed_length 0 (coil_length-port_length)/2]; +stop=[coil_rad+feed_length 0 (coil_length+port_length)/2]; +%start(3)=(coil_length-port_length)/2;stop(3)=(coil_length+port_length)/2; +CSX = AddCylinder(CSX,'resist',5 ,start,stop,wire_rad); + +CSX = AddExcitation(CSX,'excite',0,[0 0 1]); +CSX = AddCylinder(CSX,'excite', 0 ,start,stop,wire_rad); + +%% define voltage calc boxes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%voltage calc +CSX = AddProbe(CSX,'ut1',0); +CSX = AddBox(CSX,'ut1', 0 ,stop,start); + +%current calc +CSX = AddProbe(CSX,'it1',1); +start(3) = coil_length/2+mesh_size;stop(3) = coil_length/2+mesh_size; +start(1) = start(1)-2;start(2) = start(2)-2; +stop(1) = stop(1)+2;stop(2) = stop(2)+2; +CSX = AddBox(CSX,'it1', 0 ,start,stop); + +%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = AddDump(CSX,'Et_'); +start = [mesh.x(1) , 0 , mesh.z(1)]; +stop = [mesh.x(end) , 0 , mesh.z(end)]; +CSX = AddBox(CSX,'Et_',0 , start,stop); + +CSX = AddDump(CSX,'Ht_','DumpType',1); +start = [mesh.x(1) , 0 , mesh.z(1)]; +stop = [mesh.x(end) , 0 , mesh.z(end)]; +CSX = AddBox(CSX,'Ht_',0 , start,stop); + +%% Write openEMS compatoble xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + +%% run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts); + +%% postproc & do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +U = ReadUI('ut1','tmp/'); +I = ReadUI('it1','tmp/'); + +Z = U.FD{1}.val./I.FD{1}.val; +f = U.FD{1}.f; +L = imag(Z)./(f*2*pi); +R = real(Z); +ind = find(f<f_max); + +subplot(2,1,1); +plot(f(ind)*1e-6,L(ind)*1e9,'Linewidth',2); +xlabel('frequency (MHz)'); +ylabel('coil inductance (nH)'); +grid on; +subplot(2,1,2); +plot(f(ind)*1e-6,R(ind),'Linewidth',2); +hold on +plot(f(ind)*1e-6,imag(Z(ind)),'r','Linewidth',2); +xlabel('frequency (MHz)'); +ylabel('resistance (Ohm)'); +grid on; +legend( {'real','imaginary'}, 'location', 'northwest' ) + +figure +plot(U.TD{1}.t/1e-6,U.TD{1}.val,'Linewidth',2); +xlabel('time (us)'); +ylabel('amplitude (V)'); +grid on; diff --git a/openEMS/matlab/examples/other/LumpedElement.m b/openEMS/matlab/examples/other/LumpedElement.m new file mode 100644 index 0000000..d44094c --- /dev/null +++ b/openEMS/matlab/examples/other/LumpedElement.m @@ -0,0 +1,158 @@ +% +% EXAMPLE / other / lumped elements +% +% This example demonstrates how to: +% - use lumped elements +% +% +% Tested with +% - Matlab 2009b +% - openEMS v0.0.21-3 +% +% (C) 2010 Thorsten Liebig <thorsten.liebig@uni-due.de> + +close all +clear +clc + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +f_max = 100e6; +f_excite = 300e6; +SimBox = 100; +mesh_size = 2; + +Lumped.R = 1000; +Lumped.C = 10e-12; + +% the parasitice inductance of the feeding has to be deduced with a R=0 +% simulation +parasitic_L = 63e-9; + +%% define openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +openEMS_opts = ''; +% openEMS_opts = [openEMS_opts ' --debug-material']; +% openEMS_opts = [openEMS_opts ' --debug-boxes']; +% openEMS_opts = [openEMS_opts ' --debug-operator']; + +Sim_Path = 'tmp'; +Sim_CSX = 'lumped.xml'; + +[status, message, messageid] = rmdir(Sim_Path,'s'); +[status,message,messageid] = mkdir(Sim_Path); + +%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitFDTD(30000,1e-6); +FDTD = SetGaussExcite(FDTD,f_excite/2,f_excite/2); +BC = [1 1 1 1 1 1]; +FDTD = SetBoundaryCond(FDTD,BC); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX(); +mesh.x = SmoothMeshLines([-SimBox/2,+SimBox/2],mesh_size); +mesh.y = SmoothMeshLines([-SimBox/2,+SimBox/2],mesh_size); +mesh.z = SmoothMeshLines([-SimBox/2,+SimBox/2],mesh_size); +CSX = DefineRectGrid(CSX, 1e-3,mesh); + + +%% create structure +% insert curve port +start = [ 10 -10 0]; +stop = [ 10 10 0]; +CSX = AddCurvePort(CSX,0,1,100,start,stop,'excite'); + +% insert lumped element +CSX = AddLumpedElement( CSX, 'Capacitor', 1, 'C', Lumped.C, 'R', Lumped.R); +start = [ -14 -4 -4]; +stop = [ -6 4 4]; +CSX = AddBox( CSX, 'Capacitor', 0, start, stop ); + +% insert feeding wire +CSX = AddMetal(CSX,'metal'); +%first point +points(1,1) = -10; +points(2,1) = 4; +points(3,1) = 0; +%second point +points(1,2) = -10; +points(2,2) = 15; +points(3,2) = 0; +%3 point +points(1,end+1) = 10; +points(2,end) = 15; +points(3,end) = 0; +%4 point +points(1,end+1) = 10; +points(2,end) = 10; +points(3,end) = 0; +CSX = AddCurve(CSX,'metal', 10, points); + +points(2,:) = -1*points(2,:); +CSX = AddCurve(CSX,'metal', 10, points); + +%% Write openEMS compatoble xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + +% CSXGeomPlot([Sim_Path '/' Sim_CSX]); + +%% run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts); + +%% postproc & do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +f = linspace(1e6,f_max,1001); +w = 2*pi*f; +% read currents and voltages +U = ReadUI('port_ut1','tmp/',f); +I = ReadUI('port_it1','tmp/',f); + +% calculate analytic impedance +if (Lumped.R>=0) + Z_a = Lumped.R*(1-1i*w*Lumped.C*Lumped.R)./(1+(w*Lumped.C*Lumped.R).^2); +else + Z_a = -1i./(w*Lumped.C); +end + +% calculate numerical impedance +Z = U.FD{1}.val./I.FD{1}.val; + +% remove parasitic feeding effects +Z = Z - 1i*w*parasitic_L; + +L = imag(Z)./w; +C = -1./(w.*imag(Z)); +C(find(C<0)) = nan; +L(find(L<0)) = nan; +R = real(Z); + +subplot(2,1,1); +plot(f*1e-6,C*1e12,'Linewidth',2); +xlabel('frequency (MHz)'); +ylabel('capacitance (pF)'); +grid on; +subplot(2,1,2); +plot(f*1e-6,L*1e9,'Linewidth',2); +xlabel('frequency (MHz)'); +ylabel('inductance (nH)'); +grid on; + +figure(); +plot(f*1e-6,R,'Linewidth',2); +hold on +plot(f*1e-6,imag(Z),'r--','Linewidth',2); + +plot(f*1e-6,real(Z_a),'g-.','Linewidth',1); +plot(f*1e-6,imag(Z_a),'m--','Linewidth',1); + +xlabel('frequency (MHz)'); +ylabel('resistance (Ohm)'); +grid on; +legend( '\Re\{Z\}','\Im\{Z\}','\Re\{Z_{analytisch}\}','\Im\{Z_{analytisch}\}', 'location', 'northeast' ) + +figure(); +errorR = (R-real(Z_a))./R*100; +errorX = (imag(Z)-imag(Z_a))./imag(Z)*100; +plot(f*1e-6,errorR,'Linewidth',2); +hold on +grid on; +plot(f*1e-6,errorX,'r--','Linewidth',2); +xlabel('frequency (MHz)'); +ylabel('error (%)'); diff --git a/openEMS/matlab/examples/other/Metamaterial_PlaneWave_Drude.m b/openEMS/matlab/examples/other/Metamaterial_PlaneWave_Drude.m new file mode 100644 index 0000000..db72b3e --- /dev/null +++ b/openEMS/matlab/examples/other/Metamaterial_PlaneWave_Drude.m @@ -0,0 +1,147 @@ +%%%%%%%%%%%%%%%%%%%%%%% +% example demonstrating double drude meta-material +% +% tested with openEMS v0.0.28 +% +% author: Thorsten Liebig @ 2010,2012 +%%%%%%%%%%%%%%%%%%%%%%% + +close all +clear +clc + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +postproc_only = 0; %set to 1 if the simulation is already done + +Settings = []; +Settings.LogFile = 'openEMS.log'; + +pic_size = round([1400 1400/4]); %define the animation picture size + +%simulation domain setup (in mm) +length = 500; +width = 10; +mesh_res = 0.5; % mesh resolution +height = 3*mesh_res; % hight is ony 3 lines with PEC (top/bottom) --> quasi 2D + +%FDTD setup +f0 = 5e9; %center frequency +f_BW = f0/sqrt(2); %bandwidth +MTM.eps_R = 1; +MTM.mue_R = 1; +MTM.f0 = f0; %plasma frequency of the drude material +MTM.relaxTime = 5e-9; %relaxation time (smaller number results in greater losses, set to 0 to disable) +MTM.length = 250; %length of the metamaterial +N_TS = 5e4; %number of timesteps +endCriteria = 1e-5; %stop simulation if signal is at -50dB + +%constants +physical_constants + +%% define openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +openEMS_opts = '-vvv'; + +Sim_Path = 'MTM_PW_Drude'; +Sim_CSX = 'MTM_PW_Drude.xml'; + +if (postproc_only==0) + + if (exist(Sim_Path,'dir')) + rmdir(Sim_Path,'s'); + end + mkdir(Sim_Path); + + %% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% + FDTD = InitFDTD(N_TS,endCriteria,'OverSampling',10); + FDTD = SetGaussExcite(FDTD,0,2*f0); + BC = [1 1 0 0 2 2]; + FDTD = SetBoundaryCond(FDTD,BC); + + %% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + CSX = InitCSX(); + mesh.x = -width/2 : mesh_res : width/2; + mesh.y = -height/2 : mesh_res : height/2; + mesh.z = -length/2 : mesh_res : length/2; + CSX = DefineRectGrid(CSX, 1e-3,mesh); + + %% apply the plane wave excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + start=[-width/2 -height/2 ,mesh.z(3)]; + stop=[width/2 height/2 mesh.z(3)]; + CSX = AddExcitation(CSX,'excite',0,[0 1 0]); % excite E_y + CSX = AddBox(CSX,'excite',0 ,start,stop); + + %% apply drude material %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + CSX = AddLorentzMaterial(CSX,'drude'); + CSX = SetMaterialProperty(CSX,'drude','Epsilon',MTM.eps_R,'EpsilonPlasmaFrequency',MTM.f0,'EpsilonRelaxTime',MTM.relaxTime); + CSX = SetMaterialProperty(CSX,'drude','Mue',MTM.mue_R,'MuePlasmaFrequency',MTM.f0,'MueRelaxTime',MTM.relaxTime); + start=[mesh.x(1) mesh.y(1) -MTM.length/2]; + stop =[mesh.x(end) mesh.y(end) MTM.length/2]; + CSX = AddBox(CSX,'drude', 10 ,start,stop); + + %% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + CSX = AddDump(CSX,'Et','FileType',1,'SubSampling','10,10,1'); + start = [mesh.x(2) ,0 , mesh.z(1)]; + stop = [mesh.x(end-1) , 0 , mesh.z(end)]; + CSX = AddBox(CSX,'Et',0 , start,stop); + + %% Write openEMS compatoble xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + + %% run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts, Settings); + +end + +%% plot the drude type material dependency +f = linspace(0.1*f0,2*f0,501); +w = 2*pi*f; +epsr = MTM.eps_R * (1 - (2*pi*MTM.f0)^2./( w.^2 - 1j*w./MTM.relaxTime )); +muer = MTM.mue_R * (1 - (2*pi*MTM.f0)^2./( w.^2 - 1j*w./MTM.relaxTime )); +plot(f,real(epsr),'Linewidth',2); +hold on +grid on +plot(f,imag(epsr),'r--','Linewidth',2); +plot(f,real(muer),'c-.','Linewidth',2); +plot(f,imag(muer),'m-.','Linewidth',2); +ylim([-10 MTM.eps_R]) +% l=legend('\Re \epsilon_r','\Im \epsilon_r','\Re \mue_r','\Im \mue_r'); +l=legend('$\Re\{\varepsilon_r\}$','$\Im\{\varepsilon_r\}$','$\Re\{\mu_r\}$','$\Im\{\mu_r\}$'); +set(l,'Interpreter','latex','Fontsize',12) + +%% plot E-fields +freq = [f0/sqrt(2) f0 f0*sqrt(2)]; +field = ReadHDF5FieldData([Sim_Path '/Et.h5']); +mesh_h5 = ReadHDF5Mesh([Sim_Path '/Et.h5']); + +ET = ReadUI('et',Sim_Path); +ef = DFT_time2freq(ET.TD{1}.t,ET.TD{1}.val,freq); + +field_FD = GetField_TD2FD(field, freq); + +mesh.x = linspace(-500,500,numel(mesh_h5.lines{1})); %make animation wider... +mesh.y = mesh_h5.lines{2}; +mesh.z = mesh_h5.lines{3}; + +[X Z] = meshgrid(mesh.x,mesh.z); +X = X'; +Z = Z'; + +for n=1:numel(field_FD.FD.values) + Ec{n} = squeeze(field_FD.FD.values{n}/ef(n)); +end + +%% +figure('Position',[10 100 pic_size(1) pic_size(2)]); +phase = linspace(0,2*pi,21); +disp('press CTRL+C to stop animation'); +while (1) + for ph = phase(1:end-1) + for n=1:numel(Ec) + subplot(1,numel(Ec),n) + E = real(Ec{n}.*exp(1j*ph)); + surf(X,Z,E(:,:,2)); + title(['f_0 = ' num2str(freq(n)*1e-9) ' GHz']) + end + pause(0.1); + end +end diff --git a/openEMS/matlab/examples/other/PML_reflection_analysis.m b/openEMS/matlab/examples/other/PML_reflection_analysis.m new file mode 100644 index 0000000..f243a83 --- /dev/null +++ b/openEMS/matlab/examples/other/PML_reflection_analysis.m @@ -0,0 +1,196 @@ +% +% fake-PML parallel plate waveguide example +% +% this example analyzes the reflection coefficient of a vacuum-pml +% interface +% + +% +% currently this example uses a normal material with a certain conductivity +% profile and not a pml +% + +close all +% clear +clc + +physical_constants + + +postprocessing_only = 0; + + + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +drawingunit = 1e-6; % specify everything in um + +length = 10000; +epr = 1; + +mesh_res = [200 200 200]; +max_timesteps = 100000; +min_decrement = 1e-6; +f_max = 8e9; + +%% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitFDTD( max_timesteps, min_decrement ); +FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2 ); +BC = [0 0 1 1 0 0]; +FDTD = SetBoundaryCond( FDTD, BC ); + +%% mesh grading +N_pml = 8; +pml_delta = cumsum(mesh_res(1) * 1.0 .^ (1:N_pml)); +% pml_delta = cumsum([200 200 200 200 200]); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX(); +mesh.x = 0 : mesh_res(1) : length; +mesh.x = [mesh.x(1) - fliplr(pml_delta), mesh.x]; +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 ); + +%% fake pml %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +g = 2; % 2..3 +R0 = 1e-6; % requested analytical reflection coefficient +Zm = sqrt(MUE0/(EPS0*epr)); % calculate reflection for substrate/pml interface +delta = pml_delta(end) * drawingunit; +deltal = mean(diff(pml_delta)) * drawingunit; +kappa0 = -log(R0)*log(g)/( 2*Zm*deltal*(g^(delta/deltal)-1) ); + +% kappa0 = 1.05; +CSX = AddMaterial( CSX, 'pml_xmin' ); +CSX = SetMaterialProperty( CSX, 'pml_xmin', 'Epsilon', epr ); +CSX = SetMaterialProperty( CSX, 'pml_xmin', 'Kappa', kappa0 ); +CSX = SetMaterialProperty( CSX, 'pml_xmin', 'Sigma', kappa0 * MUE0/(EPS0*epr) ); +CSX = SetMaterialWeight( CSX, 'pml_xmin', 'Kappa', [num2str(g) '^((abs(x-100)-' num2str(abs(mesh.x(N_pml+1))) ')/(' num2str(deltal) '/' num2str(drawingunit) '))'] ); % g^(rho/deltal)*kappa0 +CSX = SetMaterialWeight( CSX, 'pml_xmin', 'Sigma', [num2str(g) '^((abs(x-100)-' num2str(abs(mesh.x(N_pml+1))) ')/(' num2str(deltal) '/' num2str(drawingunit) '))'] ); +start = [mesh.x(1), mesh.y(1), mesh.z(1)]; +stop = [100, mesh.y(end), mesh.z(end)]; +CSX = AddBox( CSX, 'pml_xmin', 1, start, stop ); + +figure +x = [-fliplr(pml_delta) 50]; +plot( x, kappa0 * g.^((abs(x-50)-abs(mesh.x(N_pml+1)))./(deltal/drawingunit)) ,'x-'); +xlabel( 'x / m' ); +ylabel( 'kappa' ); +figure +title( 'conductivity profile inside the material' ); + +%% 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 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 ); +idx = interp1( mesh.x, 1:numel(mesh.x), length*1/3, 'nearest' ); +start = [mesh.x(idx), mesh.y(3), mesh.z(1)]; +stop = [mesh.x(idx), 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 ); +idx = interp1( mesh.x, 1:numel(mesh.x), length*1/3, 'nearest' ); +start = [mesh.x(idx), mesh.y(1), mesh.z(3)]; +stop = [mesh.x(idx), 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 = 'PML_reflection_analysis.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, 1601 ); +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); % line impedance +Z = Ef ./ Hf; % impedance at measurement plane +gamma = (Z - Z0) ./ (Z + Z0); + +plot( f/1e9, 20*log10(abs(gamma)),'Linewidth',2); +xlabel('frequency (GHz)'); +ylabel('reflection coefficient gamma (dB)'); +grid on; +title( 'Reflection Coefficient' ); + +if exist('ref_1','var') + hold on + plot( f/1e9, ref_1,'--','Linewidth',2, 'Color', [1 0 0]); + hold off +end +ref_1 = 20*log10(abs(gamma)); diff --git a/openEMS/matlab/examples/other/PlaneWave.m b/openEMS/matlab/examples/other/PlaneWave.m new file mode 100644 index 0000000..5e0d3e8 --- /dev/null +++ b/openEMS/matlab/examples/other/PlaneWave.m @@ -0,0 +1,69 @@ +close all +clear +clc + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +length = 5000; +width = 300; +height = 200; +mesh_res = 15; +abs_length = mesh_res*10; + +EPS0 = 8.85418781762e-12; +MUE0 = 1.256637062e-6; + +%% define openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +openEMS_opts = ''; +% openEMS_opts = [openEMS_opts ' --disable-dumps']; +% openEMS_opts = [openEMS_opts ' --debug-material']; +openEMS_opts = [openEMS_opts ' --engine=fastest']; + +Sim_Path = 'tmp'; +Sim_CSX = 'plane_wave.xml'; + +if (exist(Sim_Path,'dir')) + rmdir(Sim_Path,'s'); +end +mkdir(Sim_Path); + +%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitFDTD(5000,1e-5,'OverSampling',10); +FDTD = SetGaussExcite(FDTD,0.5e9,0.5e9); +BC = [1 1 0 0 2 2]; +FDTD = SetBoundaryCond(FDTD,BC); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX(); +mesh.x = -width/2 : mesh_res : width/2; +mesh.y = -height/2 : mesh_res : height/2; +mesh.z = 0 : mesh_res : length; +CSX = DefineRectGrid(CSX, 1e-3,mesh); + +%% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +start=[-width/2 -height/2 mesh.z(3)]; +stop =[ width/2 height/2 mesh.z(3)]; +CSX = AddExcitation(CSX,'excite',0,[0 1 0]); +CSX = AddBox(CSX,'excite',0 ,start,stop); + +%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = AddDump(CSX,'Et','FileType',1,'SubSampling','4,4,4'); +start = [mesh.x(1) , mesh.y(1) , mesh.z(1)]; +stop = [mesh.x(end) , mesh.y(end) , mesh.z(end)]; +CSX = AddBox(CSX,'Et',0 , start,stop); + +CSX = AddDump(CSX,'Ht','DumpType',1,'FileType',1,'SubSampling','4,4,4','DumpMode',2); +CSX = AddBox(CSX,'Ht',0,start,stop); + +%% Write openEMS compatoble xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + +%% run openEMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts); + +%% do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +PlotArgs.slice = {mesh.x(round(end/2)) mesh.y(round(end/2)) mesh.z(round(end/2))}; +PlotArgs.pauseTime=0.01; +PlotArgs.component=1; +PlotArgs.Limit = 'auto'; + +PlotHDF5FieldData('tmp/Ht.h5',PlotArgs) diff --git a/openEMS/matlab/examples/other/gauss_excitation_test.m b/openEMS/matlab/examples/other/gauss_excitation_test.m new file mode 100644 index 0000000..a7e166e --- /dev/null +++ b/openEMS/matlab/examples/other/gauss_excitation_test.m @@ -0,0 +1,72 @@ +% +% this script evaluates the same gaussian excitation function, as openEMS does +% + +clear +close all +clc + +f0 = 0e9; +fc = 10e9; +dT = 1e-12; % sample time-step + + +sigma = 1/sqrt(8/9)/pi/fc; +t0 = sqrt(18)/sqrt(8/9)/pi/fc; + +len = 2 * 9/(2*pi*fc) / dT; % gauss length + +for n=1:len + t_(n) = (n-1)*dT; + ex(n) = cos(2*pi*f0*((n-1)*dT - 9/(2*pi*fc))) .* exp(-((t_(n)-t0)/sigma)^2/2); +end + +plot(t_/1e-9,ex) +xlabel( 'time (ns)' ); +ylabel( 'amplitude' ); + + +disp( ['Amplitude at t=0: ' num2str(20*log10(abs(ex(1))/1)) ' dB'] ); + +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) ) + +% 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_fft), 'r' ) +hold on + +if (f0==0) + Fw = sigma*sqrt(2*pi)*exp(-0.5*(sigma*2*pi*f).^2); + plot( f/1e9, 2*abs(Fw), 'g--' ) + legend('dft','fft','analytic') +else + legend('dft','fft') +end + +xlim([0 max(f)/1e9]) + +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)' ); + + + diff --git a/openEMS/matlab/examples/other/resistance_sheet.m b/openEMS/matlab/examples/other/resistance_sheet.m new file mode 100644 index 0000000..b63c00e --- /dev/null +++ b/openEMS/matlab/examples/other/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'} ); diff --git a/openEMS/matlab/examples/transmission_lines/CPW_Line.m b/openEMS/matlab/examples/transmission_lines/CPW_Line.m new file mode 100644 index 0000000..6dac636 --- /dev/null +++ b/openEMS/matlab/examples/transmission_lines/CPW_Line.m @@ -0,0 +1,125 @@ +% +% Tutorials / CPW_Line +% +% Describtion at: +% +% Tested with +% - Octave 3.8.1 +% - openEMS v0.0.32 +% +% (C) 2014 Thorsten Liebig <thorsten.liebig@gmx.de> + +close all +clear +clc + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +physical_constants; +unit = 1e-6; % specify everything in um +CPW_length = 40000; +CPW_port_length = 10000; +CPW_width = 1000; +CPW_gap = 140; +substrate_thickness = 512; +substrate_width = 5000 +substrate_epr = 3.66; +f_max = 10e9; +air_spacing = 7000 + +% use a finite line CPW waveguide +if 1 + feed_R = 50; + pml_add_cells = [8 8 8 8 8 8]; + feed_shift_cells = 0; + x_spacing = air_spacing; +else % or use a waveguide with start/end in a pml + feed_R = inf; % CPW ends in a pml --> disable termination resitance + feed_shift_cells = 10; % CPW ends in an 8 cells thick pml --> shift feed 10 cells + pml_add_cells = [0 0 8 8 8 8]; % do not add air-space in x-direction + x_spacing = 0; % do not add air-space in x-direction +end + +%% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitFDTD('EndCriteria', 1e-4); +FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2 ); +BC = [2 2 2 2 2 2]; +FDTD = SetBoundaryCond( FDTD, BC ); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX(); +resolution = c0/(f_max*sqrt(substrate_epr))/unit /30; % resolution of lambda/50 +edge_res = 40; +mesh.x = SmoothMeshLines( [0 CPW_length/2 CPW_length/2+x_spacing], resolution, 1.5 ,0 ); +mesh.x = unique(sort([-mesh.x mesh.x])); +mesh.y = SmoothMeshLines( [CPW_width/2+[-edge_res/3 +edge_res/3*2] CPW_gap+CPW_width/2+[-edge_res/3*2 +edge_res/3]], edge_res , 1.5 ,0); +mesh.y = SmoothMeshLines( [0 mesh.y], edge_res*2, 1.3 ,0); +mesh.y = SmoothMeshLines( [0 mesh.y substrate_width/2 substrate_width/2+air_spacing], resolution, 1.3 ,0); +mesh.y = unique(sort([-mesh.y mesh.y])); +mesh.z = SmoothMeshLines( [-air_spacing linspace(0,substrate_thickness,5) substrate_thickness+air_spacing], resolution ); + +mesh = AddPML(mesh, pml_add_cells); +CSX = DefineRectGrid( CSX, unit, mesh ); + +%% substrate +CSX = AddMaterial( CSX, 'RO4350B' ); +CSX = SetMaterialProperty( CSX, 'RO4350B', 'Epsilon', substrate_epr ); +start = [-CPW_length/2, -substrate_width/2, 0]; +stop = [+CPW_length/2, +substrate_width/2, substrate_thickness]; +CSX = AddBox( CSX, 'RO4350B', 0, start, stop ); + +%% +CSX = AddMetal( CSX, 'CPW_PORT' ); + +%% CPW port, with the measurement plane at the end of each port +portstart = [ -CPW_length/2 , -CPW_width/2, substrate_thickness]; +portstop = [ -CPW_length/2+CPW_port_length, CPW_width/2, substrate_thickness]; +[CSX,port{1}] = AddCPWPort( CSX, 999, 1, 'CPW_PORT', portstart, portstop, CPW_gap, 'x', [0 1 0], 'ExcitePort', true, 'FeedShift', feed_shift_cells*resolution, 'MeasPlaneShift', CPW_port_length, 'Feed_R', feed_R); + +portstart = [ CPW_length/2 , -CPW_width/2, substrate_thickness]; +portstop = [ CPW_length/2-CPW_port_length, CPW_width/2, substrate_thickness]; +[CSX,port{2}] = AddCPWPort( CSX, 999, 2, 'CPW_PORT', portstart, portstop, CPW_gap, 'x', [0 1 0], 'MeasPlaneShift', CPW_port_length, 'Feed_R', feed_R); + +%% CPW +CSX = AddMetal( CSX, 'CPW'); +start = [ -CPW_length/2+CPW_port_length, -CPW_width/2, substrate_thickness]; +stop = [ +CPW_length/2-CPW_port_length, CPW_width/2, substrate_thickness]; +CSX = AddBox(CSX, 'CPW', 999, start, stop); + +%% CPW grounds +CSX = AddMetal( CSX, 'GND' ); +start = [-CPW_length/2, -CPW_width/2-CPW_gap, substrate_thickness]; +stop = [+CPW_length/2, -substrate_width/2 , substrate_thickness]; +CSX = AddBox(CSX, 'GND', 999, start, stop); + +start = [-CPW_length/2, +CPW_width/2+CPW_gap, substrate_thickness]; +stop = [+CPW_length/2, +substrate_width/2 , substrate_thickness]; +CSX = AddBox(CSX, 'GND', 999, start, stop); + +%% write/show/run the openEMS compatible xml-file +Sim_Path = 'tmp'; +Sim_CSX = 'CPW.xml'; + +[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory +[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder + +WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); +CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); +RunOpenEMS( Sim_Path, Sim_CSX ); + +%% post-processing +close all +f = linspace( 1e6, f_max, 1601 ); +port = calcPort( port, Sim_Path, f, 'RefImpedance', 50); + +s11 = port{1}.uf.ref./ port{1}.uf.inc; +s21 = port{2}.uf.ref./ port{1}.uf.inc; + +plot(f/1e9,20*log10(abs(s11)),'k-','LineWidth',2); +hold on; +grid on; +plot(f/1e9,20*log10(abs(s21)),'r--','LineWidth',2); +legend('S_{11}','S_{21}'); +ylabel('S-Parameter (dB)','FontSize',12); +xlabel('frequency (GHz) \rightarrow','FontSize',12); + + diff --git a/openEMS/matlab/examples/transmission_lines/Finite_Stripline.m b/openEMS/matlab/examples/transmission_lines/Finite_Stripline.m new file mode 100644 index 0000000..da2b4a3 --- /dev/null +++ b/openEMS/matlab/examples/transmission_lines/Finite_Stripline.m @@ -0,0 +1,91 @@ +% example demonstrating the use of a stripline terminated by a resistance +% (c) 2013 Thorsten Liebig + +close all +clear +clc + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +physical_constants; +unit = 1e-6; % specify everything in um +SL_length = 50000; +SL_width = 520; +SL_height = 500; +substrate_thickness = SL_height; +substrate_epr = 3.66; +f_max = 7e9; + +Air_Spacer = 20000; + +%% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitFDTD(); +FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2 ); +BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; +FDTD = SetBoundaryCond( FDTD, BC ); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX(); +resolution = c0/(f_max*sqrt(substrate_epr))/unit /50; % resolution of lambda/50 +mesh.x = SmoothMeshLines( [-SL_length/2 0 SL_length/2], resolution, 1.5 ,0 ); +mesh.y = SmoothMeshLines( [0 SL_width/2+[-resolution/3 +resolution/3*2]/4], resolution/4 , 1.5 ,0); +mesh.y = SmoothMeshLines( [-10*SL_width -mesh.y mesh.y 10*SL_width], resolution, 1.3 ,0); +mesh.z = linspace(0,substrate_thickness,5); +mesh.z = sort(unique([mesh.z -mesh.z])); + +%% substrate +CSX = AddMaterial( CSX, 'RO4350B' ); +CSX = SetMaterialProperty( CSX, 'RO4350B', 'Epsilon', substrate_epr ); +start = [mesh.x(1), mesh.y(1), mesh.z(1)]; +stop = [mesh.x(end), mesh.y(end), mesh.z(end)]; +CSX = AddBox( CSX, 'RO4350B', 0, start, stop ); + +%% add air spacer +mesh.x = [mesh.x mesh.x(1)-Air_Spacer mesh.x(end)+Air_Spacer]; +mesh.y = [mesh.y mesh.y(1)-Air_Spacer mesh.y(end)+Air_Spacer]; +mesh.z = [mesh.z mesh.z(1)-Air_Spacer mesh.z(end)+Air_Spacer]; +mesh = SmoothMesh(mesh, c0/f_max/unit/20); +CSX = DefineRectGrid( CSX, unit, mesh ); + +%% SL port +CSX = AddMetal( CSX, 'PEC' ); +portstart = [ -SL_length/2, -SL_width/2, 0]; +portstop = [ 0, SL_width/2, 0]; +[CSX,port{1}] = AddStripLinePort( CSX, 999, 1, 'PEC', portstart, portstop, SL_height, 'x', [0 0 -1], 'ExcitePort', true, 'Feed_R', 50, 'MeasPlaneShift', SL_length/3); + +portstart = [+SL_length/2, -SL_width/2, 0]; +portstop = [0 , SL_width/2, 0]; +[CSX,port{2}] = AddStripLinePort( CSX, 999, 2, 'PEC', portstart, portstop, SL_height, 'x', [0 0 -1], 'MeasPlaneShift', SL_length/3, 'Feed_R', 50); + +% bottom PEC plane +CSX = AddBox(CSX, 'PEC', 999, [-SL_length/2 -10*SL_width -SL_height],[+SL_length/2 +10*SL_width -SL_height]); +% top PEC plane +CSX = AddBox(CSX, 'PEC', 999, [-SL_length/2 -10*SL_width SL_height],[+SL_length/2 +10*SL_width SL_height]); + +%% write/show/run the openEMS compatible xml-file +Sim_Path = ['tmp_' mfilename]; +Sim_CSX = 'stripline.xml'; + +[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory +[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder + +WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); +CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); +RunOpenEMS( Sim_Path, Sim_CSX ); + +%% post-processing +close all +f = linspace( 1e6, f_max, 1601 ); +port = calcPort( port, Sim_Path, f, 'RefImpedance', 50); + +s11 = port{1}.uf.ref./ port{1}.uf.inc; +s21 = port{2}.uf.ref./ port{1}.uf.inc; + +plot(f/1e9,20*log10(abs(s11)),'k-','LineWidth',2); +hold on; +grid on; +plot(f/1e9,20*log10(abs(s21)),'r--','LineWidth',2); +legend('S_{11}','S_{21}'); +ylabel('S-Parameter (dB)','FontSize',12); +xlabel('frequency (GHz) \rightarrow','FontSize',12); +ylim([-50 2]); + diff --git a/openEMS/matlab/examples/transmission_lines/MSL.m b/openEMS/matlab/examples/transmission_lines/MSL.m new file mode 100644 index 0000000..b6ec0b3 --- /dev/null +++ b/openEMS/matlab/examples/transmission_lines/MSL.m @@ -0,0 +1,185 @@ +% +% EXAMPLE / microstrip / MSL +% +% Microstrip line on air "substrate" in z-direction. +% +% This example demonstrates: +% - simple microstrip geometry +% - characteristic impedance +% - material grading function +% - geometric priority concept +% +% +% Tested with +% - Matlab 2009b +% - Octave 3.3.52 +% - openEMS v0.0.14 +% +% (C) 2010 Thorsten Liebig <thorsten.liebig@uni-due.de> + +close all +clear +clc + +%% setup the simulation +physical_constants; +unit = 1e-3; % all length in mm + +% geometry +abs_length = 100; % absorber length +length = 600; +width = 400; +height = 200; +MSL_width = 50; +MSL_height = 10; + +%% prepare simulation folder +Sim_Path = 'tmp'; +Sim_CSX = 'msl.xml'; +[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory +[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder + +%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% +max_timesteps = 2000; +min_decrement = 1e-5; % equivalent to -50 dB +f0 = 2e9; % center frequency +fc = 1e9; % 10 dB corner frequency (in this case 1e9 Hz - 3e9 Hz) +FDTD = InitFDTD( max_timesteps, min_decrement ); +FDTD = SetGaussExcite( FDTD, f0, fc ); +BC = {'PMC' 'PMC' 'PEC' 'PMC' 'PEC' 'PEC'}; +FDTD = SetBoundaryCond( FDTD, BC ); + +%% setup CSXCAD geometry & mesh +% very simple mesh +CSX = InitCSX(); +resolution = c0/(f0+fc)/unit /15; % resolution of lambda/15 +mesh.x = SmoothMeshLines( [-width/2, width/2, -MSL_width/2, MSL_width/2], resolution ); % create smooth lines from fixed lines +mesh.y = SmoothMeshLines( [linspace(0,MSL_height,5) MSL_height+1 height], resolution ); +mesh.z = SmoothMeshLines( [0 length], resolution ); +CSX = DefineRectGrid( CSX, unit, mesh ); + +%% create MSL +% attention! the skin effect is not simulated, because the MSL is +% discretized with only one cell! +CSX = AddMaterial( CSX, 'copper' ); +CSX = SetMaterialProperty( CSX, 'copper', 'Kappa', 56e6 ); +start = [-MSL_width/2, MSL_height, 0]; +stop = [ MSL_width/2, MSL_height+1, length]; +priority = 100; % the geometric priority is set to 100 +CSX = AddBox( CSX, 'copper', priority, start, stop ); + +%% add excitation below the strip +start = [-MSL_width/2, 0 , mesh.z(1)]; +stop = [ MSL_width/2, MSL_height, mesh.z(1)]; +CSX = AddExcitation( CSX, 'excite', 0, [0 -1 0] ); +CSX = AddBox( CSX, 'excite', 0, start, stop ); + +%% fake pml +% this "pml" is a normal material with graded losses +% electric and magnetic losses are related to give low reflection +% for normally incident TEM waves +finalKappa = 1/abs_length^2; +finalSigma = finalKappa*MUE0/EPS0; +CSX = AddMaterial( CSX, 'fakepml' ); +CSX = SetMaterialProperty( CSX, 'fakepml', 'Kappa', finalKappa ); +CSX = SetMaterialProperty( CSX, 'fakepml', 'Sigma', finalSigma ); +CSX = SetMaterialWeight( CSX, 'fakepml', 'Kappa', ['pow(z-' num2str(length-abs_length) ',2)'] ); +CSX = SetMaterialWeight( CSX, 'fakepml', 'Sigma', ['pow(z-' num2str(length-abs_length) ',2)'] ); +start = [mesh.x(1) mesh.y(1) length-abs_length]; +stop = [mesh.x(end) mesh.y(end) length]; +% the geometric priority is set to 0, which is lower than the priority +% of the MSL, thus the MSL (copper) has precendence +priority = 0; +CSX = AddBox( CSX, 'fakepml', priority, start, stop ); + +%% define dump boxes +start = [mesh.x(1), MSL_height/2, mesh.z(1)]; +stop = [mesh.x(end), MSL_height/2, mesh.z(end)]; +CSX = AddDump( CSX, 'Et_', 'DumpMode', 2 ); % cell interpolated +CSX = AddBox( CSX, 'Et_', 0, start, stop ); +CSX = AddDump( CSX, 'Ht_', 'DumpType', 1, 'DumpMode', 2 ); % cell interpolated +CSX = AddBox( CSX, 'Ht_', 0, start, stop ); + +%% define voltage calc box +% voltage calc boxes will automatically snap to the next mesh-line +CSX = AddProbe( CSX, 'ut1', 0 ); +zidx = interp1( mesh.z, 1:numel(mesh.z), length/2, 'nearest' ); +start = [0 MSL_height mesh.z(zidx)]; +stop = [0 0 mesh.z(zidx)]; +CSX = AddBox( CSX, 'ut1', 0, start, stop ); +% add a second voltage probe to compensate space offset between voltage and +% current +CSX = AddProbe( CSX, 'ut2', 0 ); +start = [0 MSL_height mesh.z(zidx+1)]; +stop = [0 0 mesh.z(zidx+1)]; +CSX = AddBox( CSX, 'ut2', 0, start, stop ); + +%% define current calc box +% current calc boxes will automatically snap to the next dual mesh-line +CSX = AddProbe( CSX, 'it1', 1 ); +xidx1 = interp1( mesh.x, 1:numel(mesh.x), -MSL_width/2, 'nearest' ); +xidx2 = interp1( mesh.x, 1:numel(mesh.x), MSL_width/2, 'nearest' ); +xdelta = diff(mesh.x); +yidx1 = interp1( mesh.y, 1:numel(mesh.y), MSL_height, 'nearest' ); +yidx2 = interp1( mesh.y, 1:numel(mesh.y), MSL_height+1, 'nearest' ); +ydelta = diff(mesh.y); +zdelta = diff(mesh.z); +start = [mesh.x(xidx1)-xdelta(xidx1-1)/2, mesh.y(yidx1)-ydelta(yidx1-1)/2, mesh.z(zidx)+zdelta(zidx)/2]; +stop = [mesh.x(xidx2)+xdelta(xidx2)/2, mesh.y(yidx2)+ydelta(yidx2)/2, mesh.z(zidx)+zdelta(zidx)/2]; +CSX = AddBox( CSX, 'it1', 0, start, stop ); + +%% write openEMS compatible xml-file +WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); + +%% show the structure +CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); + +%% run openEMS +openEMS_opts = ''; +openEMS_opts = [openEMS_opts ' --engine=fastest']; +% openEMS_opts = [openEMS_opts ' --debug-material']; +% openEMS_opts = [openEMS_opts ' --debug-boxes']; +RunOpenEMS( Sim_Path, Sim_CSX, openEMS_opts ); + +%% postprocess +freq = linspace( f0-fc, f0+fc, 501 ); +U = ReadUI( {'ut1','ut2','et'}, 'tmp/', freq ); % time domain/freq domain voltage +I = ReadUI( 'it1', 'tmp/', freq ); % time domain/freq domain current (half time step offset is corrected) + +% plot time domain voltage +figure +[ax,h1,h2] = plotyy( U.TD{1}.t/1e-9, U.TD{1}.val, U.TD{3}.t/1e-9, U.TD{3}.val ); +set( h1, 'Linewidth', 2 ); +set( h1, 'Color', [1 0 0] ); +set( h2, 'Linewidth', 2 ); +set( h2, 'Color', [0 0 0] ); +grid on +title( 'time domain voltage' ); +xlabel( 'time t / ns' ); +ylabel( ax(1), 'voltage ut1 / V' ); +ylabel( ax(2), 'voltage et / V' ); +% now make the y-axis symmetric to y=0 (align zeros of y1 and y2) +y1 = ylim(ax(1)); +y2 = ylim(ax(2)); +ylim( ax(1), [-max(abs(y1)) max(abs(y1))] ); +ylim( ax(2), [-max(abs(y2)) max(abs(y2))] ); + +% calculate characteristic impedance +% arithmetic mean of ut1 and ut2 -> voltage in the middle of ut1 and ut2 +U = (U.FD{1}.val + U.FD{2}.val) / 2; +Z = U ./ I.FD{1}.val; + +% plot characteristic impedance +figure +plot( freq/1e6, real(Z), 'k-', 'Linewidth', 2 ); +hold on +grid on +plot( freq/1e6, imag(Z), 'r--', 'Linewidth', 2 ); +title( 'characteristic impedance of MSL' ); +xlabel( 'frequency f / MHz' ); +ylabel( 'characteristic impedance Z / Ohm' ); +legend( 'real', 'imag' ); + +%% visualize electric and magnetic fields +% you will find vtk dump files in the simulation folder (tmp/) +% use paraview to visualize them diff --git a/openEMS/matlab/examples/transmission_lines/MSL_Losses.m b/openEMS/matlab/examples/transmission_lines/MSL_Losses.m new file mode 100644 index 0000000..95cad7b --- /dev/null +++ b/openEMS/matlab/examples/transmission_lines/MSL_Losses.m @@ -0,0 +1,102 @@ +% +% examples / microstrip / MSL_Losses +% +% This example demonstrates how to model sheet conductor losses +% +% Tested with +% - Matlab 2013a / Octave 3.8.1+ +% - openEMS v0.0.32 +% +% (C) 2012-2014 Thorsten Liebig <thorsten.liebig@gmx.de> + +close all +clear +clc + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +physical_constants; +unit = 1e-6; % specify everything in um +MSL.length = 10000; +MSL.port_dist = 5000; +MSL.width = 225; +MSL.conductivity = 41e6; +MSL.thickness = 35e-6; + +substrate.thickness = 250; +substrate.epr = 9.8; + +f_start = 0e9; +f_stop = 25e9; + +lambda = c0/f_stop; + +%% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitFDTD('endCriteria',1e-4); +FDTD = SetGaussExcite(FDTD,0.5*(f_start+f_stop),0.5*(f_stop-f_start)); +BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PEC' 'PML_8'}; +FDTD = SetBoundaryCond( FDTD, BC ); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX(); +resolution = c0/(f_stop*sqrt(substrate.epr))/unit /20; +mesh.x = SmoothMeshLines( [-MSL.length*0.5-MSL.port_dist 0 MSL.length*0.5+MSL.port_dist], resolution, 1.3 ,0 ); +mesh.y = SmoothMeshLines2( [0 MSL.width/2], resolution/6 , 1.3); +mesh.y = SmoothMeshLines( [-0.5*lambda/unit -mesh.y mesh.y 0.5*lambda/unit], resolution, 1.4); +mesh.z = SmoothMeshLines( [linspace(0,substrate.thickness,10) 0.5*lambda/unit], resolution ); +CSX = DefineRectGrid( CSX, unit, mesh ); + +%% substrate +CSX = AddMaterial( CSX, 'RO4350B' ); +CSX = SetMaterialProperty( CSX, 'RO4350B', 'Epsilon', substrate.epr ); +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 ports and lossy line +CSX = AddConductingSheet( CSX, 'gold', MSL.conductivity, MSL.thickness ); +portstart = [ mesh.x(1), -MSL.width/2, substrate.thickness]; +portstop = [ mesh.x(1)+MSL.port_dist, MSL.width/2, 0]; +[CSX, port{1}] = AddMSLPort( CSX, 999, 1, 'gold', portstart, portstop, 0, [0 0 -1], 'ExcitePort', true, 'FeedShift', 10*resolution, 'MeasPlaneShift', MSL.port_dist); + +portstart = [mesh.x(end), -MSL.width/2, substrate.thickness]; +portstop = [mesh.x(end)-MSL.port_dist, MSL.width/2, 0]; +[CSX, port{2}] = AddMSLPort( CSX, 999, 2, 'gold', portstart, portstop, 0, [0 0 -1], 'MeasPlaneShift', MSL.port_dist ); + +start = [mesh.x(1)+MSL.port_dist, -MSL.width/2, substrate.thickness]; +stop = [mesh.x(end)-MSL.port_dist, MSL.width/2, substrate.thickness]; +CSX = AddBox(CSX,'gold',500,start,stop); + +%% write/show/run the openEMS compatible xml-file +Sim_Path = 'tmp'; +Sim_CSX = 'msl.xml'; + +[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory +[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder + +WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); +CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); +RunOpenEMS( Sim_Path, Sim_CSX ,''); + +%% post-processing +close all +f = linspace( f_start, f_stop, 1601 ); +port = calcPort(port, Sim_Path, f, 'RefImpedance', 50); + +s11 = port{1}.uf.ref./ port{1}.uf.inc; +s21 = port{2}.uf.ref./ port{1}.uf.inc; + +plot(f/1e9,-20*log10(abs(s21)),'r--','LineWidth',2); +grid on; +hold on; +ylabel('-|S_21| (dB)','Interpreter','None'); +xlabel('frequency (GHz)'); + +%% plot 35um thickness loss model curve +% values extracted from http://wcalc.sourceforge.net/cgi-bin/microstrip.cgi +model.f = [1 2 2.5 3 4 5 7.5 10 12.5 15 17.5 20 25 ]; % frequency in GHz +model.loss = [3.0 4.2 4.7 5.2 5.9 6.6 8.1 9.38 10.5 11.5 12.4 13.2 14.65]; % loss in db/m + +plot(model.f, model.loss * MSL.length * unit ,'k-','LineWidth',1); +legend('FDTD simulated attenuation','t=35um, loss model by E. Hammerstad & F. Bekkadal','Location','NorthWest'); + + diff --git a/openEMS/matlab/examples/transmission_lines/Stripline.m b/openEMS/matlab/examples/transmission_lines/Stripline.m new file mode 100644 index 0000000..71fd774 --- /dev/null +++ b/openEMS/matlab/examples/transmission_lines/Stripline.m @@ -0,0 +1,78 @@ +% example demonstrating the use of a stripline terminated by the pml +% (c) 2013 Thorsten Liebig + +close all +clear +clc + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +physical_constants; +unit = 1e-6; % specify everything in um +SL_length = 50000; +SL_width = 520; +SL_height = 500; +substrate_thickness = SL_height; +substrate_epr = 3.66; +f_max = 7e9; + +%% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitFDTD(); +FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2 ); +BC = {'PML_8' 'PML_8' 'PMC' 'PMC' 'PEC' 'PEC'}; +FDTD = SetBoundaryCond( FDTD, BC ); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX(); +resolution = c0/(f_max*sqrt(substrate_epr))/unit /50; % resolution of lambda/50 +mesh.x = SmoothMeshLines( [-SL_length/2 0 SL_length/2], resolution, 1.5 ,0 ); +mesh.y = SmoothMeshLines( [0 SL_width/2+[-resolution/3 +resolution/3*2]/4], resolution/4 , 1.5 ,0); +mesh.y = SmoothMeshLines( [-10*SL_width -mesh.y mesh.y 10*SL_width], resolution, 1.3 ,0); +mesh.z = linspace(0,substrate_thickness,5); +mesh.z = sort(unique([mesh.z -mesh.z])); +CSX = DefineRectGrid( CSX, unit, mesh ); + +%% substrate +CSX = AddMaterial( CSX, 'RO4350B' ); +CSX = SetMaterialProperty( CSX, 'RO4350B', 'Epsilon', substrate_epr ); +start = [mesh.x(1), mesh.y(1), mesh.z(1)]; +stop = [mesh.x(end), mesh.y(end), mesh.z(end)]; +CSX = AddBox( CSX, 'RO4350B', 0, start, stop ); + +%% SL port +CSX = AddMetal( CSX, 'PEC' ); +portstart = [ mesh.x(1), -SL_width/2, 0]; +portstop = [ 0, SL_width/2, 0]; +[CSX,port{1}] = AddStripLinePort( CSX, 999, 1, 'PEC', portstart, portstop, SL_height, 'x', [0 0 -1], 'ExcitePort', true, 'FeedShift', 10*resolution, 'MeasPlaneShift', SL_length/3); + +portstart = [mesh.x(end), -SL_width/2, 0]; +portstop = [0 , SL_width/2, 0]; +[CSX,port{2}] = AddStripLinePort( CSX, 999, 2, 'PEC', portstart, portstop, SL_height, 'x', [0 0 -1], 'MeasPlaneShift', SL_length/3 ); + +%% write/show/run the openEMS compatible xml-file +Sim_Path = ['tmp_' mfilename]; +Sim_CSX = 'stripline.xml'; + +[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory +[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder + +WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); +CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); +RunOpenEMS( Sim_Path, Sim_CSX ); + +%% post-processing +close all +f = linspace( 1e6, f_max, 1601 ); +port = calcPort( port, Sim_Path, f, 'RefImpedance', 50); + +s11 = port{1}.uf.ref./ port{1}.uf.inc; +s21 = port{2}.uf.ref./ port{1}.uf.inc; + +plot(f/1e9,20*log10(abs(s11)),'k-','LineWidth',2); +hold on; +grid on; +plot(f/1e9,20*log10(abs(s21)),'r--','LineWidth',2); +legend('S_{11}','S_{21}'); +ylabel('S-Parameter (dB)','FontSize',12); +xlabel('frequency (GHz) \rightarrow','FontSize',12); +ylim([-50 2]); + diff --git a/openEMS/matlab/examples/transmission_lines/directional_coupler.m b/openEMS/matlab/examples/transmission_lines/directional_coupler.m new file mode 100644 index 0000000..8ba86f2 --- /dev/null +++ b/openEMS/matlab/examples/transmission_lines/directional_coupler.m @@ -0,0 +1,261 @@ +function directional_coupler +% +% EXAMPLE / microstrip / directional_coupler +% +% Stacked directional coupler in microstrip technology. +% +% This example demonstrates: +% - simple microstrip geometry +% - S-parameter calculation using the ypar-method +% - display of coupler parameters +% - display of S11 (smith chart) +% +% +% Tested with +% - Matlab 2010b +% - Octave 3.2.4 +% - openEMS v0.0.17 +% +% (C) 2010 Sebastian Held <sebastian.held@gmx.de> + +clear +close all +clc + +% sim settings +showStructure = 1; +runSimulation = 1; + +for n=1:4 + if n > 1, showStructure = 0; end + ports{n} = sim( n, showStructure, runSimulation ); +end +postprocess( ports ); + + + + +function ports = sim( simnr, showStructure, runSimulation ) +physical_constants + +% setup the simulation +drawingunit = 1e-6; % specify everything in um +Sim_Path = ['tmp' int2str(simnr)]; +Sim_CSX = 'tmp.xml'; +f_max = 100e6; +lambda = c0/f_max; + +% specify the coupler +pcb1.w = 147000; +pcb1.h = 54500; +pcb1.t = 1524; +pcb1.epr = 3; +msl1.w = 135000; +msl1.h = 2800; +pcb2.w = 107000; +pcb2.h = 14000; +pcb2.t = 1524; +pcb2.epr = 3; +msl2.w = 95000; +msl2.h = 4000; + + +CSX = InitCSX(); + +% create the mesh +mesh.x = [-pcb1.w/2 pcb1.w/2 -pcb2.w/2 pcb2.w/2 -msl1.w/2 msl1.w/2 -msl2.w/2 msl2.w/2]; +mesh.x = [mesh.x linspace(-msl2.w/2,-msl2.w/2+msl2.h, 5) linspace(msl2.w/2,msl2.w/2-msl2.h, 5)]; +mesh.y = [-pcb1.h/2 pcb1.h/2 -pcb2.h/2 pcb2.h/2 -msl1.h/2 msl1.h/2 -msl2.h/2 msl2.h/2]; +mesh.z = [linspace(0,pcb1.t,5) linspace(pcb1.t,pcb1.t+pcb2.t,5)]; +mesh.z = [mesh.z mesh.z(end)+10*(mesh.z(end)-mesh.z(1))]; % add space above pcb +res = lambda/sqrt(max([pcb1.epr,pcb2.epr])) / 20 / drawingunit; +mesh.x = SmoothMeshLines2(mesh.x,res); +mesh.y = SmoothMeshLines2(mesh.y,res); +mesh.z = SmoothMeshLines2(mesh.z,res); +mesh = AddPML( mesh, [8 8 8 8 8 8] ); % add space for PML +CSX = DefineRectGrid( CSX, drawingunit, mesh ); + +%% create the structure + +% microstrip +CSX = AddMetal( CSX, 'PEC' ); +start = [-msl1.w/2, -msl1.h/2, pcb1.t]; +stop = [ msl1.w/2, msl1.h/2, pcb1.t]; +priority = 100; % the geometric priority is set to 100 +CSX = AddBox( CSX, 'PEC', priority, start, stop ); + +% ground plane +CSX = AddMetal( CSX, 'PEC_ground' ); +start = [-pcb1.w/2, -pcb1.h/2, 0]; +stop = [ pcb1.w/2, pcb1.h/2, 0]; +CSX = AddBox( CSX, 'PEC_ground', priority, start, stop ); + +% substrate 1 +start = [-pcb1.w/2, -pcb1.h/2, 0]; +stop = [ pcb1.w/2, pcb1.h/2, pcb1.t]; +priority = 10; +CSX = AddMaterial( CSX, 'substrate1' ); +CSX = SetMaterialProperty( CSX, 'substrate1', 'Epsilon', pcb1.epr ); +CSX = AddBox( CSX, 'substrate1', priority, start, stop ); + +% substrate 2 +start = [-pcb2.w/2, -pcb2.h/2, pcb1.t]; +stop = [ pcb2.w/2, pcb2.h/2, pcb1.t+pcb2.t]; +priority = 10; +CSX = AddMaterial( CSX, 'substrate2' ); +CSX = SetMaterialProperty( CSX, 'substrate2', 'Epsilon', pcb2.epr ); +CSX = AddBox( CSX, 'substrate2', priority, start, stop ); + +% stripline +start = [-msl2.w/2, -msl2.h/2, pcb1.t+pcb2.t]; +stop = [ msl2.w/2, msl2.h/2, pcb1.t+pcb2.t]; +priority = 100; +CSX = AddBox( CSX, 'PEC', priority, start, stop ); + +% connections +start = [-msl2.w/2, -msl2.h/2, pcb1.t+pcb2.t]; +stop = [-msl2.w/2+msl2.h, -pcb2.h/2, pcb1.t+pcb2.t]; +priority = 100; +CSX = AddBox( CSX, 'PEC', priority, start, stop ); +start = [ msl2.w/2, -msl2.h/2, pcb1.t+pcb2.t]; +stop = [ msl2.w/2-msl2.h, -pcb2.h/2, pcb1.t+pcb2.t]; +priority = 100; +CSX = AddBox( CSX, 'PEC', priority, start, stop ); + +%% ports +% this project needs 4 simulations +for n=1:4 + portexcite{n} = []; +end +portexcite{simnr} = 'excite'; + +% port 1: input port +start = [-msl1.w/2, 0, pcb1.t]; +stop = [-msl1.w/2, 0, 0]; +[CSX ports{1}] = AddCurvePort( CSX, 999, 1, 50, start, stop, portexcite{1} ); +% port 2: output port +start = [msl1.w/2, 0, pcb1.t]; +stop = [msl1.w/2, 0, 0]; +[CSX ports{2}] = AddCurvePort( CSX, 999, 2, 50, start, stop, portexcite{2} ); +% port 3: coupled port +start = [-msl2.w/2+msl2.h/2, -pcb2.h/2, pcb1.t+pcb2.t]; +stop = [-msl2.w/2+msl2.h/2, -pcb2.h/2, 0]; +[CSX ports{3}] = AddCurvePort( CSX, 999, 3, 50, start, stop, portexcite{3} ); +% port 4: isolated port +start = [msl2.w/2-msl2.h/2, -pcb2.h/2, pcb1.t+pcb2.t]; +stop = [msl2.w/2-msl2.h/2, -pcb2.h/2, 0]; +[CSX ports{4}] = AddCurvePort( CSX, 999, 4, 50, start, stop, portexcite{4} ); + +%% setup FDTD parameters & excitation function +max_timesteps = 50000; +min_decrement = 1e-6; +FDTD = InitFDTD( max_timesteps, min_decrement ); +FDTD = SetGaussExcite( FDTD, 0, f_max ); +BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'}; +BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % faster +FDTD = SetBoundaryCond( FDTD, BC ); + +%% Write openEMS compatible xml-file +if runSimulation + [dummy,dummy,dummy] = rmdir(Sim_Path,'s'); +end +[dummy,dummy,dummy] = mkdir(Sim_Path); +WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + +if showStructure + CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); +end + +%% run openEMS +openEMS_opts = ''; +openEMS_opts = [openEMS_opts ' --engine=fastest']; +% openEMS_opts = [openEMS_opts ' --debug-material']; +% openEMS_opts = [openEMS_opts ' --debug-boxes']; +if runSimulation + RunOpenEMS( Sim_Path, Sim_CSX, openEMS_opts ); +end + + + + +function postprocess( ports ) +f = linspace( 0, 100e6, 201 ); +Y = calc_ypar( f, ports{1}, 'tmp' ); +R = 50; +S = y2s(Y,R); + +% insertion loss +IL_dB = -20 * log10(abs(squeeze(S(2,1,:)))); + +% coupling factor +CF_dB = -20 * log10(abs(squeeze(S(3,1,:)))); + +% isolation +I_dB = -20 * log10(abs(squeeze(S(4,1,:)))); + +% directivity +D_dB = -20 * log10(abs(squeeze(S(4,1,:) ./ S(3,1,:)))); + +figure +plot( f, [IL_dB CF_dB I_dB D_dB] ); +legend( {'insertion loss','coupling factor','isolation','directivity'} ); +title( ['performance of the coupler for a termination resistance of R=' num2str(R)] ); +grid on + +smithchart +S11 = squeeze(S(1,1,:)); +plot( real(S11), imag(S11) ); +legend( 'S_{11}' ); +title( ['performance of the coupler for a termination resistance of R=' num2str(R)] ); +axis( [-1 1 -1 1] ); + + + +function smithchart +% 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 + color = [.6 .6 .6]; + h = plot( sin(0:0.01:2*pi), cos(0:0.01:2*pi), 'Color', color ); + hg = hggroup; + set( h,'Parent',hg ); + hold on + plot( hg, 0.25+0.75*sin(0:0.01:2*pi), 0.75*cos(0:0.01:2*pi), 'Color', color ); + plot( hg, 0.5+0.5*sin(0:0.01:2*pi), 0.5*cos(0:0.01:2*pi), 'Color', color ); + plot( hg, 0.75+0.25*sin(0:0.01:2*pi), 0.25*cos(0:0.01:2*pi), 'Color', color ); + plot( hg, [-1 1], [0 0], 'Color', color ); + axis equal + axis off +end + + +function s = y2s(y, ZL) +% S = y2s(Y, ZL) +% +% Admittance to Scattering transformation +% for square matrices at multiple frequencies +% +% ZL defaults to 50 Ohm + +if nargin < 2 + ZL = 50; +end + +if size(size(y),2) > 2 + nF = size(y,3); +else + nF = 1; +end + +I = diag(ones(1, size(y,2)))/ZL; + +for i=1:nF + %s(:,:,i) = inv(I+y(:,:,i)) * (I-y(:,:,i)); + s(:,:,i) = (I+y(:,:,i)) \ (I-y(:,:,i)); +end diff --git a/openEMS/matlab/examples/waveguide/Circ_Waveguide.m b/openEMS/matlab/examples/waveguide/Circ_Waveguide.m new file mode 100644 index 0000000..9ee860e --- /dev/null +++ b/openEMS/matlab/examples/waveguide/Circ_Waveguide.m @@ -0,0 +1,207 @@ +% +% EXAMPLE / waveguide / circular waveguide +% +% This example demonstrates how to: +% - setup a circular waveguide +% - use analytic functions for waveguide excitations and voltage/current +% calculations +% +% +% Tested with +% - Matlab 2009b +% - openEMS v0.0.17 +% +% (C) 2010 Thorsten Liebig <thorsten.liebig@uni-due.de> + +close all +clear +clc + +%% switches & options... +postprocessing_only = 0; +use_pml = 0; % use pml boundaries instead of mur +openEMS_opts = ''; +% openEMS_opts = [openEMS_opts ' --disable-dumps']; + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +numTS = 1e5; %number of timesteps +length = 1000; %length of the waveguide +unit = 1e-3; %drawing unit used +rad = 300; %radius of the circular waveguide +mesh_res = [10 10 15]; %desired mesh resolution + +%excitation +f0 = 350e6; %center frequency +f0_BW = 25e6; %bandwidth: 10dB cut-off frequency + +physical_constants + +%% TE11 mode definitions (Pozar 3rd edition) +p11 = 1.841; +kc = p11 / rad /unit; +k = 2*pi*f0/C0; +fc = C0*kc/2/pi; +beta = sqrt(k^2 - kc^2); +n_eff = (beta/k); + +kc = kc*unit; %functions must be defined in drawing units +func_Er = [ num2str(-1/kc^2) '/rho*cos(a)*j1(' num2str(kc) '*rho)']; +func_Ea = [ num2str(1/kc) '*sin(a)*0.5*(j0(' num2str(kc) '*rho)-jn(2,' num2str(kc) '*rho))']; +func_Ex = ['(' func_Er '*cos(a) - ' func_Ea '*sin(a) )*(rho<' num2str(rad) ')']; +func_Ey = ['(' func_Er '*sin(a) + ' func_Ea '*cos(a) )*(rho<' num2str(rad) ')']; + +func_Ha = [ num2str(-1/kc^2,'%14.13f') '/rho*cos(a)*j1(' num2str(kc,'%14.13f') '*rho)']; +func_Hr = [ '-1*' num2str(1/kc,'%14.13f') '*sin(a)*0.5*(j0(' num2str(kc,'%14.13f') '*rho)-jn(2,' num2str(kc,'%14.13f') '*rho))']; +func_Hx = ['(' func_Hr '*cos(a) - ' func_Ha '*sin(a) )*(rho<' num2str(rad) ')']; +func_Hy = ['(' func_Hr '*sin(a) + ' func_Ha '*cos(a) )*(rho<' num2str(rad) ')']; + +%% define files and path %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Sim_Path = 'tmp'; +Sim_CSX = 'Circ_WG.xml'; + +if (postprocessing_only==0) + [status, message, messageid] = rmdir(Sim_Path,'s'); + [status, message, messageid] = mkdir(Sim_Path); +end + +%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitFDTD(numTS,1e-6,'OverSampling',5); +FDTD = SetGaussExcite(FDTD,f0,f0_BW); +BC = {'PEC','PEC','PEC','PEC','PEC','MUR'}; +if (use_pml>0) + BC = {'PEC','PEC','PEC','PEC','PEC','PML_8'}; +end +FDTD = SetBoundaryCond(FDTD,BC,'MUR_PhaseVelocity',C0 / n_eff); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX(); +mesh.x = -mesh_res(1)/2-rad:mesh_res(1):rad+mesh_res(1)/2; +mesh.y = -mesh_res(2)/2-rad:mesh_res(2):rad+mesh_res(2)/2; +mesh.z = 0 : mesh_res(3) : length; +CSX = DefineRectGrid(CSX, 1e-3,mesh); + +start = [0,0,0]; +stop = [0,0,length]; + +%%% fill everything with copper, priority 0 +CSX = AddMetal(CSX,'copper'); +% CSX = SetMaterialProperty(CSX,'copper','Kappa',56e6); +CSX = AddBox(CSX,'copper',0,[mesh.x(1) mesh.y(1) mesh.z(1)],[mesh.x(end) mesh.y(end) mesh.z(end)]); + +%%% cut out an air cylinder as circular waveguide... priority 5 +CSX = AddMaterial(CSX,'air'); +CSX = SetMaterialProperty(CSX,'air','Epsilon',1); +CSX = AddCylinder(CSX,'air', 5 ,start,stop,rad); + +CSX = AddExcitation(CSX,'excite',0,[1 1 0]); +weight{1} = func_Ex; +weight{2} = func_Ey; +weight{3} = 0; +CSX = SetExcitationWeight(CSX, 'excite', weight ); +CSX = AddCylinder(CSX,'excite', 5 ,[0 0 -0.1],[0 0 0.1],rad); + +%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = AddDump(CSX,'Et_','SubSampling','2,2,2','FileType',0,'DumpMode',2); +start = [mesh.x(1) , 0 , mesh.z(1)]; +stop = [mesh.x(end), 0 , mesh.z(end)]; +CSX = AddBox(CSX,'Et_',0 , start,stop); + +CSX = AddDump(CSX,'Ht_','SubSampling','2,2,2','DumpType',1,'FileType',0,'DumpMode',2); +CSX = AddBox(CSX,'Ht_',0,start,stop); + +%% define voltage calc boxes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%voltage calc +start = [mesh.x(1) mesh.y(1) mesh.z(10)]; +stop = [mesh.x(end) mesh.y(end) mesh.z(10)]; +CSX = AddProbe(CSX, 'ut1', 10, 1, [], 'ModeFunction',{func_Ex,func_Ey,0}); +CSX = AddBox(CSX, 'ut1', 0 ,start,stop); +CSX = AddProbe(CSX,'it1', 11, 1, [], 'ModeFunction',{func_Hx,func_Hy,0}); +CSX = AddBox(CSX,'it1', 0 ,start,stop); + +start = [mesh.x(1) mesh.y(1) mesh.z(end-10)]; +stop = [mesh.x(end) mesh.y(end) mesh.z(end-10)]; +CSX = AddProbe(CSX, 'ut2', 10, 1, [], 'ModeFunction',{func_Ex,func_Ey,0}); +CSX = AddBox(CSX, 'ut2', 0 ,start,stop); +CSX = AddProbe(CSX,'it2', 11, 1, [], 'ModeFunction',{func_Hx,func_Hy,0}); +CSX = AddBox(CSX,'it2', 0 ,start,stop); + +port_dist = mesh.z(end-10) - mesh.z(10); + +%% Write openEMS +if (postprocessing_only==0) + WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + + RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts); +end + +%% do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +freq = linspace(f0-f0_BW,f0+f0_BW,201); +U = ReadUI({'ut1','ut2'},[Sim_Path '/'],freq); +I = ReadUI({'it1','it2'},[Sim_Path '/'],freq); +Exc = ReadUI('et',Sim_Path,freq); + +k = 2*pi*freq/C0; +kc = p11 / rad /unit; +beta = sqrt(k.^2 - kc^2); + +ZL_a = Z0*k./beta ; + +uf1 = U.FD{1}.val./Exc.FD{1}.val; +uf2 = U.FD{2}.val./Exc.FD{1}.val; +if1 = I.FD{1}.val./Exc.FD{1}.val; +if2 = I.FD{2}.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; +plot(freq,20*log10(abs(s11)),'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 +figure() +ZL = uf1./if1; +plot(freq,real(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,ZL_a,'g-.','Linewidth',2); +legend('\Re\{ZL\}','\Im\{ZL\}','ZL-analytic','Location','Best'); + +% beta compare +figure() +da = angle(uf1_inc)-angle(uf2_inc); +da = mod(da,2*pi); +beta_12 = (da)/port_dist/unit; +plot(freq,beta_12,'Linewidth',2); +xlim([freq(1) freq(end)]); +xlabel('frequency (Hz)'); +ylabel('\beta (m^{-1})'); +grid on; +hold on; +plot(freq,beta,'g--','Linewidth',2); +legend('\beta-FDTD','\beta-analytic','Location','Best'); + +%% visualize electric & magnetic fields +disp('you will find vtk dump files in the simulation folder (tmp/)') +disp('use paraview to visulaize them'); diff --git a/openEMS/matlab/examples/waveguide/Circ_Waveguide_CylinderCoords.m b/openEMS/matlab/examples/waveguide/Circ_Waveguide_CylinderCoords.m new file mode 100644 index 0000000..a767c0a --- /dev/null +++ b/openEMS/matlab/examples/waveguide/Circ_Waveguide_CylinderCoords.m @@ -0,0 +1,204 @@ +% +% EXAMPLE / waveguide / circular waveguide cylindrical coordinates +% +% This example demonstrates how to: +% - use cylindrical coordinates +% - setup a circular waveguide defined by the boundary conditions of the +% cylindrical coordinate system +% - use analytic functions for waveguide excitations and voltage/current +% calculations +% +% +% Tested with +% - Matlab 2009b +% - openEMS v0.0.17 +% +% (C) 2010 Thorsten Liebig <thorsten.liebig@uni-due.de> + +close all +clear +clc + +%% switches & options... +postprocessing_only = 0; +use_pml = 0; % use pml boundaries instead of mur +use_MultiGrid = 1; % disable multi-grid for this example +openEMS_opts = ''; +% openEMS_opts = [openEMS_opts ' --disable-dumps']; + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +numTS = 1e5; %number of timesteps +length = 1000; %length of the waveguide +unit = 1e-3; %drawing unit used +rad = 300; %radius of the circular waveguide +mesh_res = [10 nan 15]; %desired mesh resolution +N_alpha = 50; %mesh lines in azimuth direction + +MultiGrid_Level = [50]; % define multigrid radii (if enabled) + +%excitation +f0 = 350e6; %center frequency +f0_BW = 25e6; %bandwidth: 10dB cut-off frequency + +physical_constants + +%% TE11 mode definitions (Pozar 3rd edition) +p11 = 1.841; +kc = p11 / rad /unit; +k = 2*pi*f0/C0; +fc = C0*kc/2/pi; +beta = sqrt(k^2 - kc^2); +n_eff = (beta/k); + +kc = kc*unit; %functions must be defined in drawing units +func_Er = [ num2str(-1/kc^2,15) '/rho*cos(a)*j1(' num2str(kc,15) '*rho)']; +func_Ea = [ num2str(1/kc,15) '*sin(a)*0.5*(j0(' num2str(kc,15) '*rho)-jn(2,' num2str(kc,15) '*rho))']; +func_Ha = [ num2str(-1/kc^2,'%14.13f') '/rho*cos(a)*j1(' num2str(kc,'%14.13f') '*rho)']; +func_Hr = [ '-1*' num2str(1/kc,'%14.13f') '*sin(a)*0.5*(j0(' num2str(kc,'%14.13f') '*rho)-jn(2,' num2str(kc,'%14.13f') '*rho))']; + +%% define files and path %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Sim_Path = 'tmp'; +Sim_CSX = 'Circ_WG_CC.xml'; + +if (postprocessing_only==0) + [status, message, messageid] = rmdir(Sim_Path,'s'); + [status, message, messageid] = mkdir(Sim_Path); +end + +%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if (use_MultiGrid==0) + FDTD = InitCylindricalFDTD(numTS,1e-5,'OverSampling',10); +else + mg_str = num2str(MultiGrid_Level,'%d,'); %create comma-separated string + N_alpha = round(N_alpha * 2^numel(MultiGrid_Level)); + FDTD = InitCylindricalFDTD(numTS,1e-5,'OverSampling',10,'MultiGrid',mg_str(1:end-1)); +end +FDTD = SetGaussExcite(FDTD,f0,f0_BW); +BC = {'PEC','PEC','PEC','PEC','PEC','MUR'}; +if (use_pml>0) + BC = {'PEC','PEC','PEC','PEC','PEC','PML_8'}; +end +FDTD = SetBoundaryCond(FDTD,BC,'MUR_PhaseVelocity',C0 / n_eff); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX('CoordSystem',1); +mesh.x = 0:mesh_res(1):rad; +%define an odd number of lines in alpha-direction +mesh.y = linspace(-pi,pi,N_alpha+mod(N_alpha+1,2))+pi/2; +mesh.z = 0 : mesh_res(3) : length; +CSX = DefineRectGrid(CSX, unit,mesh); + +%% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = AddExcitation(CSX,'excite',0,[1 1 0]); +weight{1} = func_Er; +weight{2} = func_Ea; +weight{3} = 0; +CSX = SetExcitationWeight(CSX, 'excite', weight ); +start = [mesh.x(1) mesh.y(1) mesh.z(1)]; +stop = [mesh.x(end) mesh.y(end) mesh.z(1)]; +CSX = AddBox(CSX,'excite', 5 ,start,stop); + +%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = AddDump(CSX,'Et_','FileType',0,'DumpMode',2,'SubSampling','2,2,2'); +start = [mesh.x(1) , 0 , mesh.z(1)]; +stop = [mesh.x(end), 0 , mesh.z(end)]; +CSX = AddBox(CSX,'Et_',0 , start,stop); + +CSX = AddDump(CSX,'Ht','FileType',0,'DumpType',1,'DumpMode',2,'SubSampling','2,2,2'); +CSX = AddBox(CSX,'Ht',0 , start,stop); + +%% define voltage calc boxes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +start = [mesh.x(1) mesh.y(1) mesh.z(10)]; +stop = [mesh.x(end) mesh.y(end) mesh.z(10)]; +CSX = AddProbe(CSX, 'ut1', 10, 1, [], 'ModeFunction',{func_Er,func_Ea,0}); +CSX = AddBox(CSX, 'ut1', 0 ,start,stop); +CSX = AddProbe(CSX,'it1', 11, 1, [], 'ModeFunction',{func_Hr,func_Ha,0}); +CSX = AddBox(CSX,'it1', 0 ,start,stop); + +start = [mesh.x(1) mesh.y(1) mesh.z(end-10)]; +stop = [mesh.x(end) mesh.y(end) mesh.z(end-10)]; +CSX = AddProbe(CSX, 'ut2', 10, 1, [], 'ModeFunction',{func_Er,func_Ea,0}); +CSX = AddBox(CSX, 'ut2', 0 ,start,stop); +CSX = AddProbe(CSX,'it2', 11, 1, [], 'ModeFunction',{func_Hr,func_Ha,0}); +CSX = AddBox(CSX,'it2', 0 ,start,stop); + +port_dist = mesh.z(end-10) - mesh.z(10); + +%% Write openEMS +if (postprocessing_only==0) + WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + + RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts); +end + +%% do the plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +freq = linspace(f0-f0_BW,f0+f0_BW,201); +U = ReadUI({'ut1','ut2'},[Sim_Path '/'],freq); +I = ReadUI({'it1','it2'},[Sim_Path '/'],freq); +Exc = ReadUI('et',Sim_Path,freq); + +k = 2*pi*freq/C0; +kc = p11 / rad /unit; +beta = sqrt(k.^2 - kc^2); + +ZL_a = Z0*k./beta ; + +uf1 = U.FD{1}.val./Exc.FD{1}.val; +uf2 = U.FD{2}.val./Exc.FD{1}.val; +if1 = I.FD{1}.val./Exc.FD{1}.val; +if2 = I.FD{2}.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; +plot(freq,20*log10(abs(s11)),'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 +figure() +ZL = uf1./if1; +plot(freq,real(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,ZL_a,'g-.','Linewidth',2); +legend('\Re\{ZL\}','\Im\{ZL\}','ZL-analytic','Location','Best'); + +%% beta compare +figure() +da = angle(uf1_inc)-angle(uf2_inc); +da = mod(da,2*pi); +beta_12 = (da)/port_dist/unit; +plot(freq,beta_12,'Linewidth',2); +xlim([freq(1) freq(end)]); +xlabel('frequency (Hz)'); +ylabel('\beta (m^{-1})'); +grid on; +hold on; +plot(freq,beta,'g--','Linewidth',2); +legend('\beta-FDTD','\beta-analytic','Location','Best'); + +%% visualize electric & magnetic fields +disp('you will find vtk dump files in the simulation folder (tmp/)') +disp('use paraview to visulaize them'); diff --git a/openEMS/matlab/examples/waveguide/Coax.m b/openEMS/matlab/examples/waveguide/Coax.m new file mode 100644 index 0000000..3db4d6a --- /dev/null +++ b/openEMS/matlab/examples/waveguide/Coax.m @@ -0,0 +1,120 @@ +% +% EXAMPLE / waveguide / coaxial cable +% +% This example demonstrates how to: +% - setup a coaxial waveguide +% - use analytic functions for waveguide excitations and voltage/current +% calculations +% +% +% Tested with +% - Matlab 2009b +% - openEMS v0.0.17 +% +% (C) 2010 Thorsten Liebig <thorsten.liebig@uni-due.de> + +close all +clear +clc + +%% switches & options... +postprocessing_only = 0; +use_pml = 0; % use pml boundaries instead of mur +openEMS_opts = ''; +% openEMS_opts = [openEMS_opts ' --disable-dumps']; + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +numTS = 5000; %number of timesteps +length = 1000; %length of the waveguide +unit = 1e-3; %drawing unit used +coax_rad_i = 100; %inner radius +coax_rad_ai = 230; %inner radius of outer cladding +coax_rad_aa = 240; %outer radius of outer cladding +mesh_res = [5 5 5]; %desired mesh resolution + +physical_constants; + +%excitation +f0 = 0.5e9; +epsR = 1; + +%% create sim path %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Sim_Path = 'tmp'; +Sim_CSX = 'coax.xml'; + +if (postprocessing_only==0) + [status, message, messageid] = rmdir(Sim_Path,'s'); + [status, message, messageid] = mkdir(Sim_Path); +end + +%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitFDTD(numTS,1e-5); +FDTD = SetGaussExcite(FDTD,f0,f0); +BC = {'PEC','PEC','PEC','PEC','MUR','MUR'}; +if (use_pml>0) + BC = {'PEC','PEC','PEC','PEC','PML_8','PML_8'}; +end +FDTD = SetBoundaryCond(FDTD,BC); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX(); +mesh.x = -coax_rad_aa : mesh_res(1) : coax_rad_aa; +mesh.y = mesh.x; +mesh.z = SmoothMeshLines([0 length], mesh_res(3)); +CSX = DefineRectGrid(CSX, unit, mesh); + +%%% coax +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) ); + +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); +start = [mesh.x(1) , 0 , mesh.z(1)]; +stop = [mesh.x(end) , 0 , mesh.z(end)]; +CSX = AddBox(CSX,'Et_',0 , start,stop); + +CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',2); +CSX = AddBox(CSX,'Ht_',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); +port = calcPort(port, Sim_Path, freq); + +%% plot s-parameter +figure +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)'); + +%% plot line-impedance comparison +figure() +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(port{1}.ZL),'r--','Linewidth',2); +plot(freq,ZL_a,'g-.','Linewidth',2); +legend('\Re\{ZL\}','\Im\{ZL\}','ZL-analytic','Location','Best'); diff --git a/openEMS/matlab/examples/waveguide/Coax_CylinderCoords.m b/openEMS/matlab/examples/waveguide/Coax_CylinderCoords.m new file mode 100644 index 0000000..5e4606c --- /dev/null +++ b/openEMS/matlab/examples/waveguide/Coax_CylinderCoords.m @@ -0,0 +1,190 @@ +% +% EXAMPLE / waveguide / coaxial cable using cylindrical coordinates +% +% This example demonstrates how to: +% - use cylindrical coordinates +% - setup a coaxial waveguide +% - use analytic functions for waveguide excitations and voltage/current +% calculations +% +% +% Tested with +% - Matlab 2009b +% - openEMS v0.0.17 +% +% (C) 2010 Thorsten Liebig <thorsten.liebig@uni-due.de> + +close all +clear +clc + +%% switches & options... +postprocessing_only = 0; +use_pml = 0; % use pml boundaries instead of mur +openEMS_opts = ''; +% openEMS_opts = [openEMS_opts ' --disable-dumps']; + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +numTS = 1e5; %number of timesteps +length = 1000; %length of the waveguide +unit = 1e-3; %drawing unit used +coax_rad_i = 100; %inner radius +coax_rad_a = 230; %outer radius +mesh_res = [10 nan 10]; %desired mesh resolution +N_alpha = 71; %mesh lines in azimuth direction + +physical_constants; + +%excitation +f0 = 0.5e9; +epsR = 1; + +%% create sim path %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Sim_Path = 'tmp'; +Sim_CSX = 'coax.xml'; + +if (postprocessing_only==0) + [status, message, messageid] = rmdir(Sim_Path,'s'); + [status, message, messageid] = mkdir(Sim_Path); +end + +%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitCylindricalFDTD(numTS,1e-5); +FDTD = SetGaussExcite(FDTD,f0,f0); +BC = {'PEC','PEC','PEC','PEC','PEC','MUR'}; +if (use_pml>0) + BC = {'PEC','PEC','PEC','PEC','PEC','PML_8'}; +end +FDTD = SetBoundaryCond(FDTD,BC); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX('CoordSystem',1); +mesh.x = coax_rad_i : mesh_res(1) : coax_rad_a; +mesh.y = linspace(0,2*pi,N_alpha); +mesh.z = 0 : mesh_res(3) : length; +CSX = DefineRectGrid(CSX, unit, mesh); + +%% material +CSX = AddMaterial(CSX,'fill'); +CSX = SetMaterialProperty(CSX,'fill','Epsilon',epsR); +start = [mesh.x(1) mesh.y(1) 0]; +stop = [mesh.x(end) mesh.y(end) length]; +CSX = AddBox(CSX,'fill',0 ,start,stop); + +%% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = AddExcitation(CSX,'excite',0,[1 0 0]); +weight{1} = '1/rho'; +weight{2} = 0; +weight{3} = 0; +CSX = SetExcitationWeight(CSX, 'excite', weight ); +start = [coax_rad_i mesh.y(1) 0]; +stop = [coax_rad_a mesh.y(end) 0]; +CSX = AddBox(CSX,'excite',0 ,start,stop); + +%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = AddDump(CSX,'Et_','DumpMode',0); +start = [mesh.x(1) , 0 , mesh.z(1)]; +stop = [mesh.x(end) , 0 , mesh.z(end)]; +CSX = AddBox(CSX,'Et_',0 , start,stop); + +CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',0); +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_a 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_a 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_a); +start = [ 0 mesh.z(1) mesh.z(9) ]; +stop = [ mid mesh.z(end) mesh.z(9) ]; +CSX = AddBox(CSX,'it1a', 0 ,start,stop); +CSX = AddProbe(CSX,'it1b',1); +start = [ 0 mesh.z(1) mesh.z(10) ]; +stop = [ mid mesh.z(end) mesh.z(10) ]; +CSX = AddBox(CSX,'it1b', 0 ,start,stop); + +CSX = AddProbe(CSX,'it2a',1); +start = [ 0 mesh.z(1) mesh.z(end-11) ]; +stop = [ mid mesh.z(end) mesh.z(end-11) ]; +CSX = AddBox(CSX,'it2a', 0 ,start,stop); +CSX = AddProbe(CSX,'it2b',1); +start = [ 0 mesh.z(1) mesh.z(end-10) ]; +stop = [ mid mesh.z(end) mesh.z(end-10) ]; +CSX = AddBox(CSX,'it2b', 0 ,start,stop); + +%% Write openEMS +if (postprocessing_only==0) + WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,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); + +%% plot voltages +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_a/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; +plot(freq,20*log10(abs(s11)),'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 +figure() +ZL = uf1./if1; +plot(freq,real(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,ZL_a,'g-.','Linewidth',2); +legend('\Re\{ZL\}','\Im\{ZL\}','ZL-analytic','Location','Best'); diff --git a/openEMS/matlab/examples/waveguide/Coax_Cylindrical_MG.m b/openEMS/matlab/examples/waveguide/Coax_Cylindrical_MG.m new file mode 100644 index 0000000..84a1668 --- /dev/null +++ b/openEMS/matlab/examples/waveguide/Coax_Cylindrical_MG.m @@ -0,0 +1,155 @@ +close all +clear +clc + +%example for an cylindrical mesh, modeling a coaxial cable +% this example is using a multi-grid approach + + +%% setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Settings = []; +Settings.LogFile = 'openEMS.log'; + +physical_constants + +f0 = 0.5e9; +epsR = 1; %material filling + +length = 1000; +port_dist = length/2; +rad_i = 10; %inner radius +rad_a = 200; %outer radius +partial = 0.5; %e.g. 0.5 means only one half of a coax, should be <1 or change boundary cond. +max_mesh = 10 / sqrt(epsR); +max_alpha = max_mesh; +N_alpha = ceil(rad_a * 2*pi * partial / max_alpha); + +%make it even... +N_alpha = N_alpha + mod(N_alpha,2); +%make sure it is multiple of 4, needed for 2 multi-grid steps +N_alpha = ceil((N_alpha)/4) *4 + 1; + +openEMS_opts = ''; +% openEMS_opts = [openEMS_opts ' --disable-dumps']; +% openEMS_opts = [openEMS_opts ' --debug-material']; +% openEMS_opts = [openEMS_opts ' --numThreads=1']; + +def_refSimu = 0; % do a reference simulation without the multi-grid + +%% setup done %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +if (def_refSimu>0) + Sim_Path = 'tmp_ref'; +else + Sim_Path = 'tmp'; +end +Sim_CSX = 'coax.xml'; + +if (exist(Sim_Path,'dir')) + rmdir(Sim_Path,'s'); +end +mkdir(Sim_Path); + +%setup FDTD parameter +if (def_refSimu>0) + FDTD = InitCylindricalFDTD(1e5,1e-5,'OverSampling',5 ); +else + FDTD = InitCylindricalFDTD(1e5,1e-5,'OverSampling',5 ,'MultiGrid','60,120'); +end +FDTD = SetGaussExcite(FDTD,f0,f0); +BC = [0 0 1 1 2 2]; +FDTD = SetBoundaryCond(FDTD,BC); + +mesh_res = [max_mesh 2*pi*partial/N_alpha max_mesh]; + +%setup CSXCAD geometry +CSX = InitCSX(); +mesh.x = SmoothMeshLines([rad_i rad_a],mesh_res(1)); +mesh.y = linspace(-pi*partial,pi*partial,N_alpha); +mesh.z = SmoothMeshLines([0 port_dist length],mesh_res(3)); +CSX = DefineRectGrid(CSX, 1e-3,mesh); + +start = [rad_i mesh.y(1) mesh.z(3)]; +stop = [rad_a mesh.y(end) mesh.z(3)]; + +CSX = AddExcitation(CSX,'excite',0,[1 0 0]); +weight{1} = '1/rho'; +weight{2} = 0; +weight{3} = 0; +CSX = SetExcitationWeight(CSX, 'excite', weight ); +CSX = AddBox(CSX,'excite',0 ,start,stop); + + +start = [mesh.x(1) mesh.y(1) mesh.z(1)]; +stop = [mesh.x(end) mesh.y(end) mesh.z(end)]; +CSX = AddMaterial(CSX,'material'); +CSX = SetMaterialProperty(CSX,'material','Epsilon',epsR); +CSX = AddBox(CSX,'material',0 ,start,stop); + +%dump +CSX = AddDump(CSX,'Et_rz_','DumpMode',0); +start = [mesh.x(1) 0 mesh.z(1)]; +stop = [mesh.x(end) 0 mesh.z(end)]; +CSX = AddBox(CSX,'Et_rz_',0 , start,stop); + +CSX = AddDump(CSX,'Ht_rz_','DumpType',1,'DumpMode',0); +CSX = AddBox(CSX,'Ht_rz_',0 , start,stop); + +CSX = AddDump(CSX,'Et_','DumpType',0,'DumpMode',0); +start = [mesh.x(1) mesh.y(1) length/2]; +stop = [mesh.x(end) mesh.y(end) length/2]; +CSX = AddBox(CSX,'Et_',0,start,stop); + +CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',0); +start = [mesh.x(1) mesh.y(1) length/2]; +stop = [mesh.x(end) mesh.y(end) length/2]; +CSX = AddBox(CSX,'Ht_',0,start,stop); + +% voltage calc (take a voltage average to be at the same spot as the +% current calculation) +CSX = AddProbe(CSX,'ut1_1',0); +start = [ rad_i 0 port_dist ];stop = [ rad_a 0 port_dist ]; +CSX = AddBox(CSX,'ut1_1', 0 ,start,stop); +CSX = AddProbe(CSX,'ut1_2',0); +start = [ rad_i 0 port_dist+mesh_res(3) ];stop = [ rad_a 0 port_dist+mesh_res(3) ]; +CSX = AddBox(CSX,'ut1_2', 0 ,start,stop); + +% current calc +CSX = AddProbe(CSX,'it1',1); +mid = 75; +start = [ 0 mesh.y(1) port_dist+mesh_res(3)/2 ];stop = [ mid mesh.y(end) port_dist+mesh_res(3)/2 ]; +CSX = AddBox(CSX,'it1', 0 ,start,stop); + +%% Write openEMS compatoble xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + +RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts, Settings) + +%% +close all +freq = linspace(0,2*f0,201); +UI = ReadUI({'ut1_1','ut1_2','it1'},Sim_Path,freq); +u_f = (UI.FD{1}.val + UI.FD{2}.val)/2; %averaging voltages to fit current +i_f = UI.FD{3}.val / partial; + +% plot(UI.TD{1}.t,UI.TD{1}.val); +% grid on; +% +% figure +% plot(UI.TD{3}.t,UI.TD{3}.val); +% grid on; + +%plot Z_L compare +figure +ZL = Z0/2/pi/sqrt(epsR)*log(rad_a/rad_i); %analytic line-impedance of a coax +plot(UI.FD{1}.f,ZL*ones(size(u_f)),'g','Linewidth',3); +hold on; +grid on; +Z = u_f./i_f; +plot(UI.FD{1}.f,real(Z),'k--','Linewidth',2); +plot(UI.FD{1}.f,imag(Z),'r-','Linewidth',2); +xlim([0 2*f0]); +legend('Z_L - analytic','\Re\{Z\} - FDTD','\Im\{Z\} - FDTD','Location','Best'); + + diff --git a/openEMS/matlab/examples/waveguide/Rect_Waveguide.m b/openEMS/matlab/examples/waveguide/Rect_Waveguide.m new file mode 100644 index 0000000..a9601ad --- /dev/null +++ b/openEMS/matlab/examples/waveguide/Rect_Waveguide.m @@ -0,0 +1,240 @@ +% +% EXAMPLE / waveguide / Rect_Waveguide +% +% This example demonstrates: +% - waveguide mode excitation +% - waveguide mode matching +% - pml absorbing boundaries +% +% +% Tested with +% - Matlab 2009b +% - openEMS v0.0.17 +% +% (C) 2010 Thorsten Liebig <thorsten.liebig@gmx.de> + +close all +clear +clc + +%% switches +postproc_only = 0; + +%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +physical_constants; +unit = 1e-3; %drawing unit in mm +numTS = 50000; %max. number of timesteps + +% waveguide dimensions +length = 1000; +a = 1000; %waveguide width +b = 600; %waveguide heigth + +%waveguide TE-mode definition +m = 1; +n = 0; + +mesh_res = [10 10 10]; + +%% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%% +f_start = 175e6; +f_stop = 500e6; + +% dump special frequencies to vtk, use paraview (www.paraview.org) to +% animate this dumps over phase +vtk_dump_freq = [200e6 300e6 500e6]; + +freq = linspace(f_start,f_stop,201); + +k = 2*pi*freq/c0; +kc = sqrt((m*pi/a/unit)^2 + (n*pi/b/unit)^2); +fc = c0*kc/2/pi; %cut-off frequency +beta = sqrt(k.^2 - kc^2); %waveguide phase-constant +ZL_a = k * Z0 ./ beta; %analytic waveguide impedance + +disp([' Cutoff frequencies for this mode and wavguide is: ' num2str(fc/1e6) ' MHz']); + +if (f_start<fc) + warning('openEMS:example','f_start is smaller than the cutoff-frequency, this may result in a long simulation... '); +end + +%% mode functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% by David M. Pozar, Microwave Engineering, third edition, page 113 +func_Ex = [num2str( n/b/unit) '*cos(' num2str(m*pi/a) '*x)*sin(' num2str(n*pi/b) '*y)']; +func_Ey = [num2str(-m/a/unit) '*sin(' num2str(m*pi/a) '*x)*cos(' num2str(n*pi/b) '*y)']; + +func_Hx = [num2str(m/a/unit) '*sin(' num2str(m*pi/a) '*x)*cos(' num2str(n*pi/b) '*y)']; +func_Hy = [num2str(n/b/unit) '*cos(' num2str(m*pi/a) '*x)*sin(' num2str(n*pi/b) '*y)']; + +%% define and openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +openEMS_opts = ''; +% openEMS_opts = [openEMS_opts ' --disable-dumps']; +% openEMS_opts = [openEMS_opts ' --debug-material']; +% openEMS_opts = [openEMS_opts ' --engine=basic']; + +Settings = []; +Settings.LogFile = 'openEMS.log'; + +Sim_Path = 'tmp'; +Sim_CSX = 'rect_wg.xml'; + +if (postproc_only==0) + [status, message, messageid] = rmdir(Sim_Path,'s'); + [status, message, messageid] = mkdir(Sim_Path); +end + +%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%% +FDTD = InitFDTD(numTS,1e-5,'OverSampling',6); +FDTD = SetGaussExcite(FDTD,0.5*(f_start+f_stop),0.5*(f_stop-f_start)); +BC = [0 0 0 0 0 3]; +FDTD = SetBoundaryCond(FDTD,BC); + +%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = InitCSX(); +mesh.x = SmoothMeshLines([0 a], mesh_res(1)); +mesh.y = SmoothMeshLines([0 b], mesh_res(2)); +mesh.z = SmoothMeshLines([0 length], mesh_res(3)); +CSX = DefineRectGrid(CSX, unit,mesh); + +%% apply the excitation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +start=[mesh.x(1) mesh.y(1) mesh.z(1) ]; +stop =[mesh.x(end) mesh.y(end) mesh.z(1) ]; +CSX = AddExcitation(CSX,'excite',0,[1 1 0]); +weight{1} = func_Ex; +weight{2} = func_Ey; +weight{3} = 0; +CSX = SetExcitationWeight(CSX,'excite',weight); +CSX = AddBox(CSX,'excite',0 ,start,stop); + +%% voltage and current definitions using the mode matching probes %%%%%%%%% +%port 1 +start = [mesh.x(1) mesh.y(1) mesh.z(15)]; +stop = [mesh.x(end) mesh.y(end) mesh.z(15)]; +CSX = AddProbe(CSX, 'ut1', 10, 1, [], 'ModeFunction',{func_Ex,func_Ey,0}); +CSX = AddBox(CSX, 'ut1', 0 ,start,stop); +CSX = AddProbe(CSX,'it1', 11, 1, [], 'ModeFunction',{func_Hx,func_Hy,0}); +CSX = AddBox(CSX,'it1', 0 ,start,stop); + +%port 2 +start = [mesh.x(1) mesh.y(1) mesh.z(end-15)]; +stop = [mesh.x(end) mesh.y(end) mesh.z(end-15)]; +CSX = AddProbe(CSX, 'ut2', 10, 1, [], 'ModeFunction',{func_Ex,func_Ey,0}); +CSX = AddBox(CSX, 'ut2', 0 ,start,stop); +CSX = AddProbe(CSX,'it2', 11, 1, [], 'ModeFunction',{func_Hx,func_Hy,0}); +CSX = AddBox(CSX,'it2', 0 ,start,stop); + +port_dist = mesh.z(end-15) - mesh.z(15); + +%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +CSX = AddDump(CSX,'Et','FileType',1,'SubSampling','4,4,2'); +start = [mesh.x(1) mesh.y(1) mesh.z(1)]; +stop = [mesh.x(end) mesh.y(end) mesh.z(end)]; +CSX = AddBox(CSX,'Et',0 , start,stop); + +CSX = AddDump(CSX,'Ht','DumpType',1,'FileType',1,'SubSampling','4,4,2'); +CSX = AddBox(CSX,'Ht',0,start,stop); + +%% Write openEMS compatoble xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if (postproc_only==0) + WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + + RunOpenEMS(Sim_Path, Sim_CSX, openEMS_opts, Settings) +end + +%% postproc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +U = ReadUI({'ut1','ut2'},[Sim_Path '/'],freq); +I = ReadUI({'it1','it2'},[Sim_Path '/'],freq); +Exc = ReadUI('et',Sim_Path,freq); + +uf1 = U.FD{1}.val./Exc.FD{1}.val; +uf2 = U.FD{2}.val./Exc.FD{1}.val; +if1 = I.FD{1}.val./Exc.FD{1}.val; +if2 = I.FD{2}.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; +plot(freq,20*log10(abs(s11)),'Linewidth',2); +xlim([freq(1) freq(end)]); +% ylim([-40 5]); +grid on; +hold on; +plot(freq,20*log10(abs(s21)),'r','Linewidth',2); +legend('s11','s21','Location','SouthEast'); +ylabel('s-para (dB)'); +xlabel('freq (Hz)'); + +%% compare analytic and numerical wave-impedance %%%%%%%%%%%%%%%%%%%%%%%%%% +ZL = uf1./if1; +figure() +plot(freq,real(ZL),'Linewidth',2); +hold on; +grid on; +plot(freq,imag(ZL),'r--','Linewidth',2); +plot(freq,ZL_a,'g-.','Linewidth',2); +ylabel('ZL (\Omega)'); +xlabel('freq (Hz)'); +xlim([freq(1) freq(end)]); +legend('\Re(Z_L)','\Im(Z_L)','Z_L analytic','Location','Best'); + +%% beta compare +figure() +da = unwrap(angle(uf1_inc./uf2_inc)) ; +% da = mod(da,2*pi); +beta_12 = (da)/port_dist/unit; +plot(freq,beta_12,'Linewidth',2); +xlim([freq(1) freq(end)]); +xlabel('frequency (Hz)'); +ylabel('\beta (m^{-1})'); +grid on; +hold on; +plot(freq,beta,'g--','Linewidth',2); +legend('\beta-FDTD','\beta-analytic','Location','Best'); + +%% Plot the field dumps %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +dump_file = [Sim_Path '/Et.h5']; +figure() +PlotArgs.slice = {a/2*unit b/2*unit 0}; +PlotArgs.pauseTime=0.01; +PlotArgs.component=0; +PlotArgs.Limit = 'auto'; +PlotHDF5FieldData(dump_file, PlotArgs) + +%% dump frequency to vtk %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% cleanup and create dump folder +vtk_path = [Sim_Path '/vtk']; +[status, message, messageid] = rmdir(vtk_path,'s'); +[status, message, messageid] = mkdir(vtk_path); + +disp('Dumping to vtk files... this may take a minute...') +% define interpolation mesh +mesh_interp{1}=mesh.x * unit; +mesh_interp{2}=b/2 * unit; +mesh_interp{3}=mesh.z * unit; +[field mesh_FD] = ReadHDF5Dump(dump_file,'Interpolation',mesh_interp,'Frequency',vtk_dump_freq); + +% dump animated phase to vtk +for n=1:numel(vtk_dump_freq) + phase = linspace(0,360,21); + phase = phase(1:end-1); + for ph = phase + filename = [vtk_path '/E_xz_f=' num2str(vtk_dump_freq(n)) '_p' num2str(ph) '.vtk']; + Dump2VTK(filename,real(field.FD.values{n}.*exp(1j*ph/180*pi)),mesh_FD,'E-Field'); + end + + filename = [vtk_path '/E_xz_f=' num2str(vtk_dump_freq(n)) '_mag.vtk']; + Dump2VTK(filename,abs(field.FD.values{n}),mesh_FD,'E-Field'); +end + +disp('done... you can open and visualize the vtk-files using Paraview (www.paraview.org)!') |