From 4690e6a5b902defd5505402a64e7d88873aa8b3d Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Thu, 1 Mar 2012 13:52:15 +0100 Subject: [PATCH] matlab Tutorials: update and new helical antenna tutorial --- matlab/CreateNF2FFBox.m | 1 + matlab/Tutorials/CRLH_LeakyWaveAnt.m | 93 ++++---------- matlab/Tutorials/Helical_Antenna.m | 185 +++++++++++++++++++++++++++ 3 files changed, 211 insertions(+), 68 deletions(-) create mode 100644 matlab/Tutorials/Helical_Antenna.m diff --git a/matlab/CreateNF2FFBox.m b/matlab/CreateNF2FFBox.m index e14df50..1a87b37 100644 --- a/matlab/CreateNF2FFBox.m +++ b/matlab/CreateNF2FFBox.m @@ -21,6 +21,7 @@ function [CSX nf2ff] = CreateNF2FFBox(CSX, name, start, stop, varargin) % % example: % see Tutorials/Simple_Patch_Antenna.m +% see Tutorials/Helical_Antenna.m % % See also CalcNF2FF % diff --git a/matlab/Tutorials/CRLH_LeakyWaveAnt.m b/matlab/Tutorials/CRLH_LeakyWaveAnt.m index 6b4c5c4..b973b94 100644 --- a/matlab/Tutorials/CRLH_LeakyWaveAnt.m +++ b/matlab/Tutorials/CRLH_LeakyWaveAnt.m @@ -43,9 +43,9 @@ Air_Spacer = 30000; f_start = 1e9; f_stop = 6e9; -f_rad = (1.9:0.1:4.2)*1e9; - -Plot_3D_Rad_Pattern = 1; %this may take a long time! > 30min +% frequencies to calculate the 3D radiation pattern +f_rad = (1.9:0.05:4.2)*1e9; +nf2ff_resolution = c0/max(f_rad)/unit/15; %% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%% FDTD = InitFDTD(100000, 1e-3); @@ -101,11 +101,10 @@ portstart = [ feed_length+(N_Cells*CRLH.LL)/2 , -CRLH.LW/2, substratelines(end)] portstop = [ +(N_Cells*CRLH.LL)/2, CRLH.LW/2, 0]; [CSX,portstruct{2}] = AddMSLPort( CSX, 999, 2, 'PEC', portstart, portstop, 0, [0 0 -1], 'MeasPlaneShift', feed_length/2, 'Feed_R', 50 ); - %% nf2ff calc start = [mesh.x(1) mesh.y(1) mesh.z(1) ] + 10*resolution; stop = [mesh.x(end) mesh.y(end) mesh.z(end)] - 10*resolution; -[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop); +[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop, 'OptResolution', nf2ff_resolution); %% write/show/run the openEMS compatible xml-file Sim_Path = 'tmp_CRLH_LeakyWave'; @@ -139,72 +138,30 @@ ylim([-40 2]); drawnow -%% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -theta = (0:3:359) - 180; -phi = [0 90]; - -disp( 'calculating far field at phi=[0 90] deg...' ); - -nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_rad, theta*pi/180, phi*pi/180, 'Verbose',1); - -%% -% prepare figures -figure(10) -hold on; -grid on; -xlabel( 'theta (deg)' ); -ylabel( 'directivity (dBi)'); -title('phi = 0°'); -ylim([-20 10]); -figure(11) -hold on; -grid on; -xlabel( 'theta (deg)' ); -ylabel( 'directivity (dBi)'); -title('phi = 90°'); -ylim([-20 10]); -line_styles = {'b-','g:','r-.','c--','m-','y:','k-.'}; - -for n=1:numel(f_rad) - f_res = f_rad(n); - - % display power and directivity - disp( ['frequency: f = ' num2str(f_res/1e9) ' GHz']); - disp( ['radiated power: Prad = ' num2str(nf2ff.Prad(n)) ' Watt']); - disp( ['directivity: Dmax = ' num2str(nf2ff.Dmax(n)) ' (' num2str(10*log10(nf2ff.Dmax(n))) ' dBi)'] ); - - % normalized directivity - D_log = 20*log10(nf2ff.E_norm{n}/max(max(nf2ff.E_norm{n}))); - % directivity - D_log = D_log + 10*log10(nf2ff.Dmax(n)); - - figure(10) - plot( nf2ff.theta, D_log(:,1) ,line_styles{1+mod(n-1,numel(line_styles))}); - hold on; - - figure(11) - plot( nf2ff.theta, D_log(:,2) ,line_styles{1+mod(n-1,numel(line_styles))} ); - hold on; -end - -%% -figure() -plot(f_rad,nf2ff.Dmax,'b-*','Linewidth',2) -grid on -xlabel( 'frequency' ); -ylabel( 'directivity (dBi)'); - -%% -if (Plot_3D_Rad_Pattern==0) - return -end - %% calculate 3D pattern -phi = 0:3:360; -theta = 0:3:180; +phi = 0:2:360; +theta = 0:2:180; disp( 'calculating 3D far field pattern...' ); -nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_rad, theta*pi/180, phi*pi/180,'Verbose',2, 'Outfile','3D_Pattern.h5'); +nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_rad, theta*pi/180, phi*pi/180, 'Outfile','3D_Pattern.h5', 'Mode', 1,'Verbose',1); + +%% +P_in = 0.5*port{1}.uf.inc .* conj( port{1}.if.inc ); % accepted antenna feed power +P_in = interp1(f, P_in, f_rad); + +figure() + +[AX,H1,H2] = plotyy(f_rad/1e9,nf2ff.Dmax',f_rad/1e9,100*nf2ff.Prad'./real(P_in),'plot'); +grid on +xlabel( 'frequency (GHz)' ); +set(get(AX(1),'Ylabel'),'String','directivity (dBi)') +set(get(AX(2),'Ylabel'),'String','radiation efficiency (%)') +set(H1,'Linewidth',2) +set(H2,'Linewidth',2) +set(H1,'Marker','*') +set(H2,'Marker','s') + +drawnow %% disp( 'dumping 3D far field pattern to vtk, use Paraview to visualize...' ); diff --git a/matlab/Tutorials/Helical_Antenna.m b/matlab/Tutorials/Helical_Antenna.m new file mode 100644 index 0000000..8208c40 --- /dev/null +++ b/matlab/Tutorials/Helical_Antenna.m @@ -0,0 +1,185 @@ +% +% Tutorials / helical antenna +% +% Describtion at: +% http://openems.de/index.php/Tutorial:_Helical_Antenna +% +% Tested with +% - Matlab 2011a / Octave 3.4.3 +% - openEMS v0.0.27 +% +% (C) 2012 Thorsten Liebig + +close all +clear +clc + +post_proc_only = 0; + +close all + +%% setup the simulation +physical_constants; +unit = 1e-3; % all length in mm + +f0 = 2.4e9; % center frequency, frequency of interest! +lambda0 = round(c0/f0/unit); % wavelength in mm +fc = 0.5e9; % 20 dB corner frequency + +Helix.radius = 20; % --> diameter is ~ lambda/pi +Helix.turns = 10; % --> expected gain is G ~ 4 * 10 = 40 (16dBi) +Helix.pitch = 30; % --> pitch is ~ lambda/4 +Helix.wire_rad = 1; + +gnd.radius = lambda0/2; + +% feeding +feed.width = 2; %feeding port width +feed.heigth = 2; +feed.R = 120; %feed impedance + +% size of the simulation box +SimBox = [1 1 1.5]*2*lambda0; + +%% setup FDTD parameter & excitation function +FDTD = InitFDTD( 30000 ); +FDTD = SetGaussExcite( FDTD, f0, fc ); +BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'PML_8'}; % boundary conditions +FDTD = SetBoundaryCond( FDTD, BC ); + +%% setup CSXCAD geometry & mesh +max_res = floor(c0 / (f0+fc) / unit / 20); % cell size: lambda/20 +CSX = InitCSX(); + +mesh.x = [-SimBox(1)/2-gnd.radius -Helix.radius:Helix.wire_rad:Helix.radius SimBox(1)/2+gnd.radius]; +mesh.x = SmoothMeshLines( mesh.x, max_res, 1.4); % create a smooth mesh between specified fixed mesh lines + +mesh.y = mesh.x; + +mesh.z = unique([-SimBox(3)/2 0:Helix.wire_rad:(Helix.turns*Helix.pitch+feed.heigth+Helix.wire_rad) (feed.heigth+Helix.wire_rad+Helix.turns*Helix.pitch)+SimBox(3)/2 ]); +mesh.z = SmoothMeshLines( mesh.z, max_res, 1.4 ); + +CSX = DefineRectGrid( CSX, unit, mesh ); + +%% create helix using the wire primitive +CSX = AddMetal( CSX, 'helix' ); % create a perfect electric conductor (PEC) + +ang = linspace(0,2*pi,21); +coil_x = Helix.radius*cos(ang); +coil_y = Helix.radius*sin(ang); +coil_z = ang/2/pi*Helix.pitch; + +helix.x=[]; +helix.y=[]; +helix.z=[]; +zpos = feed.heigth+Helix.wire_rad; +for n=0:Helix.turns-1 + helix.x = [helix.x coil_x]; + helix.y = [helix.y coil_y]; + helix.z = [helix.z coil_z+zpos]; + zpos = zpos + Helix.pitch; +end +clear p +p(1,:) = helix.x; +p(2,:) = helix.y; +p(3,:) = helix.z; +CSX = AddWire(CSX, 'helix', 0, p, Helix.wire_rad); +start = [Helix.radius-feed.width/2 -feed.width/2 feed.heigth]; +stop = [Helix.radius+feed.width/2 +feed.width/2 feed.heigth+2*Helix.wire_rad]; +CSX = AddBox(CSX,'helix',0,start,stop); + +%% create ground (same size as substrate) +CSX = AddMetal( CSX, 'gnd' ); % create a perfect electric conductor (PEC) +start = [0 0 -0.1]; +stop = [0 0 0.1]; +CSX = AddCylinder(CSX,'gnd',10,start,stop,gnd.radius); + +%% apply the excitation & resist as a current source +start = [Helix.radius-feed.width/2 -feed.width/2 0]; +stop = [Helix.radius+feed.width/2 +feed.width/2 feed.heigth]; +[CSX] = AddLumpedPort(CSX, 5 ,1 ,feed.R, start, stop, [0 0 1], 'excite'); + +%%nf2ff calc +start = [mesh.x(11) mesh.y(11) mesh.z(11)]; +stop = [mesh.x(end-10) mesh.y(end-10) mesh.z(end-10)]; +[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop, 'OptResolution', lambda0/15); + +%% prepare simulation folder +Sim_Path = 'tmp_Helical_Ant'; +Sim_CSX = 'Helix_Ant.xml'; + +if (post_proc_only==0) + [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); +end + +%% postprocessing & do the plots +freq = linspace( f0-fc, f0+fc, 501 ); +U = ReadUI( {'port_ut1','et'}, Sim_Path, freq ); % time domain/freq domain voltage +I = ReadUI( 'port_it1', Sim_Path, freq ); % time domain/freq domain current (half time step is corrected) + +% 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 * feed.R); +if_inc = 0.5*(I.FD{1}.val + U.FD{1}.val / feed.R); +uf_ref = U.FD{1}.val - uf_inc; +if_ref = if_inc - I.FD{1}.val; +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*uf_inc .* conj( if_inc ); % accepted antenna feed power + +drawnow + +%% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%find resonance frequncy from s11 +f_res = f0; + +% get accepted antenna power at frequency f0 +P_in_0 = interp1(freq, P_in, f0); + +% calculate the far field at phi=0 degrees and at phi=90 degrees +thetaRange = unique([0:0.5:90 90:180]); +phiRange = (0:2:360) - 180; +disp( 'calculating far field at phi=[0 90] deg...' ); + +nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, thetaRange*pi/180, phiRange*pi/180,'Mode',1,'Outfile','3D_Pattern.h5','Verbose',1); + +theta_HPBW = thetaRange(find(nf2ff.E_norm{1}(:,1)