diff --git a/matlab/examples/microstrip/directional_coupler.m b/matlab/examples/microstrip/directional_coupler.m
new file mode 100644
index 0000000..15e0fb0
--- /dev/null
+++ b/matlab/examples/microstrip/directional_coupler.m
@@ -0,0 +1,235 @@
+function directional_coupler
+%
+% EXAMPLE / microstrip / directional_coupler
+%
+% Stacked directional coupler in microstrip technology.
+%
+% This example demonstrates:
+%  - simple microstrip geometry
+%  - S-parameter calculation using the ypar-method
+%  - display of coupler parameters
+%  - display of S11 (smith chart)
+% 
+%
+% Tested with
+%  - Matlab 2010b
+%  - Octave 3.3.54
+%  - openEMS v0.0.17
+%
+% (C) 2010 Sebastian Held <sebastian.held@gmx.de>
+
+clear
+close all
+clc
+
+% sim settings
+showStructure = 1;
+runSimulation = 1;
+
+for n=1:4
+    if n > 1, showStructure = 0; end
+    ports{n} = sim( n, showStructure, runSimulation );
+end
+postprocess( ports );
+
+
+
+
+function ports = sim( simnr, showStructure, runSimulation )
+physical_constants
+
+% setup the simulation
+drawingunit = 1e-6; % specify everything in um
+Sim_Path = ['tmp' int2str(simnr)];
+Sim_CSX  = 'tmp.xml';
+f_max    = 100e6;
+lambda   = c0/f_max;
+
+% specify the coupler
+pcb1.w   = 147000;
+pcb1.h   = 54500;
+pcb1.t   = 1524;
+pcb1.epr = 3;
+msl1.w   = 135000;
+msl1.h   = 2800;
+pcb2.w   = 107000;
+pcb2.h   = 14000;
+pcb2.t   = 1524;
+pcb2.epr = 3;
+msl2.w   = 95000;
+msl2.h   = 4000;
+
+
+CSX = InitCSX();
+
+% create the mesh
+mesh.x = [-pcb1.w/2 pcb1.w/2 -pcb2.w/2 pcb2.w/2 -msl1.w/2 msl1.w/2 -msl2.w/2 msl2.w/2];
+mesh.x = [mesh.x linspace(-msl2.w/2,-msl2.w/2+msl2.h, 5) linspace(msl2.w/2,msl2.w/2-msl2.h, 5)];
+mesh.y = [-pcb1.h/2 pcb1.h/2 -pcb2.h/2 pcb2.h/2 -msl1.h/2 msl1.h/2 -msl2.h/2 msl2.h/2];
+mesh.z = [linspace(0,pcb1.t,5) linspace(pcb1.t,pcb1.t+pcb2.t,5)];
+mesh.z = [mesh.z mesh.z(end)+10*(mesh.z(end)-mesh.z(1))]; % add space above pcb
+res = lambda/sqrt(max([pcb1.epr,pcb2.epr])) / 20 / drawingunit;
+mesh.x = SmoothMeshLines2(mesh.x,res);
+mesh.y = SmoothMeshLines2(mesh.y,res);
+mesh.z = SmoothMeshLines2(mesh.z,res);
+mesh = AddPML( mesh, [8 8 8 8 8 8] ); % add space for PML
+CSX = DefineRectGrid( CSX, drawingunit, mesh );
+
+%% create the structure
+
+% microstrip
+CSX = AddMetal( CSX, 'PEC' );
+start = [-msl1.w/2, -msl1.h/2, pcb1.t];
+stop  = [ msl1.w/2,  msl1.h/2, pcb1.t];
+priority = 100; % the geometric priority is set to 100
+CSX = AddBox( CSX, 'PEC', priority, start, stop );
+
+% ground plane
+CSX = AddMetal( CSX, 'PEC_ground' );
+start = [-pcb1.w/2, -pcb1.h/2, 0];
+stop  = [ pcb1.w/2,  pcb1.h/2, 0];
+CSX = AddBox( CSX, 'PEC_ground', priority, start, stop );
+
+% substrate 1
+start = [-pcb1.w/2, -pcb1.h/2, 0];
+stop  = [ pcb1.w/2,  pcb1.h/2, pcb1.t];
+priority = 10;
+CSX = AddMaterial( CSX, 'substrate1' );
+CSX = SetMaterialProperty( CSX, 'substrate1', 'Epsilon', pcb1.epr );
+CSX = AddBox( CSX, 'substrate1', priority, start, stop );
+
+% substrate 2
+start = [-pcb2.w/2, -pcb2.h/2, pcb1.t];
+stop  = [ pcb2.w/2,  pcb2.h/2, pcb1.t+pcb2.t];
+priority = 10;
+CSX = AddMaterial( CSX, 'substrate2' );
+CSX = SetMaterialProperty( CSX, 'substrate2', 'Epsilon', pcb2.epr );
+CSX = AddBox( CSX, 'substrate2', priority, start, stop );
+
+% stripline
+start = [-msl2.w/2, -msl2.h/2, pcb1.t+pcb2.t];
+stop  = [ msl2.w/2,  msl2.h/2, pcb1.t+pcb2.t];
+priority = 100;
+CSX = AddBox( CSX, 'PEC', priority, start, stop );
+
+% connections
+start = [-msl2.w/2,        -msl2.h/2, pcb1.t+pcb2.t];
+stop  = [-msl2.w/2+msl2.h, -pcb2.h/2, pcb1.t+pcb2.t];
+priority = 100;
+CSX = AddBox( CSX, 'PEC', priority, start, stop );
+start = [ msl2.w/2,        -msl2.h/2, pcb1.t+pcb2.t];
+stop  = [ msl2.w/2-msl2.h, -pcb2.h/2, pcb1.t+pcb2.t];
+priority = 100;
+CSX = AddBox( CSX, 'PEC', priority, start, stop );
+
+%% ports
+% this project needs 4 simulations
+for n=1:4
+    portexcite{n} = [];
+end
+portexcite{simnr} = 'excite';
+
+% port 1: input port
+start = [-msl1.w/2, 0, pcb1.t];
+stop  = [-msl1.w/2, 0, 0];
+[CSX ports{1}] = AddCurvePort( CSX, 1, 50, start, stop, portexcite{1} );
+% port 2: output port
+start = [msl1.w/2, 0, pcb1.t];
+stop  = [msl1.w/2, 0, 0];
+[CSX ports{2}] = AddCurvePort( CSX, 2, 50, start, stop, portexcite{2} );
+% port 3: coupled port
+start = [-msl2.w/2+msl2.h/2, -pcb2.h/2, pcb1.t+pcb2.t];
+stop  = [-msl2.w/2+msl2.h/2, -pcb2.h/2, 0];
+[CSX ports{3}] = AddCurvePort( CSX, 3, 50, start, stop, portexcite{3} );
+% port 4: isolated port
+start = [msl2.w/2-msl2.h/2, -pcb2.h/2, pcb1.t+pcb2.t];
+stop  = [msl2.w/2-msl2.h/2, -pcb2.h/2, 0];
+[CSX ports{4}] = AddCurvePort( CSX, 4, 50, start, stop, portexcite{4} );
+
+%% setup FDTD parameters & excitation function
+max_timesteps = 50000;
+min_decrement = 1e-6;
+FDTD = InitFDTD( max_timesteps, min_decrement );
+FDTD = SetGaussExcite( FDTD, 0, f_max );
+BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'};
+BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % faster
+FDTD = SetBoundaryCond( FDTD, BC );
+
+%% Write openEMS compatible xml-file
+if runSimulation
+    [~,~,~] = rmdir(Sim_Path,'s');
+end
+[~,~,~] = mkdir(Sim_Path);
+WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX);
+
+if showStructure
+    CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
+end
+
+%% run openEMS
+openEMS_opts = '';
+openEMS_opts = [openEMS_opts ' --engine=fastest'];
+% openEMS_opts = [openEMS_opts ' --debug-material'];
+% openEMS_opts = [openEMS_opts ' --debug-boxes'];
+if runSimulation
+    RunOpenEMS( Sim_Path, Sim_CSX, openEMS_opts );
+end
+
+
+
+
+function postprocess( ports )
+f = linspace( 0, 100e6, 201 );
+Y = calc_ypar( f, ports{1}, 'tmp' );
+R = 50;
+S = y2s(Y,R);
+
+% insertion loss
+IL_dB = -20 * log10(abs(squeeze(S(2,1,:))));
+
+% coupling factor
+CF_dB = -20 * log10(abs(squeeze(S(3,1,:))));
+
+% isolation
+I_dB  = -20 * log10(abs(squeeze(S(4,1,:))));
+
+% directivity
+D_dB  = -20 * log10(abs(squeeze(S(4,1,:) ./ S(3,1,:))));
+
+figure
+plot( f, [IL_dB CF_dB I_dB D_dB] );
+legend( {'insertion loss','coupling factor','isolation','directivity'} );
+title( ['performance of the coupler for a termination resistance of R=' num2str(R)] );
+grid on
+
+smithchart
+S11 = squeeze(S(1,1,:));
+plot( real(S11), imag(S11) );
+legend( 'S_{11}' );
+title( ['performance of the coupler for a termination resistance of R=' num2str(R)] );
+axis( [-1 1 -1 1] );
+
+
+
+function smithchart
+% smith chart
+figure
+if exist( 'smith', 'file' )
+    % smith chart
+    % www.ece.rutgers.edu/~orfanidi/ewa
+    % or cmt toolbox from git.ate.uni-duisburg.de
+    smith
+else
+    % poor man smith chart
+    color = [.6 .6 .6];
+    h = plot( sin(0:0.01:2*pi), cos(0:0.01:2*pi), 'Color', color );
+    hg = hggroup;
+    set( h,'Parent',hg );
+    hold on
+    plot( hg, 0.25+0.75*sin(0:0.01:2*pi), 0.75*cos(0:0.01:2*pi), 'Color', color );
+    plot( hg, 0.5+0.5*sin(0:0.01:2*pi), 0.5*cos(0:0.01:2*pi), 'Color', color );
+    plot( hg, 0.75+0.25*sin(0:0.01:2*pi), 0.25*cos(0:0.01:2*pi), 'Color', color );
+    plot( hg, [-1 1], [0 0], 'Color', color );
+    axis equal
+    axis off
+end