diff --git a/matlab/AnalyzeNFFF2.m b/matlab/AnalyzeNFFF2.m
new file mode 100644
index 0000000..5a09e1f
--- /dev/null
+++ b/matlab/AnalyzeNFFF2.m
@@ -0,0 +1,208 @@
+function [E_theta,E_phi,Prad,Dmax] = AnalyzeNFFF2( Sim_Path, filenames_E, filenames_H, f, theta, phi, r )
+% [E_theta,E_phi,Prad,Dmax] = AnalyzeNFFF2( Sim_Path, filenames_E, filenames_H, f, theta, phi, r )
+%
+% calculates the farfield via a near field to far field transformation
+%
+% input:
+%   Sim_Path:    simulation directory
+%   filenames_E: cell array of filenames for the time domain fields on the NFFF contour (6 E-planes; hdf5)
+%   filenames_H: cell array of filenames for the time domain fields on the NFFF contour (6 H-planes; hdf5)
+%   f:           frequency (Hz) for far field calculation
+%   theta:       (degrees) vector of discrete theta values to calculate the far field for
+%   phi:         (degrees) vector of discrete phi values to calculate the far field for
+%   r:           (optional) Radius (m) at which the E-fields are calculated (default: 1 m)
+%
+% output:
+%   E_theta:     E_theta(theta,phi); theta component of the electric field strength at radius r
+%   E_phi:       E_phi(theta,phi);     phi component of the electric field strength at radius r
+%   Prad:        time averaged radiated power
+%   Dmax:        maximum directivity
+%
+% example:
+%   see examples/NF2FF/infDipol.m
+%
+% (C) 2010 Sebastian Held <sebastian.held@gmx.de>
+
+% check arguments
+error( nargchk(7,7,nargin) );
+if ~isscalar(f)
+    error 'Currently only one frequency is supported. Call this function multiple times.'
+end
+
+
+
+% read time domain field data and transform into frequency domain
+for n=1:numel(filenames_E)
+    [Ef{n}, E_mesh{n}] = ReadHDF5Dump( [Sim_Path '/' filenames_E{n}], 'Frequency', f );
+    [Hf{n}, H_mesh{n}] = ReadHDF5Dump( [Sim_Path '/' filenames_H{n}], 'Frequency', f );
+    
+    % reshape mesh into row vector
+    mesh{n}.x = reshape( E_mesh{n}.lines{1}, 1, [] );
+    mesh{n}.y = reshape( E_mesh{n}.lines{2}, 1, [] );
+    mesh{n}.z = reshape( E_mesh{n}.lines{3}, 1, [] );
+end
+
+% create a normal vector for every plane
+% FIXME!!! this is dependent upon the order of filenames_*
+n = {};
+for a=1:6
+    temp = [(a<=2), ((a>=3)&&(a<=4)), (a>=5)];
+    n{a} = temp - 2*mod(a,2)*temp;
+end
+
+
+physical_constants
+
+k = 2*pi*f/c0;
+center = [0 0 0];
+Umax = 0;
+
+phi_idx = 0;
+for phi_deg_aufpunkt = phi
+    phi_rad_aufpunkt = phi_deg_aufpunkt/180*pi; % radiant
+    phi_idx = phi_idx + 1;
+    
+    theta_idx = 0;
+    for theta_deg_aufpunkt = theta
+        theta_rad_aufpunkt = theta_deg_aufpunkt/180*pi; % radiant
+        theta_idx = theta_idx + 1;
+
+        N_theta = 0;
+        N_phi = 0;
+        L_theta = 0;
+        L_phi = 0;
+        for a=1:6
+            [N_theta_,N_phi_,L_theta_,L_phi_] = process_plane( k, n{a}, center, mesh{a}, Ef{a}.FD.values{1}, Hf{a}.FD.values{1}, theta_rad_aufpunkt, phi_rad_aufpunkt );
+            N_theta = N_theta + N_theta_; N_phi = N_phi + N_phi_;
+            L_theta = L_theta + L_theta_; L_phi = L_phi + L_phi_;
+        end
+
+        % E-fields
+        erg_E_theta = -1i*k*exp(-1i*k*r) / (4*pi*r)*(L_phi+Z0*N_theta);
+        erg_E_phi   =  1i*k*exp(-1i*k*r) / (4*pi*r)*(L_theta-Z0*N_phi);
+
+        % output
+        E_theta(theta_idx,phi_idx) = erg_E_theta;
+        E_phi(theta_idx,phi_idx)   = erg_E_phi;
+
+        % directivity
+        U = r^2/(2*Z0) * sum(abs([erg_E_theta erg_E_phi]).^2);
+        Umax = max( [Umax U] );
+    end
+end
+
+% power
+Prad = 0;
+for a=1:6
+    [~,~,~,~,P] = process_plane( k, n{a}, center, mesh{a}, Ef{a}.FD.values{1}, Hf{a}.FD.values{1}, theta_rad_aufpunkt, phi_rad_aufpunkt );
+    Prad = Prad + P;
+end
+
+% directivity
+Dmax = 4*pi*Umax / Prad;
+
+
+% integrate over one plane
+function [N_theta,N_phi,L_theta,L_phi,Prad] = process_plane( k, n, center, mesh, E_field, H_field, theta_rad_aufpunkt, phi_rad_aufpunkt )
+% [N_theta,N_phi,L_theta,L_phi,Prad] = process_plane( k, n, center, mesh, E_field, H_field, theta_rad_aufpunkt, phi_rad_aufpunkt )
+%
+% k: wave number
+% n: normal vector of the plane
+% center: correction coordinates for the center of the antenna
+% mesh: mesh info
+% E_field: E field array ?x?x?x3
+% H_field: H field array ?x?x?x3
+
+% speed up
+sin__theta_rad_aufpunkt = sin(theta_rad_aufpunkt);
+cos__theta_rad_aufpunkt = cos(theta_rad_aufpunkt);
+sin__phi_rad_aufpunkt   = sin(phi_rad_aufpunkt);
+cos__phi_rad_aufpunkt   = cos(phi_rad_aufpunkt);
+
+if abs(n(1)) == 1
+    % x-plane
+    x = mesh.x(1);
+    [y z] = ndgrid( mesh.y, mesh.z );
+    coord1 = mesh.y.';
+    coord2 = mesh.z.';
+    Ex = squeeze( E_field(1,:,:,1) );
+    Ey = squeeze( E_field(1,:,:,2) );
+    Ez = squeeze( E_field(1,:,:,3) );
+    Hx = squeeze( H_field(1,:,:,1) );
+    Hy = squeeze( H_field(1,:,:,2) );
+    Hz = squeeze( H_field(1,:,:,3) );
+elseif abs(n(2)) == 1
+    % y-plane
+    y = mesh.y(1);
+    [x z] = ndgrid( mesh.x, mesh.z );
+    coord1 = mesh.x.';
+    coord2 = mesh.z.';
+    Ex = squeeze( E_field(:,1,:,1) );
+    Ey = squeeze( E_field(:,1,:,2) );
+    Ez = squeeze( E_field(:,1,:,3) );
+    Hx = squeeze( H_field(:,1,:,1) );
+    Hy = squeeze( H_field(:,1,:,2) );
+    Hz = squeeze( H_field(:,1,:,3) );
+elseif abs(n(3)) == 1
+    % z-plane
+    z = mesh.z(1);
+    [x y] = ndgrid( mesh.x, mesh.y );
+    coord1 = mesh.x.';
+    coord2 = mesh.y.';
+    Ex = squeeze( E_field(:,:,1,1) );
+    Ey = squeeze( E_field(:,:,1,2) );
+    Ez = squeeze( E_field(:,:,1,3) );
+    Hx = squeeze( H_field(:,:,1,1) );
+    Hy = squeeze( H_field(:,:,1,2) );
+    Hz = squeeze( H_field(:,:,1,3) );
+end    
+
+Jx = n(2) .* Hz - n(3) .* Hy;
+Jy = n(3) .* Hx - n(1) .* Hz;
+Jz = n(1) .* Hy - n(2) .* Hx;
+Mx = -n(2) .* Ez + n(3) .* Ey;
+My = -n(3) .* Ex + n(1) .* Ez;
+Mz = -n(1) .* Ey + n(2) .* Ex;
+r_cos_psi = x*sin__theta_rad_aufpunkt*cos__phi_rad_aufpunkt + y*sin__theta_rad_aufpunkt*sin__phi_rad_aufpunkt + z*cos__theta_rad_aufpunkt;
+e_fkt     = exp( +1i*k*r_cos_psi );
+N_theta = dbltrapz( ( Jx*cos__theta_rad_aufpunkt*cos__phi_rad_aufpunkt + Jy*cos__theta_rad_aufpunkt*sin__phi_rad_aufpunkt - Jz*sin__theta_rad_aufpunkt) .* e_fkt, coord1, coord2 );
+N_phi   = dbltrapz( (-Jx*sin__phi_rad_aufpunkt + Jy*cos__phi_rad_aufpunkt) .* e_fkt, coord1, coord2 );
+L_theta = dbltrapz( ( Mx*cos__theta_rad_aufpunkt*cos__phi_rad_aufpunkt + My*cos__theta_rad_aufpunkt*sin__phi_rad_aufpunkt - Mz*sin__theta_rad_aufpunkt) .* e_fkt, coord1, coord2 );
+L_phi   = dbltrapz( (-Mx*sin__phi_rad_aufpunkt + My*cos__phi_rad_aufpunkt) .* e_fkt, coord1, coord2 );
+
+if nargout > 4
+    % Prad requested
+    
+    % this is crap! recode it!
+    EH = zeros(size(Ex));
+    for i1 = 1:numel(coord1)
+        for i2 = 1:numel(coord2)
+            E = [Ex(i1,i2) Ey(i1,i2) Ez(i1,i2)];
+            H = [Hx(i1,i2) Hy(i1,i2) Hz(i1,i2)];
+            EH(i1,i2) = real( dot(cross(E,conj(H)),n) );
+        end
+    end
+    Prad = 0.5 * dbltrapz( EH, coord1, coord2 );
+end
+
+
+
+
+function Q = dbltrapz(matrix,a,b)
+%DBLTRAPZ  Trapezoidal numerical integration in two dimensions.
+%   Z = DBLTRAPZ(MATRIX,A,B) computes an approximation of the double integral
+%   of MATRIX via the trapezoidal method (with respect to A and B).  A and B must be
+%   column vectors of the same length.
+%   index like this: MATRIX(A,B)
+
+if nargin < 3, error('MATLAB:dblquad:NotEnoughInputs',...
+                     'Requires at least three inputs.'); end
+if size(a,2) ~= 1, error('column vectors required'); end
+if size(b,2) ~= 1, error('column vectors required'); end
+
+temp = zeros(size(b));
+for i = 1:length(b) 
+    temp(i) = trapz( a, matrix(:,i) );
+end
+
+Q = trapz( b, temp );
diff --git a/matlab/examples/NF2FF/infDipol.m b/matlab/examples/NF2FF/infDipol.m
new file mode 100644
index 0000000..d72d316
--- /dev/null
+++ b/matlab/examples/NF2FF/infDipol.m
@@ -0,0 +1,212 @@
+function infDipol
+%
+% infinitesimal dipole example
+%
+
+close all
+clear
+clc
+drawnow
+
+
+postprocessing_only = 0;
+
+
+
+
+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';
+
+g.f_max  = 1e9;
+g.lambda = c0/g.f_max;
+
+% setup geometry values
+g.dipole_length = g.lambda/50 /g.drawingunit;
+
+
+
+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 );
+
+% 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;
+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 = AddBox( AddDump(CSX,'Et_xn','DumpType',0,'DumpMode',2,'FileType',1), 'Et_xn', 0, s1, [s1(1) s2(2) s2(3)] );
+CSX = AddBox( AddDump(CSX,'Et_xp','DumpType',0,'DumpMode',2,'FileType',1), 'Et_xp', 0, [s2(1) s1(2) s1(3)], s2 );
+CSX = AddBox( AddDump(CSX,'Et_yn','DumpType',0,'DumpMode',2,'FileType',1), 'Et_yn', 0, s1, [s2(1) s1(2) s2(3)] );
+CSX = AddBox( AddDump(CSX,'Et_yp','DumpType',0,'DumpMode',2,'FileType',1), 'Et_yp', 0, [s1(1) s2(2) s1(3)], s2 );
+CSX = AddBox( AddDump(CSX,'Et_zn','DumpType',0,'DumpMode',2,'FileType',1), 'Et_zn', 0, s1, [s2(1) s2(2) s1(3)] );
+CSX = AddBox( AddDump(CSX,'Et_zp','DumpType',0,'DumpMode',2,'FileType',1), 'Et_zp', 0, [s1(1) s1(2) s2(3)], s2 );
+CSX = AddBox( AddDump(CSX,'Ht_xn','DumpType',1,'DumpMode',2,'FileType',1), 'Ht_xn', 0, s1, [s1(1) s2(2) s2(3)] );
+CSX = AddBox( AddDump(CSX,'Ht_xp','DumpType',1,'DumpMode',2,'FileType',1), 'Ht_xp', 0, [s2(1) s1(2) s1(3)], s2 );
+CSX = AddBox( AddDump(CSX,'Ht_yn','DumpType',1,'DumpMode',2,'FileType',1), 'Ht_yn', 0, s1, [s2(1) s1(2) s2(3)] );
+CSX = AddBox( AddDump(CSX,'Ht_yp','DumpType',1,'DumpMode',2,'FileType',1), 'Ht_yp', 0, [s1(1) s2(2) s1(3)], s2 );
+CSX = AddBox( AddDump(CSX,'Ht_zn','DumpType',1,'DumpMode',2,'FileType',1), 'Ht_zn', 0, s1, [s2(1) s2(2) s1(3)] );
+CSX = AddBox( AddDump(CSX,'Ht_zp','DumpType',1,'DumpMode',2,'FileType',1), 'Ht_zp', 0, [s1(1) s1(2) s2(3)], s2 );
+
+
+
+
+
+
+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
+
+disp( ' ' );
+disp( ' ********************************************************** ' );
+disp( ' ' );
+
+% NFFF contour
+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'};
+
+% 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] = AnalyzeNFFF2( g.Sim_Path, filenames_E, filenames_H, g.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)
+    E_phi0_far(n) = norm( [E_far_theta(n,1) E_far_phi(n,1)] );
+end
+
+% display polar plot
+figure
+polar( thetaRange/180*pi, E_phi0_far );
+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)
+    E_phi90_far(n) = norm([E_far_theta(n,2) E_far_phi(n,2)]);
+end
+
+% display polar plot
+figure
+polar( thetaRange/180*pi, E_phi90_far );
+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] = AnalyzeNFFF2( g.Sim_Path, filenames_E, filenames_H, g.f_max, 90, phiRange, r );
+
+E_theta90_far = zeros(1,numel(phiRange));
+for n=1:numel(phiRange)
+    E_theta90_far(n) = norm([E_far_theta(1,n) E_far_phi(1,n)]);
+end
+
+% display polar plot
+figure
+polar( phiRange/180*pi, E_theta90_far );
+ylabel( 'phi / deg' );
+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] = AnalyzeNFFF2( g.Sim_Path, filenames_E, filenames_H, g.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);
+x = E_far_normalized .* sin(theta) .* cos(phi);
+y = E_far_normalized .* sin(theta) .* sin(phi);
+z = E_far_normalized .* cos(theta);
+figure
+surf( x,y,z, E_far_normalized );
+axis equal
+xlabel( 'x' );
+ylabel( 'y' );
+zlabel( 'z' );