diff --git a/Matlab/axial.m b/Matlab/axial.m
index 1be3bdb..41af23a 100644
--- a/Matlab/axial.m
+++ b/Matlab/axial.m
@@ -12,7 +12,7 @@ load_constants
 % Intensity profile plots
 rho = linspace(-r_max, r_max, 500);
 figure
-I = gauss(rho, r_max, w0);
+I = gauss(rho, w0, r_max);
 I0 = max(I);
 plot(rho, I/I0, 'k')
 grid
@@ -20,12 +20,11 @@ xlabel('r, м')
 ylabel('I(r)')
 
 % Integration
-G0 = gauss_peak(r_max, w0);
-Qres_g = @(z) 2 * pi * G0 * integral2(@(beta, r) r .* gauss(r, r_max, w0) .* ...
+Qres_g = @(z) 1 / (pi * w0 ^ 2) * integral2(@(beta, r) r .* gauss(r, w0, r_max) .* ...
     iscomplex(qg_z_factor(r, z, n1, n2, Rsp, f)), 0, 2*pi, 0, r_max, ...
     'Method', 'iterated', 'AbsTol', 1e-12, 'RelTol', 1e-6);
 
-Qres_s = @(z) 2 * pi * G0 * integral2(@(beta, r) r .* gauss(r, r_max, w0) .* ...
+Qres_s = @(z) 1 / (pi * w0 ^ 2) * integral2(@(beta, r) r .* gauss(r, w0, r_max) .* ...
     iscomplex(qs_z_factor(r, z, n1, n2, Rsp, f)), 0, 2*pi, 0, r_max, ...
     'Method', 'iterated', 'AbsTol', 1e-12, 'RelTol', 1e-6);
 
diff --git a/Matlab/bessel.m b/Matlab/bessel.m
index 413aaea..281e296 100644
--- a/Matlab/bessel.m
+++ b/Matlab/bessel.m
@@ -1,4 +1,6 @@
 % Bessel beam
-function b = bessel(r, r_max, w0, P)
-    b = besselj(0, 2.405/w0 * r).^2;
+function b = bessel(r, w0, r_max)
+    ring_radius = 2.405; % radisu of the first ring of the besselj_0
+    A = 2 * ring_radius / (w0 * r_max * exp(0.5)) * 2 * pi * integral(@(r) r.*besselj(0, ring_radius/w0 * r).^2, 0, r_max);
+    b = A * besselj(0, 2.405/w0 * r).^2;
 end
\ No newline at end of file
diff --git a/Matlab/bessel_peak.m b/Matlab/bessel_peak.m
deleted file mode 100644
index 66c1834..0000000
--- a/Matlab/bessel_peak.m
+++ /dev/null
@@ -1,3 +0,0 @@
-function peak = bessel_peak(r_max, w0, P)
-    peak = P * 4.81 / (w0 * r_max * exp(0.5)) * 2 * pi * integral(@(r) r.*besselj(0, 2.405/w0 * r).^2, 0, r_max) / P0;
-end
\ No newline at end of file
diff --git a/Matlab/gauss.m b/Matlab/gauss.m
index 4e15277..02b5ae7 100644
--- a/Matlab/gauss.m
+++ b/Matlab/gauss.m
@@ -1,4 +1,5 @@
 % Gaussian TEM00 beam
-function g = gauss(r, r_max, w0)
-    g = exp(-2 * r.^2 / w0^2);
+function g = gauss(r, w0, r_max)
+    A = (1 - exp(-2*r_max.^2 / w0^2)); % the fraction of power that falls on the pupil of the micro lens
+    g = 2 * A * exp(-2 * r.^2 / w0^2);
 end
\ No newline at end of file
diff --git a/Matlab/gauss_peak.m b/Matlab/gauss_peak.m
deleted file mode 100644
index 4216cd5..0000000
--- a/Matlab/gauss_peak.m
+++ /dev/null
@@ -1,4 +0,0 @@
-function peak = gauss_peak(r_max, w0)
-    A = (1 - exp(-2*r_max.^2 / w0^2));
-    peak = 2*A / (pi * w0^2);
-end
\ No newline at end of file
diff --git a/Matlab/transverse.m b/Matlab/transverse.m
index eab9a5c..7b27416 100644
--- a/Matlab/transverse.m
+++ b/Matlab/transverse.m
@@ -12,7 +12,7 @@ load_constants
 
 % Intensity profile graphics
 rho = linspace(-r_max, r_max, 500);
-I = gauss(rho, r_max, w0);
+I = gauss(rho, w0, r_max);
 I0 = max(I);
 figure
 plot(rho, I/I0, 'k')
@@ -21,12 +21,11 @@ xlabel('r, m')
 ylabel('I(r)')
 
 % Integration
-G0 = gauss_peak(r_max, w0);
-Qres_g = @(y) G0 * integral2(@(beta, r) r .* gauss(r, r_max, w0) .* ...
+Qres_g = @(y) 1 / (pi * w0^2) * integral2(@(beta, r) r .* gauss(r, w0, r_max) .* ...
     iscomplex(qg_y_factor(beta, r, y, n1, n2, Rsp, f)), 0, 2*pi, 0, r_max, ...
     'Method', 'iterated', 'AbsTol', 1e-8, 'RelTol', 1e-6);
 
-Qres_s = @(y) G0 * integral2(@(beta, r) r .* gauss(r, r_max, w0) .* ...
+Qres_s = @(y) 1 / (pi * w0^2) * integral2(@(beta, r) r .* gauss(r, w0, r_max) .* ...
     iscomplex(qs_y_factor(beta, r, y, n1, n2, Rsp, f)), 0, 2*pi, 0, r_max, ...
     'Method', 'iterated', 'AbsTol', 1e-8, 'RelTol', 1e-6);