1%
2% Tutorials / CylindricalWave_CC
3%
4% Describtion at:
5% http://openems.de/index.php/Tutorial:_2D_Cylindrical_Wave
6%
7% Tested with
8%  - Matlab 2011a/ Octave 4.0
9%  - openEMS v0.0.33
10%
11% (C) 2011-2015 Thorsten Liebig <thorsten.liebig@gmx.de>
12
13close all
14clear
15clc
16
17%% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
18physical_constants
19mesh_res = 10;      %desired mesh resolution
20radius = 2560;      %simulation domain radius
21split = ['80,160,320,640,1280']; %radii to split the mesh into sub-grids
22split_N = 5;        %number of nested sub-grids
23heigth = mesh_res*4;
24
25f0 = 1e9;
26
27exite_offset = 1300;
28excite_angle = 45;
29
30%% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31FDTD = InitFDTD('NrTS', 100000, 'EndCriteria', 1e-4, 'CoordSystem', 1, 'MultiGrid', split);
32FDTD = SetGaussExcite(FDTD,f0,f0/2);
33BC = [0 3 0 0 0 0];             % pml in positive r-direction
34FDTD = SetBoundaryCond(FDTD,BC);
35
36%% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37% 50 mesh lines for the inner most mesh
38% increase the total number of meshlines in alpha direcion for all sub-grids
39N_alpha = 50 * 2^split_N + 1;
40
41CSX = InitCSX('CoordSystem',1);
42mesh.r = SmoothMeshLines([0 radius],mesh_res);
43mesh.a = linspace(-pi,pi,N_alpha);
44mesh.z = SmoothMeshLines([-heigth/2 0 heigth/2],mesh_res);
45CSX = DefineRectGrid(CSX, 1e-3,mesh);
46
47%% add the dipol %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48start = [exite_offset excite_angle/180*pi-0.001 -20];
49stop =  [exite_offset excite_angle/180*pi+0.001  20];
50if (exite_offset==0)
51    start(2) = mesh.a(1);
52    stop(2)  = mesh.a(1);
53end
54CSX = AddExcitation(CSX,'excite',1,[0 0 1]);
55CSX = AddBox(CSX,'excite',0 ,start,stop);
56
57%% define dump boxes... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58start = [mesh.r(1)   mesh.a(1)   0];
59stop =  [mesh.r(end-8) mesh.a(end) 0];
60
61% time domain vtk dump
62CSX = AddDump(CSX,'Et_ra','DumpType',0,'FileType',0,'SubSampling','4,10,1');
63CSX = AddBox(CSX,'Et_ra',0 , start,stop);
64
65% frequency domain hdf5 dump
66CSX = AddDump(CSX,'Ef_ra','DumpType',10,'FileType',1,'SubSampling','2,2,2','Frequency',f0);
67CSX = AddBox(CSX,'Ef_ra',0 , start,stop);
68
69%% write/run the openEMS compatible xml-file
70Sim_Path = 'tmp';
71Sim_CSX = '2D_CC_Wave.xml';
72
73[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
74[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder
75
76WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
77RunOpenEMS(Sim_Path, Sim_CSX);
78
79%%
80disp('use Paraview to visualize the vtk field dump...');
81
82%%
83[field mesh_h5] = ReadHDF5Dump([Sim_Path '/Ef_ra.h5']);
84
85r = mesh_h5.lines{1};
86a = mesh_h5.lines{2};
87a(end+1) = a(1);            %closeup mesh for visualization
88[R A] = ndgrid(r,a);
89X = R.*cos(A);
90Y = R.*sin(A);
91
92Ez = squeeze(field.FD.values{1}(:,:,1,3));
93Ez(:,end+1) = Ez(:,1);      %closeup mesh for visualization
94
95E_max = max(max(abs(Ez)));  %get maximum E_z amplitude
96
97while 1
98    for ph = linspace(0,360,41) %animate phase from 0..360 degree
99        surf(X,Y,real(Ez*exp(1j*ph*pi/180)),'EdgeColor','none')
100        caxis([-E_max E_max]/10)
101        zlim([-E_max E_max])
102        pause(0.3)
103    end
104end
105