Imported Upstream version 0.0.34

8 files changed, 1015 insertions, 0 deletions

diff --git a/openEMS/TESTSUITE/combinedtests/Coax.m b/openEMS/TESTSUITE/combinedtests/Coax.m
new file mode 100644
index 0000000..18ed5ab
--- /dev/null
+++ b/openEMS/TESTSUITE/combinedtests/Coax.m
@@ -0,0 +1,158 @@
+function pass = Coax( openEMS_options, options )
+
+physical_constants;
+
+
+ENABLE_PLOTS = 1;
+CLEANUP = 1;        % if enabled and result is PASS, remove simulation folder
+STOP_IF_FAILED = 1; % if enabled and result is FAILED, stop with error
+SILENT = 0;         % 0=show openEMS output
+
+if nargin < 1
+    openEMS_options = '';
+end
+if nargin < 2
+    options = '';
+end
+if any(strcmp( options, 'run_testsuite' ))
+    ENABLE_PLOTS = 0;
+    STOP_IF_FAILED = 0;
+    SILENT = 1;
+end
+
+% LIMITS
+upper_error = 0.03; % max +3%
+lower_error = 0.01; % max -1%
+
+% structure
+length = 1000;
+coax_rad_i  = 100;
+coax_rad_ai = 230;
+coax_rad_aa = 240;
+mesh_res = [5 5 5];
+f_start = 0;
+f_stop = 1e9;
+
+Sim_Path = 'tmp_Coax';
+Sim_CSX = 'coax.xml';
+
+[status,message,messageid]=rmdir(Sim_Path,'s');
+[status,message,messageid]=mkdir(Sim_Path);
+
+%setup FDTD parameter
+FDTD = InitFDTD(5000,1e-6);
+FDTD = SetGaussExcite(FDTD,0,f_stop);
+FDTD = SetBoundaryCond(FDTD,{'PEC','PEC','PEC','PEC','PEC','PML_8'});
+
+%setup CSXCAD geometry
+CSX = InitCSX();
+mesh.x = -2.5*mesh_res(1)-coax_rad_aa : mesh_res(1) : coax_rad_aa+2.5*mesh_res(1);
+mesh.y = mesh.x;
+mesh.z = 0 : mesh_res(3) : length;
+mesh.z = linspace(0,length,numel(mesh.z));
+CSX = DefineRectGrid(CSX, 1e-3,mesh);
+
+% create a perfect electric conductor
+CSX = AddMetal(CSX,'PEC');
+
+%%% coax
+start = [0, 0 , 0];stop = [0, 0 , length];
+CSX = AddCylinder(CSX,'PEC',1 ,start,stop,coax_rad_i); % inner conductor
+CSX = AddCylindricalShell(CSX,'PEC',0 ,start,stop,0.5*(coax_rad_aa+coax_rad_ai),(coax_rad_aa-coax_rad_ai)); % outer conductor
+
+%%% add excitation
+start(3) = 0; stop(3)=mesh_res(1)/2;
+CSX = AddExcitation(CSX,'excite',0,[1 1 0]);
+weight{1} = '(x)/(x*x+y*y)';
+weight{2} = 'y/pow(rho,2)';
+weight{3} = '0';
+CSX = SetExcitationWeight(CSX, 'excite', weight );
+CSX = AddCylindricalShell(CSX,'excite',0 ,start,stop,0.5*(coax_rad_i+coax_rad_ai),(coax_rad_ai-coax_rad_i));
+
+% %dump
+% CSX = AddDump(CSX,'Et_',0,2);
+% start = [mesh.x(1) , 0 , mesh.z(1)];
+% stop = [mesh.x(end) , 0 , mesh.z(end)];
+% CSX = AddBox(CSX,'Et_',0 , start,stop);
+% 
+% CSX = AddDump(CSX,'Ht_',1,2);
+% CSX = AddBox(CSX,'Ht_',0,start,stop);
+
+%voltage calc
+CSX = AddProbe(CSX,'ut1',0);
+start = [ coax_rad_i 0 length/2 ];stop = [ coax_rad_ai 0 length/2 ];
+CSX = AddBox(CSX,'ut1', 0 ,start,stop);
+
+%current calc
+CSX = AddProbe(CSX,'it1',1);
+% mid = 0.5*(coax_rad_i+coax_rad_ai);
+mid = coax_rad_i+3*mesh_res(1);
+start = [ -mid -mid length/2 ];stop = [ mid mid length/2 ];
+CSX = AddBox(CSX,'it1', 0 ,start,stop);
+
+%Write openEMS compatible xml-file
+WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
+
+% show structure
+% CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
+
+% run openEMS
+folder = fileparts( mfilename('fullpath') );
+Settings.LogFile = [folder '/' Sim_Path '/openEMS.log'];
+Settings.Silent = SILENT;
+RunOpenEMS( Sim_Path, Sim_CSX, openEMS_options, Settings );
+UI = ReadUI( {[Sim_Path '/ut1'], [Sim_Path '/it1']} );
+
+
+%
+% analysis
+%
+
+f = UI.FD{2}.f;
+u = UI.FD{1}.val;
+i = UI.FD{2}.val;
+
+f_idx_start = interp1( f, 1:numel(f), f_start, 'nearest' );
+f_idx_stop  = interp1( f, 1:numel(f), f_stop,  'nearest' );
+f = f(f_idx_start:f_idx_stop);
+u = u(f_idx_start:f_idx_stop);
+i = i(f_idx_start:f_idx_stop);
+
+Z = abs(u./i);
+
+% analytic formular for characteristic impedance
+Z0 = sqrt(MUE0/EPS0) * log(coax_rad_ai/coax_rad_i) / (2*pi);
+upper_limit = Z0 * (1+upper_error);
+lower_limit = Z0 * (1-lower_error);
+
+if ENABLE_PLOTS
+    upper = upper_limit * ones(1,size(Z,2));
+    lower = lower_limit * ones(1,size(Z,2));
+    Z0_plot = Z0 * ones(1,size(Z,2));
+    figure
+    plot(f/1e9,[Z;upper;lower])
+    hold on
+    plot(f/1e9,Z0_plot,'m-.','LineWidth',2)
+    hold off
+    xlabel('Frequency (GHz)')
+    ylabel('Impedance (Ohm)')
+    legend( {'sim', 'upper limit', 'lower limit', 'theoretical'} );
+end
+
+pass = check_limits( Z, upper_limit, lower_limit );
+if pass
+    disp( 'combinedtests/Coax.m (characteristic impedance):  pass' );
+else
+    disp( 'combinedtests/Coax.m (characteristic impedance):  * FAILED *' );
+end
+
+
+
+
+if pass && CLEANUP
+    rmdir( Sim_Path, 's' );
+end
+if ~pass && STOP_IF_FAILED
+    error 'test failed';
+end
+
diff --git a/openEMS/TESTSUITE/combinedtests/README b/openEMS/TESTSUITE/combinedtests/README
new file mode 100644
index 0000000..cc29c5b
--- /dev/null
+++ b/openEMS/TESTSUITE/combinedtests/README
@@ -0,0 +1,3 @@
+#
+# These scripts test the full simulator (not single features)
+#
\ No newline at end of file
diff --git a/openEMS/TESTSUITE/combinedtests/cavity.m b/openEMS/TESTSUITE/combinedtests/cavity.m
new file mode 100644
index 0000000..4123a81
--- /dev/null
+++ b/openEMS/TESTSUITE/combinedtests/cavity.m
@@ -0,0 +1,229 @@
+function pass = cavity( openEMS_options, options )
+
+physical_constants;
+
+
+ENABLE_PLOTS = 1;
+CLEANUP = 1;        % if enabled and result is PASS, remove simulation folder
+STOP_IF_FAILED = 1; % if enabled and result is FAILED, stop with error
+SILENT = 0;         % 0=show openEMS output
+
+if nargin < 1
+    openEMS_options = '';
+end
+if nargin < 2
+    options = '';
+end
+if any(strcmp( options, 'run_testsuite' ))
+    ENABLE_PLOTS = 0;
+    STOP_IF_FAILED = 0;
+    SILENT = 1;
+end
+
+% LIMITS - inside
+lower_rel_limit = 1.3e-3;    % -0.13%
+upper_rel_limit = 1.3e-3;    % +0.13%
+lower_rel_limit_TM = 2.5e-3; % -0.25%
+upper_rel_limit_TM = 0;      % +0%
+min_rel_amplitude = 0.6;     % 60%
+min_rel_amplitude_TM = 0.27; % 27%
+
+% LIMITS - outside
+outer_rel_limit = 0.02;
+max_rel_amplitude = 0.17;
+
+
+% structure
+a = 5e-2;
+b = 2e-2;
+d = 6e-2;
+if ~((b<a) && (a<d))
+    error 'correct the dimensions of the cavity'
+end
+
+f_start = 1e9;
+f_stop = 10e9;
+
+Sim_Path = 'tmp_cavity';
+Sim_CSX = 'cavity.xml';
+
+[status,message,messageid]=rmdir(Sim_Path,'s');
+[status,message,messageid]=mkdir(Sim_Path);
+
+%setup FDTD parameter
+FDTD = InitFDTD( 20000,1e-6 );
+FDTD = SetGaussExcite(FDTD,(f_stop-f_start)/2,(f_stop-f_start)/2);
+BC = [0 0 0 0 0 0]; % PEC boundaries
+FDTD = SetBoundaryCond(FDTD,BC);
+
+%setup CSXCAD geometry
+CSX = InitCSX();
+% grid_res = 2e-3;
+% mesh.x = 0:grid_res:a; %linspace(0,a,25);
+% mesh.y = 0:grid_res:b; %linspace(0,b,25);
+% mesh.z = 0:grid_res:d; %linspace(0,d,25);
+mesh.x = linspace(0,a,26);
+mesh.y = linspace(0,b,11);
+mesh.z = linspace(0,d,32);
+CSX = DefineRectGrid(CSX, 1,mesh);
+
+% excitation
+CSX = AddExcitation(CSX,'excite1',0,[1 1 1]);
+p(1,1) = mesh.x(floor(end*2/3));
+p(2,1) = mesh.y(floor(end*2/3));
+p(3,1) = mesh.z(floor(end*2/3));
+p(1,2) = mesh.x(floor(end*2/3)+1);
+p(2,2) = mesh.y(floor(end*2/3)+1);
+p(3,2) = mesh.z(floor(end*2/3)+1);
+CSX = AddCurve( CSX, 'excite1', 0, p );
+ 
+%dump
+% CSX = AddDump(CSX,'Et_',0,2);
+% pos1 = [mesh.x(1) mesh.y(1) mesh.z(1)];
+% pos2 = [mesh.x(end) mesh.y(end) mesh.z(end)];
+% CSX = AddBox(CSX,'Et_',0 , pos1,pos2);
+
+% %dump
+% CSX = AddDump(CSX,'Et2_',0,2);
+% pos1 = [mesh.x(1) mesh.y(1) mesh.z(1)];
+% pos2 = [mesh.x(end) mesh.y(1) mesh.z(end)];
+% CSX = AddBox(CSX,'Et2_',0 , pos1,pos2);
+% 
+% %dump
+% CSX = AddDump(CSX,'Et3_',0,2);
+% pos1 = [mesh.x(1) mesh.y(end-1) mesh.z(1)];
+% pos2 = [mesh.x(end) mesh.y(end-1) mesh.z(end)];
+% CSX = AddBox(CSX,'Et3_',0 , pos1,pos2);
+
+%voltage calc
+CSX = AddProbe(CSX,'ut1x',0);
+pos1 = [mesh.x(floor(end/4)) mesh.y(floor(end/2))   mesh.z(floor(end/5))];
+pos2 = [mesh.x(floor(end/4)+1) mesh.y(floor(end/2)) mesh.z(floor(end/5))];
+CSX = AddBox(CSX,'ut1x', 0 ,pos1,pos2);
+
+CSX = AddProbe(CSX,'ut1y',0);
+pos1 = [mesh.x(floor(end/4)) mesh.y(floor(end/2))   mesh.z(floor(end/5))];
+pos2 = [mesh.x(floor(end/4)) mesh.y(floor(end/2)+1) mesh.z(floor(end/5))];
+CSX = AddBox(CSX,'ut1y', 0 ,pos1,pos2);
+
+CSX = AddProbe(CSX,'ut1z',0);
+pos1 = [mesh.x(floor(end/2)) mesh.y(floor(end/2)) mesh.z(floor(end/5))];
+pos2 = [mesh.x(floor(end/2)) mesh.y(floor(end/2)) mesh.z(floor(end/5)+1)];
+CSX = AddBox(CSX,'ut1z', 0 ,pos1,pos2);
+
+%Write openEMS compatible xml-file
+WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
+
+% run openEMS
+folder = fileparts( mfilename('fullpath') );
+Settings.LogFile = [folder '/' Sim_Path '/openEMS.log'];
+Settings.Silent = SILENT;
+RunOpenEMS( Sim_Path, Sim_CSX, openEMS_options, Settings );
+UI = ReadUI( {[Sim_Path '/ut1x'], [Sim_Path '/ut1y'], [Sim_Path '/ut1z']} );
+
+
+
+%
+% analysis
+%
+
+% remove excitation from time series
+t_start = 7e-10; % FIXME to be calculated
+t_idx_start = interp1( UI.TD{1}.t, 1:numel(UI.TD{1}.t), t_start, 'nearest' );
+for n=1:numel(UI.TD)
+    UI.TD{n}.t = UI.TD{n}.t(t_idx_start:end);
+    UI.TD{n}.val = UI.TD{n}.val(t_idx_start:end);
+    [UI.FD{n}.f,UI.FD{n}.val] = FFT_time2freq( UI.TD{n}.t, UI.TD{n}.val );
+end
+
+
+f = UI.FD{1}.f;
+ux = UI.FD{1}.val;
+uy = UI.FD{2}.val;
+uz = UI.FD{3}.val;
+
+f_idx_start = interp1( f, 1:numel(f), f_start, 'nearest' );
+f_idx_stop  = interp1( f, 1:numel(f), f_stop,  'nearest' );
+f = f(f_idx_start:f_idx_stop);
+ux = ux(f_idx_start:f_idx_stop);
+uy = uy(f_idx_start:f_idx_stop);
+uz = uz(f_idx_start:f_idx_stop);
+
+% analytic formula for resonant wavenumber
+k = @(m,n,l) sqrt( (m*pi/a)^2 + (n*pi/b)^2 + (l*pi/d)^2 );
+f_TE101 = c0/(2*pi) * k(1,0,1);
+f_TE102 = c0/(2*pi) * k(1,0,2);
+f_TE201 = c0/(2*pi) * k(2,0,1);
+f_TE202 = c0/(2*pi) * k(2,0,2);
+f_TM110 = c0/(2*pi) * k(1,1,0);
+f_TM111 = c0/(2*pi) * k(1,1,1);
+
+f_TE = [f_TE101 f_TE102 f_TE201 f_TE202];
+f_TM = [f_TM110 f_TM111];
+
+% calculate frequency limits
+temp = [f_start f_TE f_stop];
+f_outer1 = [];
+f_outer2 = [];
+for n=1:numel(temp)-1
+    f_outer1 = [f_outer1 temp(n) .* (1+outer_rel_limit)];
+    f_outer2 = [f_outer2 temp(n+1) .* (1-outer_rel_limit)];
+end
+
+temp = [f_start f_TM f_stop];
+f_outer1_TM = [];
+f_outer2_TM = [];
+for n=1:numel(temp)-1
+    f_outer1_TM = [f_outer1_TM temp(n) .* (1+outer_rel_limit)];
+    f_outer2_TM = [f_outer2_TM temp(n+1) .* (1-outer_rel_limit)];
+end
+
+
+if ENABLE_PLOTS
+    figure
+    plot(f/1e9,abs(uy))
+    max1 = max(abs(uy));
+    hold on
+    plot( repmat(f_TE,2,1)/1e9, repmat([0; max1],1,numel(f_TE)), 'm-.', 'LineWidth', 2 )
+    plot( (repmat(f_TE,2,1) .* repmat(1-lower_rel_limit,2,numel(f_TE)))/1e9, repmat([0; max1],1,numel(f_TE)), 'r-', 'LineWidth', 1 )
+    plot( (repmat(f_TE,2,1) .* repmat(1+upper_rel_limit,2,numel(f_TE)))/1e9, repmat([0; max1],1,numel(f_TE)), 'r-', 'LineWidth', 1 )
+    plot( (repmat(f_TE,2,1) .* repmat([1-outer_rel_limit;1+outer_rel_limit],1,numel(f_TE)))/1e9, repmat(max1*min_rel_amplitude,2,numel(f_TE)), 'r-', 'LineWidth', 1 ) % freq limits
+    plot( [f_outer1;f_outer2]/1e9, repmat(max1*max_rel_amplitude,2,numel(f_outer1)), 'g-', 'LineWidth', 1 ) % amplitude limits
+    xlabel('Frequency (GHz)')
+    legend( {'u_y','theoretical'} )
+    title( 'TE-modes' )
+
+    figure
+    plot(f/1e9,abs(uz))
+    max1 = max(abs(uz));
+    hold on
+    plot( repmat(f_TM,2,1)/1e9, repmat([0; max1],1,numel(f_TM)), 'm-.', 'LineWidth', 2 )
+    plot( (repmat(f_TM,2,1) .* repmat(1-lower_rel_limit_TM,2,numel(f_TM)))/1e9, repmat([0; max1],1,numel(f_TM)), 'r-', 'LineWidth', 1 )
+    plot( (repmat(f_TM,2,1) .* repmat(1+upper_rel_limit_TM,2,numel(f_TM)))/1e9, repmat([0; max1],1,numel(f_TM)), 'r-', 'LineWidth', 1 )
+    plot( (repmat(f_TM,2,1) .* repmat([1-lower_rel_limit_TM;1+upper_rel_limit_TM],1,numel(f_TM)))/1e9, repmat(max1*min_rel_amplitude_TM,2,numel(f_TM)), 'r-', 'LineWidth', 1 ) % freq limits
+    plot( [f_outer1_TM;f_outer2_TM]/1e9, repmat(max1*max_rel_amplitude,2,numel(f_outer1_TM)), 'g-', 'LineWidth', 1 ) % amplitude limits
+    xlabel('Frequency (GHz)')
+    legend( {'u_z','theoretical'} )
+    title( 'TM-modes' )
+end
+
+pass1 = check_frequency( f, abs(uy), f_TE*(1+upper_rel_limit),    f_TE*(1-lower_rel_limit),    min_rel_amplitude, 'inside' );
+pass2 = check_frequency( f, abs(uz), f_TM*(1+upper_rel_limit_TM), f_TM*(1-lower_rel_limit_TM), min_rel_amplitude_TM, 'inside' );
+pass3 = check_frequency( f, abs(uy), f_outer2, f_outer1, max_rel_amplitude, 'outside' );
+pass4 = check_frequency( f, abs(uz), f_outer2_TM, f_outer1_TM, max_rel_amplitude, 'outside' );
+pass = pass1 && pass2 && pass3 && pass4;
+if pass
+    disp( 'combinedtests/cavity.m (resonance frequency):  pass' );
+else
+    disp( 'combinedtests/cavity.m (resonance frequency):  * FAILED *' );
+end
+
+
+
+
+if pass && CLEANUP
+    rmdir( Sim_Path, 's' );
+end
+if ~pass && STOP_IF_FAILED
+    error 'test failed';
+end
diff --git a/openEMS/TESTSUITE/enginetests/cavity.m b/openEMS/TESTSUITE/enginetests/cavity.m
new file mode 100644
index 0000000..03976fa
--- /dev/null
+++ b/openEMS/TESTSUITE/enginetests/cavity.m
@@ -0,0 +1,190 @@
+function pass = cavity( openEMS_options, options )
+%pass = cavity( openEMS_options, options )
+% 
+% Checks, if different engines produces identical results
+
+CLEANUP = 1;        % if enabled and result is PASS, remove simulation folder
+STOP_IF_FAILED = 1; % if enabled and result is FAILED, stop with error
+global ENABLE_PLOTS;
+ENABLE_PLOTS = 1;
+SILENT = 0;         % 0=show openEMS output
+
+if nargin < 1
+    openEMS_options = '';
+end
+if nargin < 2
+    options = '';
+end
+if any(strcmp( options, 'run_testsuite' ))
+    ENABLE_PLOTS = 0;
+    STOP_IF_FAILED = 0;
+    SILENT = 1;
+end
+% clean openEMS_options
+openEMS_options = regexprep( openEMS_options, '--engine=\w+', '' );
+
+engines = {'--engine=basic' '--engine=sse' '--engine=sse-compressed' '--engine=multithreaded'};
+% engines = [engines {'--engine=sse-compressed-linear' '--engine=multithreaded-linear'}];
+
+global Sim_Path Sim_CSX
+Sim_Path = 'tmp_cavity';
+Sim_CSX = 'cavity.xml';
+
+for n=1:numel(engines)
+    result{n} = sim( [engines{n} ' ' openEMS_options], SILENT );
+end
+
+pass = compare( result, SILENT );
+
+if pass
+    disp( 'enginetests/cavity.m (engine comparison):  pass' );
+else
+    disp( 'enginetests/cavity.m (engine comparison):  * FAILED *' );
+end
+
+if pass && CLEANUP
+    rmdir( Sim_Path, 's' );
+end
+if ~pass && STOP_IF_FAILED
+    error 'test failed'
+end
+
+return
+
+
+function result = sim( openEMS_options, SILENT )
+global Sim_Path Sim_CSX
+physical_constants;
+
+% structure
+a = 5e-2;
+b = 2e-2;
+d = 6e-2;
+if ~((b<a) && (a<d))
+    error 'correct the dimensions of the cavity'
+end
+
+f_start = 1e9;
+f_stop = 10e9;
+
+% prepare simulation dir
+[status,message,messageid] = rmdir(Sim_Path,'s');
+[status,message,messageid] = mkdir(Sim_Path);
+
+% setup FDTD parameter
+FDTD = InitFDTD( 1000, 0 );
+FDTD = SetGaussExcite(FDTD,(f_stop-f_start)/2,(f_stop-f_start)/2);
+BC = {'MUR' 'PML_8' 'PMC' 'PEC' 'PEC' 'PEC'}; % boundaries
+FDTD = SetBoundaryCond(FDTD,BC);
+
+% setup CSXCAD geometry
+CSX = InitCSX();
+mesh.x = linspace(0,a,27);
+mesh.y = linspace(0,b,11);
+mesh.z = linspace(0,d,33);
+CSX = DefineRectGrid(CSX, 1,mesh);
+
+% excitation
+CSX = AddExcitation(CSX,'excite1',0,[1 1 1]);
+p(1,1) = mesh.x(floor(end*2/3));
+p(2,1) = mesh.y(floor(end*2/3));
+p(3,1) = mesh.z(floor(end*2/3));
+p(1,2) = mesh.x(floor(end*2/3)+1);
+p(2,2) = mesh.y(floor(end*2/3)+1);
+p(3,2) = mesh.z(floor(end*2/3)+1);
+CSX = AddCurve( CSX, 'excite1', 0, p );
+
+% probes
+CSX = AddProbe( CSX, 'E_probe', 2 );
+p(1,1) = mesh.x(floor(end*1/3));
+p(2,1) = mesh.y(floor(end*1/3));
+p(3,1) = mesh.z(floor(end*1/3));
+CSX = AddPoint( CSX, 'E_probe', 0, p );
+CSX = AddProbe( CSX, 'H_probe', 3 );
+CSX = AddPoint( CSX, 'H_probe', 0, p );
+
+% material
+CSX = AddMaterial( CSX, 'RO4350B', 'Epsilon', 3.66 );
+start = [mesh.x(3) mesh.y(3) mesh.z(3)];
+stop  = [mesh.x(5) mesh.y(4) mesh.z(6)];
+CSX = AddBox( CSX, 'RO4350B', 100, start, stop );
+
+% dump
+CSX = AddDump( CSX, 'Et', 'DumpType', 0, 'DumpMode', 0, 'FileType', 1 ); % hdf5 E-field dump without interpolation
+pos1 = [mesh.x(1) mesh.y(1) mesh.z(1)];
+pos2 = [mesh.x(end) mesh.y(end) mesh.z(end)];
+CSX = AddBox( CSX, 'Et', 0, pos1, pos2 );
+
+% dump
+CSX = AddDump( CSX, 'Ht', 'DumpType', 1, 'DumpMode', 0, 'FileType', 1 ); % hdf5 H-field dump without interpolation
+pos1 = [mesh.x(1) mesh.y(1) mesh.z(1)];       % should be half a cell more than now
+pos2 = [mesh.x(end) mesh.y(end) mesh.z(end)]; % should be half a cell less than now
+CSX = AddBox( CSX, 'Ht', 0, pos1, pos2 );
+
+% Write openEMS compatible xml-file
+WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
+
+% cd to working dir and run openEMS
+folder = fileparts( mfilename('fullpath') );
+Settings.LogFile = [folder '/' Sim_Path '/openEMS.log'];
+Settings.Silent = SILENT;
+RunOpenEMS( Sim_Path, Sim_CSX, openEMS_options, Settings );
+
+% collect result
+E.mesh = ReadHDF5Mesh( [Sim_Path '/Et.h5'] );
+E.data = ReadHDF5FieldData( [Sim_Path '/Et.h5'] );
+H.mesh = ReadHDF5Mesh( [Sim_Path '/Ht.h5'] );
+H.data = ReadHDF5FieldData( [Sim_Path '/Ht.h5'] );
+result.E = E;
+result.H = H;
+result.probes = ReadUI( {'E_probe','H_probe'}, Sim_Path );
+
+
+
+function pass = compare( results, SILENT )
+pass = 0;
+% n=1: reference simulation
+for n=2:numel(results)
+    % iterate over all simulations
+    EHfields = {'E','H'};
+    for m=1:numel(EHfields)
+        % iterate over all fields (E, H)
+        EHfield = EHfields{m};
+        for o=1:numel(results{1}.(EHfield).data.TD.values)
+            % iterate over all timesteps
+            cmp_result = results{1}.(EHfield).data.TD.values{o} ~= results{n}.(EHfield).data.TD.values{o};
+            if any(cmp_result(:))
+                disp( ['compare error: n=' num2str(n) '  field=' EHfield '  timestep:' num2str(o) '=' results{1}.(EHfield).data.names{o}] );
+                disp( '  coords:' );
+                find( results{1}.(EHfield).data.TD.values{o} ~= results{n}.(EHfield).data.TD.values{o} )
+                return
+            end
+        end
+    end
+    if ~SILENT
+        disp( ['simulation ' num2str(n) ' is identical to simulation 1'] );
+    end
+end
+
+global ENABLE_PLOTS;
+if ENABLE_PLOTS
+    figure
+    l = {};
+    for n=1:numel(results)
+        plot( results{n}.probes.TD{1}.t, results{n}.probes.TD{1}.val );
+        hold all
+        l = [l ['sim ' num2str(n)]];
+    end
+    legend( l );
+
+    figure
+    l = {};
+    for n=1:numel(results)
+        plot( results{n}.probes.TD{2}.t, results{n}.probes.TD{2}.val );
+        hold all
+        l = [l ['sim ' num2str(n)]];
+    end
+    legend( l );
+end
+
+pass = 1;
diff --git a/openEMS/TESTSUITE/helperscripts/check_frequency.m b/openEMS/TESTSUITE/helperscripts/check_frequency.m
new file mode 100644
index 0000000..33a7cca
--- /dev/null
+++ b/openEMS/TESTSUITE/helperscripts/check_frequency.m
@@ -0,0 +1,31 @@
+function pass = check_frequency( f, val, f_upper, f_lower, rel_amplitude, type )
+
+pass = true;
+max1 = max(val);
+
+if numel(f_upper) ~= numel(f_lower)
+    error 'inconsistant vectors'
+end
+
+for n=1:numel(f_upper)
+    f1 = f_lower(n);
+    f2 = f_upper(n);
+    f1_idx = interp1( f, 1:numel(f), f1, 'nearest' );
+%    if f(f1_idx) < f1, f1_idx = f1_idx + 1; end
+    f2_idx = interp1( f, 1:numel(f), f2, 'nearest' );
+%    if f(f2_idx) > f2, f2_idx = f2_idx - 1; end
+
+    if strcmp( type, 'inside' )
+        if max( val(f1_idx:f2_idx) ) < max1 * rel_amplitude
+            pass = false;
+            return
+        end
+    elseif strcmp( type, 'outside' )
+        if max( val(f1_idx:f2_idx) ) > max1 * rel_amplitude
+            pass = false;
+            return
+        end
+    else
+        error 'unsupported operation'
+    end
+end
diff --git a/openEMS/TESTSUITE/helperscripts/check_limits.m b/openEMS/TESTSUITE/helperscripts/check_limits.m
new file mode 100644
index 0000000..c1cba3d
--- /dev/null
+++ b/openEMS/TESTSUITE/helperscripts/check_limits.m
@@ -0,0 +1,22 @@
+function pass = check_limits( Z, upper_limit, lower_limit )
+
+% make row vector
+if size(Z,1) ~= 1
+    Z = Z.';
+end
+
+if numel(upper_limit) == 1
+    upper_limit = upper_limit * ones(1,size(Z,2));
+end
+if numel(lower_limit) == 1
+    lower_limit = lower_limit * ones(1,size(Z,2));
+end
+
+
+pass = 1;
+if any( Z > upper_limit )
+    pass = 0;
+end
+if any( Z < lower_limit )
+    pass = 0;
+end
diff --git a/openEMS/TESTSUITE/probes/fieldprobes.m b/openEMS/TESTSUITE/probes/fieldprobes.m
new file mode 100644
index 0000000..edf7ed2
--- /dev/null
+++ b/openEMS/TESTSUITE/probes/fieldprobes.m
@@ -0,0 +1,324 @@
+function pass = fieldprobes( openEMS_options, options )
+%
+% infinitesimal dipole in free-space
+%
+% E/H-field probes are compared to hdf5 field dumps
+%
+
+pass = 1;
+
+physical_constants;
+
+
+ENABLE_PLOTS = 1;
+CLEANUP = 1;        % if enabled and result is PASS, remove simulation folder
+STOP_IF_FAILED = 1; % if enabled and result is FAILED, stop with error
+VERBOSE = 1;
+SILENT = 0;         % 0=show openEMS output
+
+if nargin < 1
+    openEMS_options = '';
+end
+if nargin < 2
+    options = '';
+end
+if any(strcmp( options, 'run_testsuite' ))
+    ENABLE_PLOTS = 0;
+    STOP_IF_FAILED = 0;
+    SILENT = 1;
+    VERBOSE = 0;
+end
+
+% LIMITS
+limit_max_time_diff  = 1e-13;
+limit_max_amp_diff   = 1e-7;  %relative amplitude difference
+limit_min_e_amp      = 5e-3;
+limit_min_h_amp      = 1e-7;
+
+
+% setup the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+drawingunit = 1e-6; % specify everything in um
+Sim_Path = 'tmp_fieldprobes';
+Sim_CSX = 'tmp.xml';
+
+f_max = 1e9;
+lambda = c0/f_max /drawingunit;
+
+% setup geometry values
+dipole_length = lambda/50;
+
+
+% setup CSXCAD geometry & mesh %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+CSX = InitCSX();
+mesh.x = -dipole_length*20:dipole_length/2:dipole_length*20;
+mesh.y = -dipole_length*20:dipole_length/2:dipole_length*20;
+mesh.z = -dipole_length*20:dipole_length/2:dipole_length*20;
+CSX = DefineRectGrid( CSX, drawingunit, mesh );
+
+% excitation
+CSX = AddExcitation( CSX, 'infDipole', 1, [0 0 1] );
+start = [0, 0, -dipole_length/2];
+stop  = [0, 0, +dipole_length/2];
+CSX = AddBox( CSX, 'infDipole', 1, start, stop );
+
+% NFFF contour
+s1 = [-4.5, -4.5, -4.5] * dipole_length/2;
+s2 = [ 4.5,  4.5,  4.5] * dipole_length/2;
+CSX = AddBox( AddDump(CSX,'Et_xn','DumpType',0,'DumpMode',0,'FileType',1), 'Et_xn', 0, s1, [s1(1) s2(2) s2(3)] );
+CSX = AddBox( AddDump(CSX,'Et_xp','DumpType',0,'DumpMode',0,'FileType',1), 'Et_xp', 0, [s2(1) s1(2) s1(3)], s2 );
+CSX = AddBox( AddDump(CSX,'Et_yn','DumpType',0,'DumpMode',0,'FileType',1), 'Et_yn', 0, s1, [s2(1) s1(2) s2(3)] );
+CSX = AddBox( AddDump(CSX,'Et_yp','DumpType',0,'DumpMode',0,'FileType',1), 'Et_yp', 0, [s1(1) s2(2) s1(3)], s2 );
+CSX = AddBox( AddDump(CSX,'Et_zn','DumpType',0,'DumpMode',0,'FileType',1), 'Et_zn', 0, s1, [s2(1) s2(2) s1(3)] );
+CSX = AddBox( AddDump(CSX,'Et_zp','DumpType',0,'DumpMode',0,'FileType',1), 'Et_zp', 0, [s1(1) s1(2) s2(3)], s2 );
+CSX = AddBox( AddDump(CSX,'Ht_xn','DumpType',1,'DumpMode',0,'FileType',1), 'Ht_xn', 0, s1, [s1(1) s2(2) s2(3)] );
+CSX = AddBox( AddDump(CSX,'Ht_xp','DumpType',1,'DumpMode',0,'FileType',1), 'Ht_xp', 0, [s2(1) s1(2) s1(3)], s2 );
+CSX = AddBox( AddDump(CSX,'Ht_yn','DumpType',1,'DumpMode',0,'FileType',1), 'Ht_yn', 0, s1, [s2(1) s1(2) s2(3)] );
+CSX = AddBox( AddDump(CSX,'Ht_yp','DumpType',1,'DumpMode',0,'FileType',1), 'Ht_yp', 0, [s1(1) s2(2) s1(3)], s2 );
+CSX = AddBox( AddDump(CSX,'Ht_zn','DumpType',1,'DumpMode',0,'FileType',1), 'Ht_zn', 0, s1, [s2(1) s2(2) s1(3)] );
+CSX = AddBox( AddDump(CSX,'Ht_zp','DumpType',1,'DumpMode',0,'FileType',1), 'Ht_zp', 0, [s1(1) s1(2) s2(3)], s2 );
+
+% E-field probes
+coords{1} = [s1(1) 0 0];
+CSX = AddPoint( AddProbe(CSX,'et1',2), 'et1', 0, coords{1} );
+coords{2} = [s2(1) 0 0];
+CSX = AddPoint( AddProbe(CSX,'et2',2), 'et2', 0, coords{2} );
+coords{3} = [0 s1(2) 0];
+CSX = AddPoint( AddProbe(CSX,'et3',2), 'et3', 0, coords{3} );
+coords{4} = [0 s2(2) 0];
+CSX = AddPoint( AddProbe(CSX,'et4',2), 'et4', 0, coords{4} );
+coords{5} = [0 0 s1(3)];
+CSX = AddPoint( AddProbe(CSX,'et5',2), 'et5', 0, coords{5} );
+coords{6} = [0 0 s2(3)];
+CSX = AddPoint( AddProbe(CSX,'et6',2), 'et6', 0, coords{6} );
+
+% H-field probes
+CSX = AddPoint( AddProbe(CSX,'ht1',3), 'ht1', 0, [s1(1) 0 0] );
+CSX = AddPoint( AddProbe(CSX,'ht2',3), 'ht2', 0, [s2(1) 0 0] );
+CSX = AddPoint( AddProbe(CSX,'ht3',3), 'ht3', 0, [0 s1(2) 0] );
+CSX = AddPoint( AddProbe(CSX,'ht4',3), 'ht4', 0, [0 s2(2) 0] );
+CSX = AddPoint( AddProbe(CSX,'ht5',3), 'ht5', 0, [0 0 s1(3)] );
+CSX = AddPoint( AddProbe(CSX,'ht6',3), 'ht6', 0, [0 0 s2(3)] );
+
+
+
+% setup FDTD parameters & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%
+max_timesteps = 10000;
+min_decrement = 1e-6;
+FDTD = InitFDTD( max_timesteps, min_decrement,'OverSampling',10 );
+FDTD = SetGaussExcite( FDTD, 0, f_max );
+BC = [2 2 2 2 2 2];
+FDTD = SetBoundaryCond( FDTD, BC );
+
+% Write openEMS compatible xml-file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+[~,~,~] = rmdir(Sim_Path,'s');
+[~,~,~] = mkdir(Sim_Path);
+WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
+
+% run openEMS
+folder = fileparts( mfilename('fullpath') );
+Settings.LogFile = [folder '/' Sim_Path '/openEMS.log'];
+Settings.Silent = SILENT;
+RunOpenEMS( Sim_Path, Sim_CSX, openEMS_options, Settings );
+
+
+%% POSTPROCESS
+filenames_E = {'Et_xn.h5','Et_xp.h5','Et_yn.h5','Et_yp.h5','Et_zn.h5','Et_zp.h5'};
+filenames_H = {'Ht_xn.h5','Ht_xp.h5','Ht_yn.h5','Ht_yp.h5','Ht_zn.h5','Ht_zp.h5'};
+
+for n=1:numel(filenames_E)
+    Et{n} = ReadHDF5FieldData( [Sim_Path '/' filenames_E{n}] );
+    E_mesh{n} = ReadHDF5Mesh( [Sim_Path '/' filenames_E{n}] );
+    Ht{n} = ReadHDF5FieldData( [Sim_Path '/' filenames_H{n}] );
+    H_mesh{n} = ReadHDF5Mesh( [Sim_Path '/' filenames_H{n}] );
+    Et_probe{n} = load( [Sim_Path '/et' num2str(n)] );
+    Ht_probe{n} = load( [Sim_Path '/ht' num2str(n)] );
+end
+
+if ENABLE_PLOTS
+    close all
+end
+
+%
+% E-fields
+%
+if VERBOSE, disp( 'extracting field components from field dumps...' ); end
+for n=1:6
+    if numel(E_mesh{n}.lines{1}) > 1
+        x_idx = interp1( E_mesh{n}.lines{1}, 1:numel(E_mesh{n}.lines{1}), coords{n}(1), 'nearest' );
+    else
+        x_idx = 1;
+    end
+    if numel(E_mesh{n}.lines{2}) > 1
+        y_idx = interp1( E_mesh{n}.lines{2}, 1:numel(E_mesh{n}.lines{2}), coords{n}(2), 'nearest' );
+    else
+        y_idx = 1;
+    end
+    if numel(E_mesh{n}.lines{3}) > 1
+        z_idx = interp1( E_mesh{n}.lines{3}, 1:numel(E_mesh{n}.lines{3}), coords{n}(3), 'nearest' );
+    else
+        z_idx = 1;
+    end
+    
+    if VERBOSE
+        disp( ['n=' num2str(n) ' coords: (' num2str(E_mesh{n}.lines{1}(x_idx)) ','...
+               num2str(E_mesh{n}.lines{2}(y_idx)) ','...
+               num2str(E_mesh{n}.lines{3}(z_idx)) ') m  indices: ('...
+               num2str(x_idx) ',' num2str(y_idx) ',' num2str(z_idx) ')'] );
+    end
+    
+    field_x = zeros(numel(Et{n}.TD.values),1);
+    field_y = zeros(numel(Et{n}.TD.values),1);
+    field_z = zeros(numel(Et{n}.TD.values),1);
+    for t=1:numel(Et{n}.TD.values)
+        field_x(t) = squeeze(Et{n}.TD.values{t}(x_idx,y_idx,z_idx,1));
+        field_y(t) = squeeze(Et{n}.TD.values{t}(x_idx,y_idx,z_idx,2));
+        field_z(t) = squeeze(Et{n}.TD.values{t}(x_idx,y_idx,z_idx,3));
+    end
+    field_t = reshape( Et{n}.TD.time, [], 1 );
+    
+    % check vector length
+    if numel(field_x) ~= size(Et_probe{n},1)
+        pass = 0;
+        disp( 'probes/fieldprobes.m (vector length):  * FAILED *' );
+        break
+    end
+    
+    % check absolute simulation time
+    if any(abs(field_t - Et_probe{n}(:,1)) > limit_max_time_diff)
+        pass = 0;
+        disp( 'probes/fieldprobes.m (time inconsistant):  * FAILED *' );
+        break
+    end
+    
+    if ENABLE_PLOTS
+        figure
+        subplot(2,3,1);
+        plot( field_t, [field_x Et_probe{n}(:,2)] );
+        subplot(2,3,2);
+        plot( field_t, [field_y Et_probe{n}(:,3)] );
+        subplot(2,3,3);
+        plot( field_t, [field_z Et_probe{n}(:,4)] );
+        subplot(2,3,4);
+        plot( field_t, (field_x - Et_probe{n}(:,2))./field_x );
+        subplot(2,3,5);
+        plot( field_t, (field_y - Et_probe{n}(:,3))./field_y );
+        subplot(2,3,6);
+        plot( field_t, (field_z - Et_probe{n}(:,4))./field_z );
+    end
+    
+    % difference
+    if any( abs( (field_x - Et_probe{n}(:,2))./field_x) > limit_max_amp_diff ) || ...
+       any( abs( (field_y - Et_probe{n}(:,3))./field_y) > limit_max_amp_diff ) || ...
+       any( abs( (field_z - Et_probe{n}(:,4))./field_z) > limit_max_amp_diff )
+        pass = 0;
+        disp( 'probes/fieldprobes.m (amplitudes differ too much):  * FAILED *' );
+        break
+    end
+    
+    % check absolute field strength of z component
+    if max(abs(field_z)) < limit_min_e_amp
+        pass = 0;
+        disp( 'probes/fieldprobes.m (amplitude of z-component too small):  * FAILED *' );
+        break
+    end
+end
+
+%
+% H-fields
+%
+if VERBOSE, disp( 'extracting field components from field dumps...' ); end
+for n=1:6
+    if numel(H_mesh{n}.lines{1}) > 1
+        x_idx = interp1( H_mesh{n}.lines{1}, 1:numel(H_mesh{n}.lines{1}), coords{n}(1), 'nearest' );
+    else
+        x_idx = 1;
+    end
+    if numel(E_mesh{n}.lines{2}) > 1
+        y_idx = interp1( H_mesh{n}.lines{2}, 1:numel(H_mesh{n}.lines{2}), coords{n}(2), 'nearest' );
+    else
+        y_idx = 1;
+    end
+    if numel(E_mesh{n}.lines{3}) > 1
+        z_idx = interp1( H_mesh{n}.lines{3}, 1:numel(H_mesh{n}.lines{3}), coords{n}(3), 'nearest' );
+    else
+        z_idx = 1;
+    end
+    
+    if VERBOSE
+        disp( ['n=' num2str(n) ' coords: (' num2str(E_mesh{n}.lines{1}(x_idx)) ','...
+               num2str(E_mesh{n}.lines{2}(y_idx)) ','...
+               num2str(E_mesh{n}.lines{3}(z_idx)) ') m  indices: ('...
+               num2str(x_idx) ',' num2str(y_idx) ',' num2str(z_idx) ')'] );
+    end
+    
+    field_x = zeros(numel(Ht{n}.TD.values),1);
+    field_y = zeros(numel(Ht{n}.TD.values),1);
+    field_z = zeros(numel(Ht{n}.TD.values),1);
+    for t=1:numel(Ht{n}.TD.values)
+        field_x(t) = squeeze(Ht{n}.TD.values{t}(x_idx,y_idx,z_idx,1));
+        field_y(t) = squeeze(Ht{n}.TD.values{t}(x_idx,y_idx,z_idx,2));
+        field_z(t) = squeeze(Ht{n}.TD.values{t}(x_idx,y_idx,z_idx,3));
+    end
+    field_t = reshape( Ht{n}.TD.time, [], 1 );
+    
+    % check vector length
+    if numel(field_x) ~= size(Ht_probe{n},1)
+        pass = 0;
+        disp( 'probes/fieldprobes.m (vector length):  * FAILED *' );
+        break
+    end
+    
+    % check absolute simulation time
+    if any(abs(field_t - Ht_probe{n}(:,1)) > limit_max_time_diff)
+        pass = 0;
+        disp( 'probes/fieldprobes.m (time inconsistant):  * FAILED *' );
+        break
+    end
+    
+    if ENABLE_PLOTS
+        figure
+        subplot(2,3,1);
+        plot( field_t, [field_x Ht_probe{n}(:,2)] );
+        subplot(2,3,2);
+        plot( field_t, [field_y Ht_probe{n}(:,3)] );
+        subplot(2,3,3);
+        plot( field_t, [field_z Ht_probe{n}(:,4)] );
+        subplot(2,3,4);
+        plot( field_t, (field_x - Ht_probe{n}(:,2))./field_x );
+        subplot(2,3,5);
+        plot( field_t, (field_y - Ht_probe{n}(:,3))./field_y );
+        subplot(2,3,6);
+        plot( field_t, (field_z - Ht_probe{n}(:,4))./field_z );
+    end
+    
+    % difference
+    if any( abs( (field_x - Ht_probe{n}(:,2))./field_x) > limit_max_amp_diff ) || ...
+       any( abs( (field_y - Ht_probe{n}(:,3))./field_y) > limit_max_amp_diff ) || ...
+       any( abs( (field_z - Ht_probe{n}(:,4))./field_z) > limit_max_amp_diff )
+        pass = 0;
+        disp( 'probes/fieldprobes.m (amplitudes differ too much):  * FAILED *' );
+        break
+    end
+    
+    % check absolute field strength of z component
+    if (max(abs(field_x)) < limit_min_h_amp) || (max(abs(field_y)) < limit_min_h_amp)
+        pass = 0;
+        disp( 'probes/fieldprobes.m (amplitude of x- or y-component too small):  * FAILED *' );
+        break
+    end
+end
+
+
+
+
+if pass
+    disp( 'probes/fieldprobes.m:  pass' );
+end
+
+    
+if pass && CLEANUP
+    rmdir( Sim_Path, 's' );
+end
+if ~pass && STOP_IF_FAILED
+    error 'test failed';
+end
diff --git a/openEMS/TESTSUITE/run_testsuite.m b/openEMS/TESTSUITE/run_testsuite.m
new file mode 100644
index 0000000..9f1707b
--- /dev/null
+++ b/openEMS/TESTSUITE/run_testsuite.m
@@ -0,0 +1,58 @@
+%
+% run the testsuite
+%
+
+clc
+clear
+close all
+drawnow
+
+if isOctave
+    confirm_recursive_rmdir(0);
+    page_screen_output(0);      % do not buffer output
+    page_output_immediately(1); % do not buffer output
+end
+
+folder = fileparts( mfilename( 'fullpath' ) );
+cd( folder );
+addpath( [folder filesep 'helperscripts'] );
+
+% openEMS options
+options = {'--engine=multithreaded', '--engine=sse-compressed', '--engine=sse', '--engine=basic'};
+
+for o=1:numel(options)
+
+    disp( [datestr(now) ' *** TESTSUITE started (options: ' options{o} ')'] );
+
+    % now list the tests
+    folders = dir();
+    for f=1:numel(folders)
+        if folders(f).isdir
+            if strcmp(folders(f).name,'.') || strcmp(folders(f).name,'..')
+                continue
+            end
+            if strcmp(folders(f).name,'helperscripts')
+                continue
+            end
+            oldpwd = pwd;
+            cd( folders(f).name );
+            scripts = dir('*.m');
+            for s=1:numel(scripts)
+                if ~scripts(s).isdir
+                    % execute function
+                    disp( [datestr(now) ' executing: ' folders(f).name '/' scripts(s).name] );
+                    [~,fname] = fileparts( scripts(s).name );
+                    if isOctave
+                        fflush(1); % flush stdout
+                    end
+                    pass = feval( fname, options{o}, 'run_testsuite' );
+                end
+            end
+            cd(oldpwd);
+        end
+    end
+end
+
+disp( '***' );
+disp( ['*** ' datestr(now) ' ALL TESTS DONE'] );
+disp( '***' );