diff --git a/tunnel_breakthrough_error.m b/tunnel_breakthrough_error.m
new file mode 100644
index 0000000..352200c
--- /dev/null
+++ b/tunnel_breakthrough_error.m
@@ -0,0 +1,19 @@
+function sigma=tunnel_breakthrough_error(n,s,sigma_angle)
+% The standard displacement error (sigma) of a geodetic station in tunneling in a
+% direction normal to the axis of the tunnel (lateral error). This error describes tunnel
+% breakthrough error (offset), and can be used in tunnel geodetic network pre-analysis.
+% Vectorized.
+% Version: 13 Jun 22
+% Input:   n - the number of the traverse legs
+%          s - the length of the open tunnel traverse [meter]
+%          sigma_angle- the standard deviation (uncertainty) of angle measurements [radians]
+% Output:  sigma - standard error (sigma) of a geodetic station in the tunnel [meter]
+% Reference: Chrzanowski, A., 1981. Optimization of the breakthrough accuracy in tunnelling surveys. Canadian Surveyor 35 (1), 5 – 16.
+% Keywords: tunneling, surveying, standard accuracy, breakthrough error
+%
+% Copyright (c) 2022, Matej Varga
+% All rights reserved.
+% Email: mvarga1989@gmail.com
+sigma = n.*(n+1).*(2.*n+1)./6.*s.^2.*sigma_angle.^2;
+
+end
diff --git a/tunnel_breakthrough_error_demo.m b/tunnel_breakthrough_error_demo.m
new file mode 100644
index 0000000..d8baa55
--- /dev/null
+++ b/tunnel_breakthrough_error_demo.m
@@ -0,0 +1,25 @@
+clear;clc;close all
+
+n= [0:20]; % number of stations
+s= 50; % leg distances between geodetic marks  in a tunnel_breakthrough_error
+sigma_angle = 10; % the standard deviation for the angles measurement in arc seconds
+
+
+%
+sigma_angle = sigma_angle/3600*pi/180; % convert angle units to radians
+
+sigma_breakthrough=tunnel_breakthrough_error(n,s,sigma_angle);
+
+figure
+plot(n,sigma_breakthrough*1000);
+grid on;
+xlabel('Number of stations');
+ylabel('Lateral displacement error [mm]');
+
+figure
+plot(n*s,sigma_breakthrough*1000);
+grid on;
+xlabel('Traverse length [meter]');
+ylabel('Lateral displacement error [mm]');
+
+