From 5b8dfba6cf251babc451f5593edded4c669495ff Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Mon, 28 Nov 2011 14:12:16 +0100 Subject: [PATCH] simplified inf dipole example --- matlab/examples/NF2FF/infDipol.m | 136 +++++++++++-------------------- 1 file changed, 46 insertions(+), 90 deletions(-) diff --git a/matlab/examples/NF2FF/infDipol.m b/matlab/examples/NF2FF/infDipol.m index 5ea8a00..ae00183 100644 --- a/matlab/examples/NF2FF/infDipol.m +++ b/matlab/examples/NF2FF/infDipol.m @@ -1,4 +1,3 @@ -function infDipol % % infinitesimal dipole example % @@ -6,110 +5,71 @@ function infDipol close all clear clc -drawnow +postprocessing_only = 0; -postprocessing_only = 1; - - - - -global g -setup - -dipole_orientation = 3; % 1,2,3: x,y,z -CSX = createStructure(dipole_orientation); -if ~postprocessing_only - writeCSX( CSX ); - CSXGeomPlot( [g.Sim_Path '/' g.Sim_CSX] ) - run; -end -postprocess; - - - - - - -function setup -global g physical_constants % setup the simulation -g.drawingunit = 1e-6; % specify everything in um -g.Sim_Path = 'tmp'; -g.Sim_CSX = 'tmp.xml'; +drawingunit = 1e-6; % specify everything in um +Sim_Path = 'tmp'; +Sim_CSX = 'tmp.xml'; -g.f_max = 1e9; -g.lambda = c0/g.f_max; +f_max = 1e9; +lambda = c0/f_max; % setup geometry values -g.dipole_length = g.lambda/50 /g.drawingunit; +dipole_length = lambda/50 /drawingunit; +dipole_orientation = 3; % 1,2,3: x,y,z -function CSX = createStructure(dipole_orientation) -global g -physical_constants CSX = InitCSX(); % create an equidistant mesh -mesh.x = -g.dipole_length*10:g.dipole_length/2:g.dipole_length*10; -mesh.y = -g.dipole_length*10:g.dipole_length/2:g.dipole_length*10; -mesh.z = -g.dipole_length*10:g.dipole_length/2:g.dipole_length*10; -mesh = AddPML( mesh, [8 8 8 8 8 8] ); % add space for PML -CSX = DefineRectGrid( CSX, g.drawingunit, 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 * -g.dipole_length/2; -stop = ex_vector * g.dipole_length/2; +start = ex_vector * -dipole_length/2; +stop = ex_vector * dipole_length/2; CSX = AddExcitation( CSX, 'infDipole', 1, ex_vector ); CSX = AddBox( CSX, 'infDipole', 1, start, stop ); % NFFF contour -s1 = [-4.5, -4.5, -4.5] * g.dipole_length/2; -s2 = [ 4.5, 4.5, 4.5] * g.dipole_length/2; -[CSX g.nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', s1, s2); +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( max_timesteps, 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); + % define openEMS options and start simulation + openEMS_opts = ''; + RunOpenEMS( Sim_Path, Sim_CSX, openEMS_opts ); +end - -function writeCSX(CSX) -global g -% setup FDTD parameters & excitation function -max_timesteps = 2000; -min_decrement = 1e-6; -FDTD = InitFDTD( max_timesteps, min_decrement, 'OverSampling',10 ); -FDTD = SetGaussExcite( FDTD, g.f_max/2, g.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(g.Sim_Path,'s'); -[~,~,~] = mkdir(g.Sim_Path); -WriteOpenEMS([g.Sim_Path '/' g.Sim_CSX],FDTD,CSX); - - - - -function run -global g -% define openEMS options and start simulation -openEMS_opts = ''; -openEMS_opts = [openEMS_opts ' --engine=fastest']; -RunOpenEMS( g.Sim_Path, g.Sim_CSX, openEMS_opts ); - - - - -function postprocess -global g - +%% post processing disp( ' ' ); disp( ' ********************************************************** ' ); disp( ' ' ); @@ -117,15 +77,13 @@ disp( ' ' ); % calculate the far field at phi=0 degrees and at phi=90 degrees thetaRange = 0:2:359; 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( g.Sim_Path, g.nf2ff, g.f_max, thetaRange, [0 90], r ); - +disp( 'calculating far field at phi=[0 90] deg..' ); +[E_far_theta,E_far_phi,Prad,Dmax] = AnalyzeNF2FF( Sim_Path, nf2ff, f_max, thetaRange, [0 90], r ); % display power and directivity disp( ['radiated power: Prad = ' num2str(Prad)] ); disp( ['directivity: Dmax = ' num2str(Dmax)] ); - % calculate the e-field magnitude for phi = 0 deg E_phi0_far = zeros(1,numel(thetaRange)); for n=1:numel(thetaRange) @@ -139,7 +97,6 @@ ylabel( 'theta / deg' ); title( ['electrical far field (V/m) @r=' num2str(r) ' m phi=0 deg'] ); legend( 'e-field magnitude', 'Location', 'BestOutside' ); - % calculate the e-field magnitude for phi = 90 deg E_phi90_far = zeros(1,numel(thetaRange)); for n=1:numel(thetaRange) @@ -153,14 +110,11 @@ ylabel( 'theta / deg' ); title( ['electrical far field (V/m) @r=' num2str(r) ' m phi=90 deg'] ); legend( 'e-field magnitude', 'Location', 'BestOutside' ); - - - % calculate the far field at theta=90 degrees phiRange = 0:2:359; r = 1; % evaluate fields at radius r -disp( 'calculating far field at theta=90 deg...' ); -[E_far_theta,E_far_phi] = AnalyzeNF2FF( g.Sim_Path, g.nf2ff, g.f_max, 90, phiRange, r ); +disp( 'calculating far field at theta=90 deg..' ); +[E_far_theta,E_far_phi] = AnalyzeNF2FF( Sim_Path, nf2ff, f_max, 90, phiRange, r ); E_theta90_far = zeros(1,numel(phiRange)); for n=1:numel(phiRange) @@ -175,14 +129,12 @@ title( ['electrical far field (V/m) @r=' num2str(r) ' m theta=90 deg'] ); legend( 'e-field magnitude', 'Location', 'BestOutside' ); - - % calculate 3D pattern phiRange = 0:15:360; thetaRange = 0:10:180; r = 1; % evaluate fields at radius r disp( 'calculating 3D far field...' ); -[E_far_theta,E_far_phi] = AnalyzeNF2FF( g.Sim_Path, g.nf2ff, g.f_max, thetaRange, phiRange, r ); +[E_far_theta,E_far_phi] = AnalyzeNF2FF( Sim_Path, nf2ff, f_max, thetaRange, phiRange, r ); E_far = sqrt( abs(E_far_theta).^2 + abs(E_far_phi).^2 ); E_far_normalized = E_far / max(E_far(:)); [theta,phi] = ndgrid(thetaRange/180*pi,phiRange/180*pi); @@ -195,3 +147,7 @@ axis equal xlabel( 'x' ); ylabel( 'y' ); zlabel( 'z' ); + +% +DumpFF2VTK([Sim_Path '/FF_pattern.vtk'],E_far_normalized, thetaRange, phiRange); +disp(['view the farfield pattern "' Sim_Path '/FF_pattern.vtk" using paraview' ]);