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