diff --git a/matlab/plotRefl.m b/matlab/plotRefl.m
new file mode 100644
index 0000000..5f88631
--- /dev/null
+++ b/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