From d0d2593ab3df0e5d076e5a9e4bad9536969cdf29 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Tue, 7 Feb 2012 10:13:50 +0100 Subject: [PATCH] Tutorials update using new nf2ff calc --- matlab/Tutorials/CRLH_LeakyWaveAnt.m | 115 +++++++++++------------- matlab/Tutorials/Conical_Horn_Antenna.m | 55 +++++------- matlab/Tutorials/Horn_Antenna.m | 52 ++++------- matlab/Tutorials/Simple_Patch_Antenna.m | 48 ++++------ 4 files changed, 110 insertions(+), 160 deletions(-) diff --git a/matlab/Tutorials/CRLH_LeakyWaveAnt.m b/matlab/Tutorials/CRLH_LeakyWaveAnt.m index 8803368..34affc2 100644 --- a/matlab/Tutorials/CRLH_LeakyWaveAnt.m +++ b/matlab/Tutorials/CRLH_LeakyWaveAnt.m @@ -6,7 +6,7 @@ % % Tested with % - Matlab 2011a / Octave 3.4.3 -% - openEMS v0.0.26 +% - openEMS v0.0.27 % % (C) 2011,2012 Thorsten Liebig @@ -44,7 +44,7 @@ f_stop = 6e9; f_rad = (1.9:0.1:4.2)*1e9; -Plot_3D_Rad_Pattern = 0; %this may take a very very long time! > 7h +Plot_3D_Rad_Pattern = 1; %this may take a long time! > 30min %% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%% FDTD = InitFDTD( 20000 ); @@ -99,7 +99,7 @@ portstop = [ -(N_Cells*CRLH.LL)/2, CRLH.LW/2, 0]; 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; @@ -107,7 +107,7 @@ stop = [mesh.x(end) mesh.y(end) mesh.z(end)] - 10*resolution; [CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop); %% write/show/run the openEMS compatible xml-file -Sim_Path = 'tmp'; +Sim_Path = 'tmp_CRLH_LeakyWave'; Sim_CSX = 'CRLH.xml'; [status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory @@ -139,43 +139,52 @@ ylim([-40 2]); drawnow %% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -thetaRange = (0:3:359) - 180; -for n=1:numel(f_rad) - f_res = f_rad(n) - % calculate the far field at phi=0 degrees and at phi=90 degrees - r = 1; % evaluate fields at radius r - disp( 'calculating far field at phi=[0 90] deg...' ); - [E_far_theta{n},E_far_phi{n},Prad(n),Dmax(n)] = AnalyzeNF2FF( Sim_Path, nf2ff, f_res, thetaRange, 0, r ); - toc -end +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, 1, 'Verbose',1); +nf2ff = ReadNF2FF(nf2ff); %% -Dlog=10*log10(Dmax); -figure -thetaRange = (0:3:359) - 180; +% 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( ['radiated power: Prad = ' num2str(Prad(n)) ' Watt']); - disp( ['directivity: Dmax = ' num2str(Dlog(n)) ' dBi'] ); + 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)'] ); - % calculate the e-field magnitude for phi = 0 deg - E_phi0_far{n} = zeros(1,numel(thetaRange)); - for m=1:numel(thetaRange) - E_phi0_far{n}(m) = norm( [E_far_theta{n}(m,1) E_far_phi{n}(m,1)] ); - end + % 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)); - E_phi0_far_log{n} = 20*log10(abs(E_phi0_far{n})/max(abs(E_phi0_far{n}))); - E_phi0_far_log{n} = E_phi0_far_log{n} + Dlog(n); + figure(10) + plot( nf2ff.theta, D_log(:,1) ,line_styles{1+mod(n-1,numel(line_styles))}); + hold on; - % display polar plot - plot( thetaRange, E_phi0_far_log{n} ,'k-' ); - xlabel( 'theta (deg)' ); - ylabel( 'directivity (dBi)'); - grid on; - ylim([-20 10]); - pause(0.5) + figure(11) + plot( nf2ff.theta, D_log(:,2) ,line_styles{1+mod(n-1,numel(line_styles))} ); + hold on; end if (Plot_3D_Rad_Pattern==0) @@ -183,39 +192,17 @@ if (Plot_3D_Rad_Pattern==0) end %% calculate 3D pattern -for n=1:numel(f_rad) - f_res = f_rad(n); - phiRange = 0:3:360; - thetaRange = 0:3:180; - r = 1; % evaluate fields at radius r - disp( 'calculating 3D far field...' ); - [E_far_theta_3D{n},E_far_phi_3D{n}] = AnalyzeNF2FF( Sim_Path, nf2ff, f_res, thetaRange, phiRange, r ); -end +phi = 0:3:360; +theta = 0:3:180; + +disp( 'calculating 3D far field pattern...' ); +nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_rad, theta*pi/180, phi*pi/180, 1, 'Verbose',2); +nf2ff = ReadNF2FF(nf2ff); %% -figure +disp( 'dumping 3D far field pattern to vtk, use Paraview to visualize...' ); for n=1:numel(f_rad) - f_res = f_rad(n); - - E_far_3D{n} = sqrt( abs(E_far_theta_3D{n}).^2 + abs(E_far_phi_3D{n}).^2 ); - E_far_normalized_3D{n} = E_far_3D{n} / max(E_far_3D{n}(:)) * max(Dmax); - - [theta,phi] = ndgrid(thetaRange/180*pi,phiRange/180*pi); - x = E_far_normalized_3D{n} .* sin(theta) .* cos(phi); - y = E_far_normalized_3D{n} .* sin(theta) .* sin(phi); - z = E_far_normalized_3D{n} .* cos(theta); - surf( x,y,z, E_far_normalized_3D{n},'EdgeColor','none'); - caxis([0 max(Dmax)]); - axis equal - xlabel( 'x' ); - xlim([-6 6]); - ylabel( 'y' ); - ylim([-6 6]); - zlabel( 'z' ); - zlim([-4 10]); - title(['f=' num2str(f_res*1e-9,3) 'GHz - D=' num2str(Dlog(n),3) 'dBi'],'FontSize',12) - pause(0.5) - - DumpFF2VTK( [Sim_Path '/FF_Pattern_' int2str(f_res/1e6) 'MHz.vtk'],E_far_normalized_3D,thetaRange,phiRange,1e-3); + E_far_normalized_3D = nf2ff.E_norm{n} / max(max(nf2ff.E_norm{n})) * nf2ff.Dmax(n); + DumpFF2VTK( [Sim_Path '/FF_Pattern_' int2str(f_rad(n)/1e6) 'MHz.vtk'],E_far_normalized_3D,theta,phi,1e-3); end diff --git a/matlab/Tutorials/Conical_Horn_Antenna.m b/matlab/Tutorials/Conical_Horn_Antenna.m index d28bdb6..822aaed 100644 --- a/matlab/Tutorials/Conical_Horn_Antenna.m +++ b/matlab/Tutorials/Conical_Horn_Antenna.m @@ -6,7 +6,7 @@ % % Tested with % - Matlab 2011a / Octave 3.4.3 -% - openEMS v0.0.26 +% - openEMS v0.0.27 % % (C) 2011,2012 Thorsten Liebig @@ -70,7 +70,7 @@ if (f_start @@ -161,10 +161,10 @@ CSX = AddBox(CSX,'it1', 0 ,start,stop); %% 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, [1 1 1 1 0 1]); +[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop, 'Directions', [1 1 1 1 0 1]); %% prepare simulation folder -Sim_Path = 'tmp'; +Sim_Path = 'tmp_Horn_Antenna'; Sim_CSX = 'horn_ant.xml'; [status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory @@ -205,58 +205,44 @@ drawnow % calculate the far field at phi=0 degrees and at phi=90 degrees thetaRange = (0:2:359) - 180; -r = 1; % evaluate fields at radius r disp( 'calculating far field at phi=[0 90] deg...' ); -[E_far_theta,E_far_phi,Prad,Dmax] = AnalyzeNF2FF( Sim_Path, nf2ff, f0, thetaRange, [0 90], r ); +nf2ff = CalcNF2FF(nf2ff, Sim_Path, f0, thetaRange*pi/180, [0 90]*pi/180); -Dlog=10*log10(Dmax); +Dlog=10*log10(nf2ff.Dmax); G_a = 4*pi*A/(c0/f0)^2; -e_a = Dmax/G_a; +e_a = nf2ff.Dmax/G_a; % display some antenna parameter -disp( ['radiated power: Prad = ' num2str(Prad) ' Watt']); +disp( ['radiated power: Prad = ' num2str(nf2ff.Prad) ' Watt']); disp( ['directivity: Dmax = ' num2str(Dlog) ' dBi'] ); disp( ['aperture efficiency: e_a = ' num2str(e_a*100) '%'] ); %% -% calculate the e-field magnitude for phi = 0 deg -E_phi0_far = zeros(1,numel(thetaRange)); -for n=1:numel(thetaRange) - E_phi0_far(n) = norm( [E_far_theta(n,1) E_far_phi(n,1)] ); -end - -E_phi0_far_log = 20*log10(abs(E_phi0_far)/max(abs(E_phi0_far))); -E_phi0_far_log = E_phi0_far_log + Dlog; +% normalized directivity +D_log = 20*log10(nf2ff.E_norm{1}/max(max(nf2ff.E_norm{1}))); +% directivity +D_log = D_log + 10*log10(nf2ff.Dmax); % display polar plot figure -plot( thetaRange, E_phi0_far_log ,'k-' ); +plot( nf2ff.theta, D_log(:,1) ,'k-' ); xlabel( 'theta (deg)' ); ylabel( 'directivity (dBi)'); grid on; hold on; - -% calculate the e-field magnitude for phi = 90 deg -E_phi90_far = zeros(1,numel(thetaRange)); -for n=1:numel(thetaRange) - E_phi90_far(n) = norm([E_far_theta(n,2) E_far_phi(n,2)]); -end - -E_phi90_far_log = 20*log10(abs(E_phi90_far)/max(abs(E_phi90_far))); -E_phi90_far_log = E_phi90_far_log + Dlog; - -% display polar plot -plot( thetaRange, E_phi90_far_log ,'r-' ); +plot( nf2ff.theta, D_log(:,2) ,'r-' ); legend('phi=0','phi=90') +drawnow + %% calculate 3D pattern phiRange = sort( unique( [-180:5:-100 -100:2.5:-50 -50:1:50 50:2.5:100 100:5:180] ) ); thetaRange = sort( unique([ 0:1:50 50:2.:100 100:5:180 ])); -r = 1; % evaluate fields at radius r + disp( 'calculating 3D far field...' ); -[E_far_theta,E_far_phi] = AnalyzeNF2FF( Sim_Path, nf2ff, f0, thetaRange, phiRange, r ); -E_far = sqrt( abs(E_far_theta).^2 + abs(E_far_phi).^2 ); -E_far_normalized = E_far / max(E_far(:)) * Dmax; +nf2ff = CalcNF2FF(nf2ff, Sim_Path, f0, thetaRange*pi/180, phiRange*pi/180, 'Verbose',2,'Outfile','nf2ff_3D.h5'); + +E_far_normalized = nf2ff.E_norm{1} / max(nf2ff.E_norm{1}(:)) * nf2ff.Dmax; [theta,phi] = ndgrid(thetaRange/180*pi,phiRange/180*pi); x = E_far_normalized .* sin(theta) .* cos(phi); diff --git a/matlab/Tutorials/Simple_Patch_Antenna.m b/matlab/Tutorials/Simple_Patch_Antenna.m index a482c17..1189469 100644 --- a/matlab/Tutorials/Simple_Patch_Antenna.m +++ b/matlab/Tutorials/Simple_Patch_Antenna.m @@ -6,7 +6,7 @@ % % Tested with % - Matlab 2011a / Octave 3.4.3 -% - openEMS v0.0.26 +% - openEMS v0.0.27 % % (C) 2010-2012 Thorsten Liebig @@ -94,7 +94,7 @@ SimBox = SimBox - max_res * 4; %reduced SimBox size for nf2ff box [CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', -SimBox/2, SimBox/2); %% prepare simulation folder -Sim_Path = 'tmp'; +Sim_Path = 'tmp_Patch_Ant'; Sim_CSX = 'patch_ant.xml'; [status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory @@ -107,12 +107,12 @@ WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX ); CSXGeomPlot( [Sim_Path '/' Sim_CSX] ); %% run openEMS -RunOpenEMS( Sim_Path, Sim_CSX ); +RunOpenEMS( Sim_Path, Sim_CSX); %% 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) +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 @@ -150,44 +150,28 @@ f_res = freq(f_res_ind); % calculate the far field at phi=0 degrees and at phi=90 degrees thetaRange = (0:2:359) - 180; -r = 1; % evaluate fields at radius r +phiRange = (0:2:359) - 180; disp( 'calculating far field at phi=[0 90] deg...' ); -[E_far_theta,E_far_phi,Prad,Dmax] = AnalyzeNF2FF( Sim_Path, nf2ff, f_res, thetaRange, [0 90], r ); -Dlog=10*log10(Dmax); +nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, thetaRange*pi/180, [0 90]*pi/180); % display power and directivity -disp( ['radiated power: Prad = ' num2str(Prad) ' Watt']); -disp( ['directivity: Dmax = ' num2str(Dlog) ' dBi'] ); -disp( ['efficiency: nu_rad = ' num2str(100*Prad./real(P_in(f_res_ind))) ' %']); +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))) ' %']); -% calculate the e-field magnitude for phi = 0 deg -E_phi0_far = zeros(1,numel(thetaRange)); -for n=1:numel(thetaRange) - E_phi0_far(n) = norm( [E_far_theta(n,1) E_far_phi(n,1)] ); -end - -E_phi0_far_log = 20*log10(abs(E_phi0_far)/max(abs(E_phi0_far))); -E_phi0_far_log = E_phi0_far_log + Dlog; +% normalized directivity +D_log = 20*log10(nf2ff.E_norm{1}/max(max(nf2ff.E_norm{1}))); +% directivity +D_log = D_log + 10*log10(nf2ff.Dmax); % display polar plot figure -plot( thetaRange, E_phi0_far_log ,'k-' ); +plot( nf2ff.theta, D_log(:,1) ,'k-' ); xlabel( 'theta (deg)' ); ylabel( 'directivity (dBi)'); grid on; hold on; - -% calculate the e-field magnitude for phi = 90 deg -E_phi90_far = zeros(1,numel(thetaRange)); -for n=1:numel(thetaRange) - E_phi90_far(n) = norm([E_far_theta(n,2) E_far_phi(n,2)]); -end - -E_phi90_far_log = 20*log10(abs(E_phi90_far)/max(abs(E_phi90_far))); -E_phi90_far_log = E_phi90_far_log + Dlog; - -% display polar plot -plot( thetaRange, E_phi90_far_log ,'r-' ); +plot( nf2ff.theta, D_log(:,2) ,'r-' ); legend('phi=0','phi=90')