New upstream version 0.0.35+ds.1

10 files changed, 655 insertions, 52 deletions

diff --git a/openEMS/matlab/DelayFidelity.m b/openEMS/matlab/DelayFidelity.m
new file mode 100644
index 0000000..3cc0ae2
--- /dev/null
+++ b/openEMS/matlab/DelayFidelity.m
@@ -0,0 +1,93 @@
+function [delay, fidelity, nf2ff_out] = DelayFidelity(nf2ff, port, path, weight_theta, weight_phi, theta, phi, f_0, f_c, varargin)

+% [delay, fidelity] = DelayFidelity(nf2ff, port, path, theta, phi, f_lo, f_hi, varargin)

+% 

+% 

+% This function calculates the time delay from the source port to the phase center of the antenna and the fidelity. 

+% The fidelity is the similarity between the excitation pulse and the radiated pulse (normalized scalar product).

+% The resolution of the delay will be equal to or better than ((f_0 + f_c)*Oversampling)^-1 when using Gaussian excitation.

+% Oversampling is an input parameter to InitFDTD. The rows of delay and fidelity correspond to theta and the columns to phi.

+% 

+% input:

+%   nf2ff: return value of CreateNF2FFBox.

+%   port: return value of AddLumpedPort

+%   path: path of the simulation results.

+%   weight_theta: weight if the E_theta component 

+%   weight_phi: eight of the E_phi component

+%     -> with both (possibly complex) parameters any polarization can be examined

+%   theta: theta values to be simulated

+%   phi: phi values to be simulated

+%   f_0: center frequency of SetGaussExcite

+%   f_c: cutoff frequency of SetGaussExcite

+%

+% variable input:

+%   'Center': phase center of the antenna for CalcNF2FF

+%   'Radius': radius for CalcNF2FF

+%   'Mode': mode CalcNF2FF

+%   

+% example:

+%   theta = [-180:10:180] * pi / 180;

+%   phi = [0, 90] * pi / 180;

+%   [delay, fidelity] = DelayFidelity2(nf2ff, port, Sim_Path, sin(tilt), cos(tilt), theta, phi, f_0, f_c, 'Mode', 1);

+%   figure

+%   polar(theta.', delay(:,1) * 3e11); % delay in mm

+%   figure

+%   polar(theta', (fidelity(:,1)-0.95)/0.05); % last 5 percent of fidelity

+% 

+% Author: Georg Michel 

+

+C0 = 299792458;

+center = [0, 0, 0];

+radius = 1;

+nf2ff_mode = 0;

+

+for n=1:2:numel(varargin)

+    if (strcmp(varargin{n},'Center')==1);

+        center = varargin{n+1};

+    elseif (strcmp(varargin{n},'Radius')==1);

+        radius = varargin{n+1};

+    elseif (strcmp(varargin{n},'Mode')==1);

+        nf2ff_mode = varargin{n+1};

+    end

+end

+

+

+port_ut = load(fullfile(path, port.U_filename));

+port_it = load(fullfile(path, port.I_filename));

+dt = port_ut(2,1) - port_ut(1,1);

+fftsize = 2^(nextpow2(size(port_ut)(1)) + 1);

+df = 1 / (dt * fftsize);

+uport = fft(port_ut(:, 2), fftsize)(1:fftsize/2+1);

+iport = fft(port_it(:, 2), fftsize)(1:fftsize/2+1);

+fport = df * (0:fftsize/2);

+f_ind = find(fport > (f_0 - f_c ) & fport < (f_0 + f_c));

+disp(["frequencies: ", num2str(numel(f_ind))]);

+exc_f = uport.' + iport.' * port.Feed_R; %excitation in freq domain

+exc_f(!f_ind) = 0;

+exc_f /= sqrt(exc_f * exc_f'); % normalization (transposing also conjugates)

+

+nf2ff = CalcNF2FF(nf2ff, path, fport(f_ind), theta, phi, ...

+        'Center', center, 'Radius', radius, 'Mode', nf2ff_mode);

+radfield = weight_theta * cell2mat(nf2ff.E_theta) + weight_phi * cell2mat(nf2ff.E_phi); % rows: theta(f1), columns: phi(f1), phi(f2), ...phi(fn)

+radfield = reshape(radfield, [length(nf2ff.theta), length(nf2ff.phi), length(nf2ff.freq)]);

+correction = reshape(exp(-2i*pi*nf2ff.r/C0*nf2ff.freq), 1,1,numel(nf2ff.freq)); %dimensions: theta, phi, frequencies

+radfield = radfield./correction; % correct for radius delay

+% normalize radfield

+radnorm = sqrt(dot(radfield, radfield, 3));

+radfield ./= radnorm;

+

+%initialize radiated field in fully populated frequency domain

+rad_f = zeros([numel(nf2ff.theta), numel(nf2ff.phi), numel(fport)]);

+rad_f(:, :, f_ind) = radfield; % assign selected frequencies

+exc_f = reshape(exc_f, [1,1,numel(exc_f)]); %make exc_f confomant with rad_f

+

+cr_f = rad_f .* conj(exc_f); % calculate cross correlation

+% calculate the cross correlation in time domain (analytic signal)

+cr = ifft(cr_f(:, :, 1:end-1), [], 3) * (numel(fport) -1); % twice the FFT normalization (sqrt^2) because product of two normalized functions 

+%search for the maxiumum of the envelope

+[fidelity, delay_ind] = max(abs(cr), [], 3);

+delay = (delay_ind - 1) * dt * 2; % double time step because of single-sided FFT

+nf2ff_out = nf2ff; %possibly needed for plotting the far field and other things

+disp(["DelayFidelity: delay resolution = ", num2str(dt*2e9), "ns"]);

+return;

+

+

diff --git a/openEMS/matlab/RunOpenEMS.m b/openEMS/matlab/RunOpenEMS.m
index a4947e3..3894adc 100644
--- a/openEMS/matlab/RunOpenEMS.m
+++ b/openEMS/matlab/RunOpenEMS.m
@@ -21,7 +21,7 @@ function RunOpenEMS(Sim_Path, Sim_File, opts, Settings)
 %             --engine=fastest         fastest available engine (default)
 %             --engine=basic           basic FDTD engine
 %             --engine=sse             engine using sse vector extensions
-%             --engine=sse_compressed  engine using compressed operator + sse vector extensions
+%             --engine=sse-compressed  engine using compressed operator + sse vector extensions
 %             --engine=MPI             engine using compressed operator + sse vector extensions + MPI parallel processing
 %             --engine=multithreaded   engine using compressed operator + sse vector extensions + MPI + multithreading
 %         --numThreads=<n>     Force use n threads for multithreaded engine
diff --git a/openEMS/matlab/RunOpenEMS_MPI.m b/openEMS/matlab/RunOpenEMS_MPI.m
index 495f7e5..779b038 100644
--- a/openEMS/matlab/RunOpenEMS_MPI.m
+++ b/openEMS/matlab/RunOpenEMS_MPI.m
@@ -87,7 +87,7 @@ end
 
 if isfield(Settings.MPI,'Hosts')
     disp(['Running remote openEMS_MPI in working dir: ' work_path]);
-    [status]  = system(['mpiexec -host ' HostList ' -n ' int2str(NrProc) ' -wdir ' work_path ' ' Settings.MPI.Binary ' ' Sim_File ' ' opts ' ' append_unix]);
+    [status]  = system(['mpiexec ' Settings.MPI.GlobalArgs ' -host ' HostList ' -n ' int2str(NrProc) ' -wdir ' work_path ' ' Settings.MPI.Binary ' ' Sim_File ' ' opts ' ' append_unix]);
 else
     disp('Running local openEMS_MPI');
     [status]  = system(['mpiexec ' Settings.MPI.GlobalArgs ' -n ' int2str(NrProc) ' ' Settings.MPI.Binary ' ' Sim_File ' ' opts ' ' append_unix]);
diff --git a/openEMS/matlab/Tutorials/CylindricalWave_CC.m b/openEMS/matlab/Tutorials/CylindricalWave_CC.m
index d55c470..b81d2d4 100644
--- a/openEMS/matlab/Tutorials/CylindricalWave_CC.m
+++ b/openEMS/matlab/Tutorials/CylindricalWave_CC.m
@@ -28,7 +28,7 @@ exite_offset = 1300;
 excite_angle = 45;
 
 %% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-FDTD = InitFDTD(100000,1e-4,'CoordSystem',1,'MultiGrid',split);
+FDTD = InitFDTD('NrTS', 100000, 'EndCriteria', 1e-4, 'CoordSystem', 1, 'MultiGrid', split);
 FDTD = SetGaussExcite(FDTD,f0,f0/2);
 BC = [0 3 0 0 0 0];             % pml in positive r-direction
 FDTD = SetBoundaryCond(FDTD,BC);
diff --git a/openEMS/matlab/Tutorials/RadarUWBTutorial.m b/openEMS/matlab/Tutorials/RadarUWBTutorial.m
new file mode 100644
index 0000000..6fd2f5b
--- /dev/null
+++ b/openEMS/matlab/Tutorials/RadarUWBTutorial.m
@@ -0,0 +1,236 @@
+% Tutorial on time delay and signal integrity for radar 

+% and UWB applications

+% 

+% Tested with

+%  - Octave 4.0

+%  - openEMS v0.0.35

+% 

+% Author: Georg Michel, 2016

+

+ clear;

+ close all;

+

+physical_constants;

+

+% --- start of configuration section ---

+

+% In radar and ultrawideband applications it is important to know the

+% delay and fidelity of RF pulses. The delay is the retardation of the

+% signal from the source to the phase center of the antenna. It is

+% composed out of linear delay, dispersion and minimum-phase

+% delay. Dispersion due to waveguides or frequency-dependent

+% permittivity and minimum-phase delay due to resonances will degrade

+% the fidelity which is the normalized similarity between excitation and

+% radiated signal. In this tutorial you can examine the performance of a

+% simple ultrawideband (UWB) monopole. The delay and fidelity of this

+% antenna are calculated and plotted. You can compare these properties

+% in different channels.

+% 

+% The Gaussian excitation is set to the same 3dB bandwidth as the

+% channels of the IEEE 802.15.4 UWB PHY. One exeption is channel4twice

+% which has the double bandwidth of channel 4. It can be seen that the

+% delay is larger and the fidelity is smaller in the vicinity of the

+% (undesired) resonances of the antenna. Note that for a real UWB system

+% the total delay and fidelity result from both the transmitting and

+% receiving antenna or twice the delay and the square of the fidelity

+% for monostatic radar. 

+% 

+% The resolution of the delay will depend on the 'Oversampling'

+% parameter to InitFDTD.  See the description of DelayFidelity

+% 

+% In the configuration section below you can uncomment the respective

+% parameter settings. As an exercise, you can examine how the permittivity 

+% of the substrate influences gain, delay and fidelity.  

+

+

+%suffix = "channel1";

+%f_0 = 3.5e9; % center frequency of the channel

+%f_c = 0.25e9 / 0.3925; % 3dB bandwidth is 0.3925 times 20dB bandwidth for Gaussian excitation

+

+%suffix = "channel2";

+%f_0 = 4.0e9; % center frequency of the channel

+%f_c = 0.25e9 / 0.3925; 

+

+%suffix = "channel3";

+%f_0 = 4.5e9; % center frequency of the channel

+%f_c = 0.25e9 / 0.3925; 

+

+suffix = "channel4";

+f_0 = 4.0e9; % center frequency of the channel

+f_c = 0.5e9 / 0.3925; 

+

+%suffix = "channel5";

+%f_0 = 6.5e9; % center frequency of the channel

+%f_c = 0.25e9 / 0.3925; 

+

+%suffix = "channel7";

+%f_0 = 6.5e9; % center frequency of the channel

+%f_c = 0.5e9 / 0.3925; 

+

+%suffix = "channel4twice"; % this is just to demonstrate the degradation of the fidelity with increasing bandwidth

+%f_0 = 4.0e9; % center frequency of the channel

+%f_c = 1e9 / 0.3925; 

+

+tilt = 45 * pi / 180; % polarization tilt angle against co-polarization (90DEG is cross polarized) 

+

+% --- end of configuration section ---

+

+% path and filename setup

+Sim_Path = 'tmp';

+Sim_CSX = 'uwb.xml';

+

+% properties of the substrate

+substrate.epsR = 4; % FR4

+substrate.height = 0.707;

+substrate.cells = 3; % thickness in cells

+

+% size of the monopole and the gap to the ground plane

+gap = 0.62; % 0.5

+patchsize = 14;

+

+% we will use millimeters

+unit = 1e-3;

+

+% set the resolution for the finer structures, e.g. the antenna gap

+fineResolution = C0 / (f_0 + f_c) / sqrt(substrate.epsR) / unit / 40;

+% set the resolution for the coarser structures, e.g. the surrounding air

+coarseResolution = C0/(f_0 + f_c) / unit / 20;

+

+

+% initialize the CSX structure

+CSX = InitCSX();

+

+% add the properties which are used to model the antenna

+CSX = AddMetal(CSX, 'Ground' );

+CSX = AddMetal(CSX, 'Patch');

+CSX = AddMetal(CSX, 'Line');

+CSX = AddMaterial(CSX, 'Substrate' );

+CSX = SetMaterialProperty(CSX, 'Substrate', 'Epsilon', substrate.epsR);

+

+% define the supstrate and sheet-like primitives for the properties

+CSX = AddBox(CSX, 'Substrate', 1, [-16, -16, -substrate.height], [16, 18, 0]);

+CSX = AddBox(CSX, 'Ground', 2, [-16, -16, -substrate.height], [16, 0, -substrate.height]);

+CSX = AddBox(CSX, 'Line', 2, [-1.15, -16, 0], [1.15, gap, 0]);

+CSX = AddBox(CSX, 'Patch', 2, [-patchsize/2, gap, 0], [patchsize/2, gap + patchsize, 0]);

+

+% setup a mesh

+mesh.x = [];

+mesh.y = [];

+

+% two mesh lines for the metal coatings of teh substrate

+mesh.z = linspace(-substrate.height, 0, substrate.cells +1);

+

+% find optimal mesh lines for the patch and ground, not yes the microstrip line

+mesh = DetectEdges(CSX, mesh, 'SetProperty',{'Patch', 'Ground'}, '2D_Metal_Edge_Res', fineResolution/2);

+

+%replace gap mesh lines which are too close by a single mesh line

+tooclose = find (diff(mesh.y) < fineResolution/4);

+if ~isempty(tooclose)

+  mesh.y(tooclose) = (mesh.y(tooclose) + mesh.y(tooclose+1))/2;

+  mesh.y(tooclose + 1) = [];

+endif

+

+% store the microstrip  edges in a temporary variable

+meshline = DetectEdges(CSX, [], 'SetProperty', 'Line', '2D_Metal_Edge_Res', fineResolution/2);

+% as well as the edges of the substrate (without 1/3 - 2/3 rule)

+meshsubstrate = DetectEdges(CSX, [], 'SetProperty', 'Substrate');

+% add only the x mesh lines of the microstrip 

+mesh.x = [mesh.x meshline.x];

+% and only the top of the substrate, the other edges are covered by the ground plane

+mesh.y = [mesh.y, meshsubstrate.y(end)]; % top of substrate

+

+% for now we have only the edges, now calculate mesh lines inbetween

+mesh = SmoothMesh(mesh, fineResolution);

+

+% add the outer boundary

+mesh.x = [mesh.x -60, 60];

+mesh.y = [mesh.y, -60, 65];

+mesh.z = [mesh.z, -46, 45];

+

+% add coarse mesh lines for the free space

+mesh = SmoothMesh(mesh, coarseResolution);

+

+% define the grid

+CSX = DefineRectGrid( CSX, unit, mesh);

+% and the feeding port

+[CSX, port] = AddLumpedPort( CSX, 999, 1, 50, [-1.15, meshline.y(2), -substrate.height], [1.15, meshline.y(2), 0], [0 0 1], true);

+

+%setup a NF2FF box for the calculation of the far field

+start = [mesh.x(10)     mesh.y(10)     mesh.z(10)];

+stop  = [mesh.x(end-9) mesh.y(end-9) mesh.z(end-9)];

+[CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop);

+

+% initialize the FDTD structure with excitation and open boundary conditions

+FDTD = InitFDTD( 'NrTs', 30000, 'EndCriteria', 1e-5, 'OverSampling', 20);

+FDTD = SetGaussExcite(FDTD, f_0, f_c );

+BC   = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'};

+FDTD = SetBoundaryCond(FDTD, BC );

+

+

+% remove old data, show structure, calculate new data

+[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory

+[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder

+

+% write the data to the working directory

+WriteOpenEMS([Sim_Path '/' Sim_CSX], FDTD, CSX);

+% show the geometry for checking

+CSXGeomPlot([Sim_Path '/' Sim_CSX]);

+% run the simulation

+RunOpenEMS( Sim_Path, Sim_CSX);

+

+% plot amplitude and phase of the reflection coefficient

+freq  = linspace(f_0-f_c, f_0+f_c, 200);

+port = calcPort(port, Sim_Path, freq);

+s11 = port.uf.ref ./ port.uf.inc;

+s11phase = unwrap(arg(s11));

+figure %("visible", "off"); % use this to plot only into files at the end of this script

+ax = plotyy( freq/1e6, 20*log10(abs(s11)), freq/1e6, s11phase);

+grid on

+title( ['reflection coefficient ', suffix, ' S_{11}']);

+xlabel( 'frequency f / MHz' );

+ylabel( ax(1), 'reflection coefficient |S_{11}|' );

+ylabel(ax(2), 'S_{11} phase (rad)');

+

+% define an azimuthal trace around the monopole

+phi = [0] * pi / 180;

+theta = [-180:10:180] * pi / 180;

+

+% calculate the delay, the fidelity and the farfield

+[delay, fidelity, nf2ff] = DelayFidelity(nf2ff, port, Sim_Path, sin(tilt), cos(tilt), theta, phi, f_0, f_c, 'Mode', 1);

+

+%plot the gain at (close to) f_0

+f_0_nearest_ind = nthargout(2, @min, abs(nf2ff.freq -f_0));

+%turn directivity into gain

+nf2ff.Dmax(f_0_nearest_ind) *= nf2ff.Prad(f_0_nearest_ind) / calcPort(port, Sim_Path, nf2ff.freq(f_0_nearest_ind)).P_inc; 

+figure %("visible", "off");

+polarFF(nf2ff, 'xaxis', 'theta', 'freq_index', f_0_nearest_ind, 'logscale', [-4, 4]);

+title(["gain ", suffix, " / dBi"]);

+

+

+% We trick polarFF into plotting the delay in mm because

+% the axes of the vanilla polar plot can not be scaled.

+plotvar = delay * C0 * 1000; 

+maxplot = 80;

+minplot = 30;

+nf2ff.Dmax(1) = 10^(max(plotvar)/10);

+nf2ff.E_norm{1} = 10.^(plotvar/20);

+figure %("visible", "off");

+polarFF(nf2ff, 'xaxis', 'theta', 'logscale', [minplot, maxplot]);

+title(["delay ", suffix, " / mm"]);

+

+% The same for the fidelity in percent.

+plotvar = fidelity * 100;

+maxplot = 100;

+minplot = 98;

+nf2ff.Dmax(1) = 10^(max(plotvar)/10);

+nf2ff.E_norm{1} = 10.^(plotvar/20);

+figure %("visible", "off");

+polarFF(nf2ff, 'xaxis', 'theta', 'logscale', [minplot, maxplot]);

+title(["fidelity ", suffix, " / %"]);

+

+% save the plots in order to compare them afer simulating the different channels

+print(1, ["s11_", suffix, ".png"]);

+print(2, ["farfield_", suffix, ".png"]);

+print(3, ["delay_mm_", suffix, ".png"]);

+print(4, ["fidelity_", suffix, ".png"]);

+return;
\ No newline at end of file
diff --git a/openEMS/matlab/Tutorials/Simple_Patch_Antenna.m b/openEMS/matlab/Tutorials/Simple_Patch_Antenna.m
index ee2a8f0..1e24c4c 100644
--- a/openEMS/matlab/Tutorials/Simple_Patch_Antenna.m
+++ b/openEMS/matlab/Tutorials/Simple_Patch_Antenna.m
@@ -1,20 +1,20 @@
-%
-% Tutorials / simple patch antenna
+%% Simple Patch Antenna Tutorial
 %
 % Describtion at:
-% http://openems.de/index.php/Tutorial:_Simple_Patch_Antenna
+% <http://openems.de/index.php/Tutorial:_Simple_Patch_Antenna>
 %
 % Tested with
 %  - Matlab 2013a / Octave 4.0
-%  - openEMS v0.0.33
+%  - openEMS v0.0.35
 %
-% (C) 2010-2015 Thorsten Liebig <thorsten.liebig@uni-due.de>
+% (C) 2010-2017 Thorsten Liebig <thorsten.liebig@uni-due.de>
+%%
 
 close all
 clear
 clc
 
-%% setup the simulation
+%% Setup the Simulation
 physical_constants;
 unit = 1e-3; % all length in mm
 
@@ -38,7 +38,7 @@ feed.R = 50;     %feed resistance
 % size of the simulation box
 SimBox = [200 200 150];
 
-%% setup FDTD parameter & excitation function
+%% Setup FDTD Parameter & Excitation Function
 f0 = 2e9; % center frequency
 fc = 1e9; % 20 dB corner frequency
 FDTD = InitFDTD( 'NrTs', 30000 );
@@ -46,7 +46,7 @@ FDTD = SetGaussExcite( FDTD, f0, fc );
 BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % boundary conditions
 FDTD = SetBoundaryCond( FDTD, BC );
 
-%% setup CSXCAD geometry & mesh
+%% Setup CSXCAD Geometry & Mesh
 CSX = InitCSX();
 
 %initialize the mesh with the "air-box" dimensions
@@ -54,13 +54,13 @@ mesh.x = [-SimBox(1)/2 SimBox(1)/2];
 mesh.y = [-SimBox(2)/2 SimBox(2)/2];
 mesh.z = [-SimBox(3)/3 SimBox(3)*2/3];
 
-%% create patch
+% Create Patch
 CSX = AddMetal( CSX, 'patch' ); % create a perfect electric conductor (PEC)
 start = [-patch.width/2 -patch.length/2 substrate.thickness];
 stop  = [ patch.width/2  patch.length/2 substrate.thickness];
 CSX = AddBox(CSX,'patch',10,start,stop); % add a box-primitive to the metal property 'patch'
 
-%% create substrate
+% Create Substrate
 CSX = AddMaterial( CSX, 'substrate' );
 CSX = SetMaterialProperty( CSX, 'substrate', 'Epsilon', substrate.epsR, 'Kappa', substrate.kappa );
 start = [-substrate.width/2 -substrate.length/2 0];
@@ -70,18 +70,18 @@ CSX = AddBox( CSX, 'substrate', 0, start, stop );
 % add extra cells to discretize the substrate thickness
 mesh.z = [linspace(0,substrate.thickness,substrate.cells+1) mesh.z];
 
-%% create ground (same size as substrate)
+% Create Ground same size as substrate
 CSX = AddMetal( CSX, 'gnd' ); % create a perfect electric conductor (PEC)
 start(3)=0;
 stop(3) =0;
 CSX = AddBox(CSX,'gnd',10,start,stop);
 
-%% apply the excitation & resist as a current source
+% Apply the Excitation & Resist as a Current Source
 start = [feed.pos 0 0];
 stop  = [feed.pos 0 substrate.thickness];
 [CSX port] = AddLumpedPort(CSX, 5 ,1 ,feed.R, start, stop, [0 0 1], true);
 
-%% finalize the mesh
+% Finalize the Mesh
 % detect all edges except of the patch
 mesh = DetectEdges(CSX, mesh,'ExcludeProperty','patch');
 % detect and set a special 2D metal edge mesh for the patch
@@ -90,35 +90,41 @@ mesh = DetectEdges(CSX, mesh,'SetProperty','patch','2D_Metal_Edge_Res', c0/(f0+f
 mesh = SmoothMesh(mesh, c0/(f0+fc)/unit/20);
 CSX = DefineRectGrid(CSX, unit, mesh);
 
-%% add a nf2ff calc box; size is 3 cells away from MUR boundary condition
+CSX = AddDump(CSX,'Hf', 'DumpType', 11, 'Frequency',[2.4e9]);
+CSX = AddBox(CSX,'Hf',10,[-substrate.width -substrate.length -10*substrate.thickness],[substrate.width +substrate.length 10*substrate.thickness]); %assign box
+
+% add a nf2ff calc box; size is 3 cells away from MUR boundary condition
 start = [mesh.x(4)     mesh.y(4)     mesh.z(4)];
 stop  = [mesh.x(end-3) mesh.y(end-3) mesh.z(end-3)];
 [CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', start, stop);
 
-%% prepare simulation folder
+%% Prepare and Run Simulation
 Sim_Path = 'tmp_Patch_Ant';
 Sim_CSX = 'patch_ant.xml';
 
+% create an empty working directory
 [status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
 [status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder
 
-%% write openEMS compatible xml-file
+% write openEMS compatible xml-file
 WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
 
-%% show the structure
+% show the structure
 CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
 
-%% run openEMS
+% run openEMS
 RunOpenEMS( Sim_Path, Sim_CSX);
 
-%% postprocessing & do the plots
+%% Postprocessing & Plots
 freq = linspace( max([1e9,f0-fc]), f0+fc, 501 );
 port = calcPort(port, Sim_Path, freq);
 
-Zin = port.uf.tot ./ port.if.tot;
-s11 = port.uf.ref ./ port.uf.inc;
+%% Smith chart port reflection
+plotRefl(port, 'threshold', -10)
+title( 'reflection coefficient' );
 
 % plot feed point impedance
+Zin = port.uf.tot ./ port.if.tot;
 figure
 plot( freq/1e6, real(Zin), 'k-', 'Linewidth', 2 );
 hold on
@@ -130,6 +136,7 @@ ylabel( 'impedance Z_{in} / Ohm' );
 legend( 'real', 'imag' );
 
 % plot reflection coefficient S11
+s11 = port.uf.ref ./ port.uf.inc;
 figure
 plot( freq/1e6, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 );
 grid on
@@ -139,7 +146,7 @@ ylabel( 'reflection coefficient |S_{11}|' );
 
 drawnow
 
-%% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% NFFF Plots
 %find resonance frequncy from s11
 f_res_ind = find(s11==min(s11));
 f_res = freq(f_res_ind);
@@ -166,7 +173,7 @@ plotFFdB(nf2ff,'xaxis','theta','param',[1 2])
 
 drawnow
 
-%%
+% Show 3D pattern
 disp( 'calculating 3D far field pattern and dumping to vtk (use Paraview to visualize)...' );
 thetaRange = (0:2:180);
 phiRange = (0:2:360) - 180;
diff --git a/openEMS/matlab/Tutorials/StripLine2MSL.m b/openEMS/matlab/Tutorials/StripLine2MSL.m
new file mode 100644
index 0000000..64eeaa1
--- /dev/null
+++ b/openEMS/matlab/Tutorials/StripLine2MSL.m
@@ -0,0 +1,133 @@
+%
+% Stripline to Microstrip Line Transition
+%
+% Describtion at:
+% <http://openems.de/index.php/Tutorial:_Stripline_to_MSL_Transition>
+%
+% Tested with
+%  - Octave 4.0
+%  - openEMS v0.0.35
+%
+% (C) 2017 Thorsten Liebig <thorsten.liebig@gmx.de>
+
+close all
+clear
+clc
+
+%% Setup the Simulation
+physical_constants;
+unit = 1e-6; % specify everything in um
+
+line_length     = 15000; % line length of strip line and microstrip line
+substrate_width = 6000;
+air_spacer      = 4000;  % air spacer above the substrate
+
+msl_width               = 500;
+msl_substrate_thickness = 254;
+
+strip_width               = 500;
+strip_substrate_thickness = 512;
+
+connect_via_rad =  500/2;
+connect_via_gap = 1250/2;
+
+substrate_epr    = 3.66;
+substrate_kappa  = 1e-3 * 2*pi*2.45e9 * EPS0*substrate_epr; % substrate losses
+
+f_max = 10e9;
+resolution = 250;
+edge_res   = 25;
+feed_shift = 2500;
+meas_shift = 5000;
+
+%% Setup FDTD Parameters & Excitation Function
+FDTD = InitFDTD();
+FDTD = SetGaussExcite( FDTD, f_max/2, f_max/2);
+BC   = {'PML_8' 'PML_8' 'PEC' 'PEC' 'PEC' 'MUR'};
+FDTD = SetBoundaryCond( FDTD, BC );
+
+%% Setup CSXCAD Geometry & Mesh
+CSX = InitCSX();
+edge_mesh  = [-1/3 2/3]*edge_res; % 1/3 - 2/3 rule for 2D metal edges
+
+mesh.x = SmoothMeshLines( [-connect_via_gap 0 connect_via_gap], 2*edge_res, 1.5 );
+mesh.x = SmoothMeshLines( [-line_length mesh.x line_length], resolution, 1.5);
+mesh.y = SmoothMeshLines( [0 msl_width/2+edge_mesh substrate_width/2], resolution/4 , 1.5);
+mesh.y = sort(unique([-mesh.y mesh.y]));
+mesh.z = SmoothMeshLines( [linspace(-strip_substrate_thickness,0,5) linspace(0,strip_substrate_thickness,5) linspace(strip_substrate_thickness,msl_substrate_thickness+strip_substrate_thickness,5) 2*strip_substrate_thickness+air_spacer] , resolution );
+CSX = DefineRectGrid( CSX, unit, mesh );
+
+% Create Substrate
+CSX = AddMaterial( CSX, 'RO4350B' );
+CSX = SetMaterialProperty( CSX, 'RO4350B', 'Epsilon', substrate_epr, 'Kappa', substrate_kappa );
+start = [mesh.x(1),   mesh.y(1),   -strip_substrate_thickness];
+stop  = [mesh.x(end), mesh.y(end), +strip_substrate_thickness+msl_substrate_thickness];
+CSX = AddBox( CSX, 'RO4350B', 0, start, stop );
+
+% Create a PEC called 'metal' and 'gnd'
+CSX = AddMetal( CSX, 'gnd' ); 
+CSX = AddMetal( CSX, 'metal' ); 
+
+% Create strip line port (incl. metal stip line)
+start = [-line_length -strip_width/2  0];
+stop  = [0            +strip_width/2  0];
+[CSX,port{1}] = AddStripLinePort( CSX, 100, 1, 'metal', start, stop, strip_substrate_thickness, 'x', [0 0 -1], 'ExcitePort', true, 'FeedShift', feed_shift, 'MeasPlaneShift', meas_shift );
+
+% Create MSL port on top
+start = [line_length  -strip_width/2 strip_substrate_thickness+msl_substrate_thickness];
+stop  = [0            +strip_width/2 strip_substrate_thickness];
+[CSX,port{2}] = AddMSLPort( CSX, 100, 2, 'metal', start, stop, 'x', [0 0 -1], 'MeasPlaneShift', meas_shift );
+
+% transitional via
+start = [0, 0, 0];
+stop  = [0, 0, strip_substrate_thickness+msl_substrate_thickness];
+CSX = AddCylinder(CSX, 'metal', 100, start, stop, connect_via_rad);
+
+% metal plane between strip line and MSL, including hole for transition
+p(1,1) = mesh.x(1);
+p(2,1) = mesh.y(1);
+p(1,2) = 0;
+p(2,2) = mesh.y(1);
+for a = linspace(-pi, pi, 21)
+  p(1,end+1) = connect_via_gap*sin(a);
+  p(2,end)   = connect_via_gap*cos(a);
+endfor
+p(1,end+1) = 0;
+p(2,end  ) = mesh.y(1);
+p(1,end+1) = mesh.x(end);
+p(2,end  ) = mesh.y(1);
+p(1,end+1) = mesh.x(end);
+p(2,end  ) = mesh.y(end);
+p(1,end+1) = mesh.x(1);
+p(2,end  ) = mesh.y(end);
+CSX = AddPolygon( CSX, 'gnd', 1, 'z', strip_substrate_thickness, p);
+
+%% Write/Show/Run the openEMS compatible xml-file
+Sim_Path = 'tmp';
+Sim_CSX = 'strip2msl.xml';
+
+[status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
+[status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder
+
+WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
+CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
+RunOpenEMS( Sim_Path, Sim_CSX );
+
+%% Post-Processing
+close all
+f = linspace( 0, f_max, 1601 );
+port = calcPort( port, Sim_Path, f, 'RefImpedance', 50);
+
+s11 = port{1}.uf.ref./ port{1}.uf.inc;
+s21 = port{2}.uf.ref./ port{1}.uf.inc;
+
+plot(f/1e9,20*log10(abs(s11)),'k-','LineWidth',2);
+hold on;
+grid on;
+plot(f/1e9,20*log10(abs(s21)),'r--','LineWidth',2);
+legend('S_{11}','S_{21}');
+ylabel('S-Parameter (dB)','FontSize',12);
+xlabel('frequency (GHz) \rightarrow','FontSize',12);
+ylim([-40 2]);
+
+
diff --git a/openEMS/matlab/h5readatt_octave.cc b/openEMS/matlab/h5readatt_octave.cc
index 8bd58d0..13e1765 100755
--- a/openEMS/matlab/h5readatt_octave.cc
+++ b/openEMS/matlab/h5readatt_octave.cc
@@ -5,13 +5,7 @@
 // this special treatment is necessary because Win32-Octave ships with a very old hdf5 version (1.6.10)
 void CloseH5Object(hid_t obj)
 {
-#if ((H5_VERS_MAJOR == 1) && (H5_VERS_MINOR == 6))
-	// try group close, than Dataset close
-	if (H5Gclose(obj)<0)
-		H5Dclose(obj);
-#else
 	H5Oclose(obj);
-#endif
 }
 
 DEFUN_DLD (h5readatt_octave, args, nargout, "h5readatt_octave(<File_Name>,<DataSet_Name>,<Attribute_Name>)")
@@ -30,7 +24,7 @@ DEFUN_DLD (h5readatt_octave, args, nargout, "h5readatt_octave(<File_Name>,<DataS
 	}
 
 	//suppress hdf5 error output
-	H5Eset_auto1(NULL, NULL);
+	H5Eset_auto2(H5E_DEFAULT, NULL, NULL);
 
 	hid_t file = H5Fopen( args(0).string_value().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT );
 	if (file==-1)
@@ -39,17 +33,7 @@ DEFUN_DLD (h5readatt_octave, args, nargout, "h5readatt_octave(<File_Name>,<DataS
 	  return retval;
 	}
 
-#if ((H5_VERS_MAJOR == 1) && (H5_VERS_MINOR == 6))
-	// this special treatment is necessary because Win32-Octave ships with a very old hdf5 version (1.6.10)
-	hid_t obj = -1;
-	//try opening the group
-	obj = H5Gopen(file, args(1).string_value().c_str());
-	//try opening the dataset if group failed
-	if (obj==-1)
-		obj = H5Dopen(file, args(1).string_value().c_str());
-#else
 	hid_t obj = H5Oopen(file, args(1).string_value().c_str(), H5P_DEFAULT);
-#endif
 
 	if (obj==-1)
 	{
@@ -59,7 +43,8 @@ DEFUN_DLD (h5readatt_octave, args, nargout, "h5readatt_octave(<File_Name>,<DataS
 	  return retval;
 	}
 
-	hid_t attr = H5Aopen_name(obj, args(2).string_value().c_str());
+	hid_t attr = H5Aopen_by_name(obj, ".", args(2).string_value().c_str(), H5P_DEFAULT, H5P_DEFAULT);
+
 	if (attr==-1)
 	{
 	  CloseH5Object(obj);
diff --git a/openEMS/matlab/plotRefl.m b/openEMS/matlab/plotRefl.m
new file mode 100644
index 0000000..5f88631
--- /dev/null
+++ b/openEMS/matlab/plotRefl.m
@@ -0,0 +1,145 @@
+function h = plotRefl(port, varargin) 

+

+%  h = plotRefl(port,varargin)

+%

+%  plot the reflection coefficient of a port into a Smith chart.

+%  left and right facing triangles mark the lower and upper cutoff

+%  frequency of the pass bands. An asterisk marks the frequnecy with 

+%  the lowest reflection.

+%

+% input:

+%   port:      port data structure. Call calcPort with an appropriate 

+%              frequency vector before calling this routine

+% 

+% output:      graphics handle for further modification of the plot.

+%

+% variable input:

+%   'precision':   - number of decimal places (floating point precision) 

+%                    for the frequency (always in MHz), default is 0

+%   'threshold':   - Threshold value (in dB) for the upper and lower

+%                    cutoff frequency, default is -3

+%   example:

+%       myport = calcPort(myport, Sim_Path, linspace(f_0-f_c, f_0+f_c, 200));

+%       plotRefl(myport);

+%

+% See also calcPort

+% 

+% openEMS matlab interface

+% -----------------------

+% author: Georg Michel

+

+%defaults

+precision = 0;

+threshold = -3;

+

+

+for n=1:2:numel(varargin)

+    if (strcmp(varargin{n},'precision')==1);

+        precision = varargin{n+1};

+    elseif (strcmp(varargin{n},'threshold')==1);

+        threshold = varargin{n+1};

+    else

+        warning('openEMS:polarFF',['unknown argument key: ''' varargin{n} '''']);

+    end

+end

+

+

+if ~isfield(port, 'uf')

+  error('Cannot plot the reflection coefficient. Please call calcPort first.');

+end

+

+s11 = port.uf.ref ./ port.uf.inc;

+ffmt = ['%.', num2str(precision), 'f'];

+

+

+

+figure; %new figure

+

+plot([-1, 1], [0, 0], 'k');

+

+

+axis ([-1.15, 1.15, -1.15, 1.15], "square");

+axis off;

+hold on

+

+ReZ = [.2; .5; 1; 2]; 

+ImZ = 1i * [1 2 5 2];

+Z = bsxfun(@plus, ReZ, linspace(-ImZ, ImZ, 256));

+Gamma = (Z-1)./(Z+1);

+plot(Gamma.', 'k');

+

+ReZ = [.5 .5 1 1 2 2 5 5 10 10];

+ImZ = 1i * [-.2; .2; -.5; .5; -1; 1; -2; 2; -5; 5];

+Z = bsxfun(@plus, linspace(0, ReZ, 256), ImZ);

+Gamma = (Z-1)./(Z+1);

+plot(Gamma.', 'k');

+

+

+angle = linspace (0, 2 * pi, 256); ReZ = [0 5 10];

+center = ReZ ./ (ReZ + 1);

+radius = 1 ./ (ReZ + 1);

+plot(bsxfun(@plus, bsxfun(@times, radius, cos(angle.')), center), bsxfun(@times, radius, sin(angle.')), 'k');

+

+

+% resistance

+ReZ = [0.2 0.5 1 2 5 10]; ImZ = zeros (1, length (ReZ));

+rho = (ReZ.^2 + ImZ.^2 - 1 + 2i * ImZ) ./ ((ReZ + 1).^2 + ImZ.^2);

+

+xoffset = [0.1 0.1 0.05 0.05 0.05 0.075]; 

+yoffset = -0.03;

+

+for idx = 1:length (ReZ)

+    text (real (rho(idx)) - xoffset(idx), ...

+        imag (rho(idx)) - yoffset, num2str (ReZ(idx)));

+end

+

+% reactance

+ReZ = [-0.06 -0.06 -0.06 -0.12 -0.5];

+ImZ = [0.2 0.5 1 2 5];

+

+

+rho = (ReZ.^2 + ImZ.^2 - 1 + 2i * ImZ) ./ ((ReZ + 1).^2 + ImZ.^2);

+

+for idx = 1:length (ImZ)

+    text (real (rho(idx)), imag (rho(idx)), [num2str(ImZ(idx)), "j"]);

+    text (real (rho(idx)), -imag (rho(idx)), [num2str(-ImZ(idx)), "j"]); end

+

+% zero

+rho = (-0.05.^2 + 0.^2 - 1) ./ ((-0.05 + 1).^2 + 0.^2);

+

+text (real (rho), imag (rho), '0');

+

+s11dB = 20*log10(abs(s11));

+

+upperind = s11dB(1:end-1) < threshold & s11dB(2:end) > threshold;

+lowerind = s11dB(1:end-1) > threshold & s11dB(2:end) < threshold;

+minind = nthargout(2, @min, s11dB);

+handle1 = plot(s11(lowerind),['<','b']);

+handle2 = plot(s11(upperind),['>','b']);

+handle3 = plot(s11(minind),['*', 'b']);

+llegend = num2str(port.f(lowerind)(1)/1e6, ffmt);

+ulegend = num2str(port.f(upperind)(1)/1e6, ffmt);

+

+if nnz(lowerind) > 1

+  for i= 2:nnz(lowerind)

+    llegend = strjoin({llegend, num2str(port.f(lowerind)(i)/1e6, ffmt)}, ', ');

+  end

+end

+

+if nnz(upperind) > 1

+  for i= 2:nnz(upperind)

+    ulegend = strjoin({ulegend, num2str(port.f(upperind)(i)/1e6, ffmt)}, ', ');

+  end

+end

+

+legend([handle1, handle2, handle3], {[llegend, " MHz"], ...

+                                     [ulegend, " MHz"], ...

+                                     [num2str(20*log10(abs(s11(minind))), "%4.0f"), ...

+                                    "dB @ ", num2str(port.f(minind)/1e6, ffmt), " MHz"]});

+h = plot(s11);

+

+if (nargout == 0)

+  clear h;

+end

+

+end
\ No newline at end of file
diff --git a/openEMS/matlab/setup.m b/openEMS/matlab/setup.m
index df56b25..e208ad4 100644
--- a/openEMS/matlab/setup.m
+++ b/openEMS/matlab/setup.m
@@ -5,7 +5,7 @@ function setup()
 %
 % openEMS matlab/octave interface
 % -----------------------
-% author: Thorsten Liebig (2011)
+% author: Thorsten Liebig (2011-2017)
 
 disp('setting up openEMS matlab/octave interface')
 
@@ -16,18 +16,22 @@ cd(dir);
 
 if isOctave()
     disp('compiling oct files')
-    fflush(stdout)
+    fflush(stdout);
     if isunix
-        [res, fn] = unix('find /usr/lib -name libhdf5.so');
-        if length(fn)>0
-            [hdf5lib_dir, hdf5lib_fn] = fileparts(fn);
+        [res, fn_so] = unix('find /usr/lib     -name libhdf5.so');
+        [res, fn_h]  = unix('find /usr/include -name hdf5.h');
+        if length(fn_so)>0 && length(fn_h)>0
+            [hdf5lib_dir, hdf5lib_fn] = fileparts(fn_so);
             disp(["HDF5 library path found at: " hdf5lib_dir])
-            mkoctfile(["-L" hdf5lib_dir ],"-lhdf5 -DH5_USE_16_API", "h5readatt_octave.cc")
+
+            [hdf5inc_dir, hdf5inc_fn] = fileparts(fn_h);
+            disp(["HDF5 include path found at: " hdf5inc_dir])
+            mkoctfile(["-L" hdf5lib_dir " -I" hdf5inc_dir],"-lhdf5", "h5readatt_octave.cc")
         else
-            mkoctfile -lhdf5 -DH5_USE_16_API h5readatt_octave.cc
+            mkoctfile -lhdf5 h5readatt_octave.cc
         end
     else
-        mkoctfile -lhdf5 -DH5_USE_16_API h5readatt_octave.cc
+        mkoctfile -lhdf5 h5readatt_octave.cc
     end
 else
     disp('Matlab does not need this function. It is Octave only.')