diff --git a/Matlab/bessel.m b/Matlab/bessel.m
index 2561504..413aaea 100644
--- a/Matlab/bessel.m
+++ b/Matlab/bessel.m
@@ -1,4 +1,4 @@
 % Bessel beam
 function b = bessel(r, r_max, w0, P)
-    b = bessel_peak(r_max, w0, P) * besselj(0, 2.405/w0 * r).^2;
+    b = besselj(0, 2.405/w0 * r).^2;
 end
\ No newline at end of file
diff --git a/Python/functions.py b/Python/functions.py
new file mode 100644
index 0000000..c3d0d68
--- /dev/null
+++ b/Python/functions.py
@@ -0,0 +1,128 @@
+import numpy as np
+from scipy import special
+from scipy import integrate
+from scipy import constants
+
+n1 = 1.3337  # index of refraction of the immersion medium
+n2 = 1.4607  # index of refraction of the fused silica at wavelength 523 nm
+NA = 1.25  # numerical aperture
+f = 2.0e-3  # objective lens focus or WD
+Rsp = 1.03e-6  # sphere radius
+P = 4.4e-3  # power of the laser
+ratio = 1.0  # the ratio of the beam radius to the aperture radius
+
+th_max = np.arcsin(NA / n1)  # maximum angle of incidence
+# radius of a Gaussian beam (1:1 with input aperture condition)
+r_max = f * np.tan(th_max)
+n = n2 / n1  # n2/n1
+w0 = ratio * r_max  # beam radius
+F0 = n1 * P / constants.speed_of_light; # resulting force
+
+
+# Angle of refraction
+def theta_r(th):
+    return np.arcsin(n1 / n2 * np.sin(th))
+
+
+# Fresnel reflectivity
+def reflectivity(th, psi):
+    theta = theta_r(th)
+    return (np.tan(th - theta) ** 2 / np.tan(th + theta) ** 2) * np.cos(psi) ** 2 + \
+           (np.sin(th - theta) ** 2 / np.sin(th + theta) ** 2) * np.sin(psi) ** 2
+
+
+# Fresnel transparency
+def transparency(th, psi):
+    return 1 - reflectivity(th, psi)
+
+
+# force factors
+def q_s(th, psi):
+    R = reflectivity(th, psi)
+    T = transparency(th, psi)
+    theta = theta_r(th)
+    return 1 + R * np.cos(2 * th) - T ** 2 * \
+        (np.cos(2 * th - 2 * theta) + R * np.cos(2*th)) / \
+        (1 + R ** 2 + 2 * R * np.cos(2*theta))
+
+
+def q_g(th, psi):
+    R = reflectivity(th, psi)
+    T = transparency(th, psi)
+    theta = theta_r(th)
+    return R * np.sin(2 * th) - T ** 2 * \
+        (np.sin(2 * th - 2 * theta) + R * np.sin(2 * th)) / \
+        (1 + R ** 2 + 2 * R * np.cos(2 * theta))
+
+
+def q_mag(th, psi):
+    return np.sqrt(q_s(th, psi) ** 2 + q_g(th, psi) ** 2)
+
+
+# Average factors (circular polarization)
+def q_s_avg(th):
+    return 0.5 * (q_s(th, 0) + q_s(th, np.pi/2))
+
+
+def q_g_avg(th):
+    return 0.5 * (q_g(th, 0) + q_g(th, np.pi/2))
+
+
+def q_mag_avg(th):
+    return np.sqrt(q_s_avg(th) ** 2 + q_g_avg(th) ** 2)
+
+
+# Angles
+def phi_i(r):
+    return np.arctan(r / f)
+
+
+def gamma(beta, r):
+    return np.arccos(np.cos(np.pi / 2 - phi_i(r)) * np.cos(beta), dtype=np.cfloat)
+
+
+def th_i_z(r, z):
+    return np.arcsin(z / Rsp * np.sin(phi_i(r)), dtype=np.cfloat)
+
+
+def th_i_y(beta, r, y):
+    return np.arcsin(y / Rsp * np.sin(gamma(beta, r)), dtype=np.cfloat)
+
+
+def q_g_z(r, z):
+    return -q_g_avg(th_i_z(r, z)) * np.sin(phi_i(r))
+
+
+def q_s_z(r, z):
+    return q_s_avg(th_i_z(r, z)) * np.cos(phi_i(r))
+
+
+def q_g_y(beta, r, y):
+    return q_g_avg(th_i_y(beta, r, y)) * np.cos(phi_i(r), dtype=np.cfloat)
+
+
+def q_s_y(beta, r, y):
+    return q_s_avg(th_i_y(beta, r, y)) * np.sin(gamma(beta, r), dtype=np.cfloat)
+
+
+# Intensity distributions
+def gauss_peak():
+    A = (1 - np.exp(-2*r_max ** 2 / w0 ** 2))
+    return 2*A / (np.pi * w0 ** 2)
+
+
+def gauss(r):
+    return np.exp(-2 * r ** 2 / w0 ** 2)
+
+
+def bessel(r):
+    return special.jv(0, 2.405 / w0 * r) ** 2
+
+
+def bessel_peak():
+    return 4.81 / (w0 * r_max * np.exp(0.5)) * 2 * np.pi * integrate.quad(lambda r: r * special.jv(0, 2.405 / w0 * r) ** 2,
+                                                                          0, r_max, epsabs=1e-12, epsrel=1e-6)[0]
+
+
+def uniform(r):
+    return P / (np.pi * r_max ** 2) * np.ones(r.shape)
diff --git a/Python/trap_forces_axial.py b/Python/trap_forces_axial.py
index 20b0542..d341d0c 100644
--- a/Python/trap_forces_axial.py
+++ b/Python/trap_forces_axial.py
@@ -5,134 +5,42 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import seaborn as sns
+from functions import *
 from scipy import integrate
-from scipy import special
-from scipy import constants
-
-n1 = 1.3337  # index of refraction of the immersion medium
-n2 = 1.4607 # index of refraction of the fused silica at wavelength 523 nm
-n = n2 / n1 # n2/n1
-NA = 1.25  # numerical aperture
-th_max = np.arcsin(NA / n1)  # maximum angle of incidence
-f = 2.0e-3  # objective lens focus or WD
-r_max = f * np.tan(th_max)  # radius of a Gaussian beam (1:1 with input aperture condition)
-Rsp = 1.03e-6  # sphere radius
-P = 4.4e-3  # power of the laser
 
 
-# angle of refraction
-def r(th):
-    return np.arcsin(n1 / n2 * np.sin(th))
-
-
-# Fresnel reflectivity
-def r_f(th, psi):
-    return (np.tan(th - r(th)) ** 2 / np.tan(th + r(th)) ** 2) * np.cos(psi) ** 2 + \
-           (np.sin(th - r(th)) ** 2 / np.sin(th + r(th)) ** 2) * np.sin(psi) ** 2
-
-
-# Fresnel transparency
-def t_f(th, psi):
-    return 1 - r_f(th, psi)
-
-
-# force factors
-def q_s(th, psi):
-    return 1 + r_f(th, psi) * np.cos(2 * th) - t_f(th, psi) ** 2 * \
-           (np.cos(2 * th - 2 * r(th)) + r_f(th, psi) * np.cos(2*th)) / \
-           (1 + r_f(th, psi) ** 2 + 2 * r_f(th, psi) * np.cos(2*r(th)))
-
-
-def q_g(th, psi):
-    return r_f(th, psi) * np.sin(2 * th) - t_f(th, psi) ** 2 * \
-           (np.sin(2 * th - 2 * r(th)) + r_f(th, psi) * np.sin(2 * th)) / \
-           (1 + r_f(th, psi) ** 2 + 2 * r_f(th, psi) * np.cos(2 * r(th)))
-
-
-def q_mag(th, psi):
-    return np.sqrt(q_s(th, psi) ** 2 + q_g(th, psi) ** 2)
-
-
-# Average factors (circular polarization
-def q_s_avg(th):
-    return 0.5 * (q_s(th, 0) + q_s(th, np.pi/2))
-
-
-def q_g_avg(th):
-    return 0.5 * (q_g(th, 0) + q_g(th, np.pi/2))
-
-
-def q_mag_avg(th):
-    return np.sqrt(q_s_avg(th) ** 2 + q_g_avg(th) ** 2)
-
-
-# Angles
-def phi(dr):
-    return np.arctan(dr / f)
-
-
-def thi(dr, dz):
-    return np.arcsin(dz / Rsp * np.sin(phi(dr)), dtype=np.cfloat)
-
-
-def q_g_z(dr, dz):
-    return -q_g_avg(thi(dr, dz)) * np.sin(phi(dr))
-
-
-def q_s_z(dr, dz):
-    return q_s_avg(thi(dr, dz)) * np.cos(phi(dr))
-
-
-# Intensity profile
-a = 1.0
-w0 = a * r_max
-
-
-def intensity_uniform():
-    return P / (np.pi * r_max ** 2)
-
-
-def intensity_gaussian_tem00(dr):
-    i_0 = P * 2 / (np.pi*w0 ** 2)
-    return i_0 * np.exp(-2 * dr ** 2 / w0 ** 2)
-
-
-def intensity_bessel(dr):
-    w0_bb = 0.5 * r_max
-    i_0 = P * 2 / (np.pi * w0 ** 2)
-    return i_0 * special.jv(0, 2.405 / w0_bb * dr) ** 2 * np.exp(- 2 * dr ** 2 / w0 ** 2)
+sns.set()
 
 
 # Intensity profile graphics
-sns.set()
-sns.set_style("darkgrid")
-
 rho = np.linspace(-r_max, r_max, 500)
 fig1 = plt.figure(1, figsize=(10, 6))
-plt.plot(rho, intensity_gaussian_tem00(rho), 'k')
+I = gauss(rho)
+I0 = np.max(I)
+plt.plot(rho, I / I0)
+plt.fill_between(rho, I / I0, 0, alpha=0.3)
 plt.xlabel('r, m', fontsize=18)
 plt.ylabel('I(r)', fontsize=18)
 
-
 # Integration
-def q_res_g(dz, func):
-    ans = integrate.quad(lambda x: x * func(x) * q_g_z(x, dz) * (~np.iscomplex(q_g_z(x, dz))).astype(float),
-                               0, r_max, epsabs=1e-12, epsrel=1e-6)
-    return 2 * np.pi * ans[0]
+def q_res_g(z, func):
+    ans = integrate.quad(lambda x: x * func(x) * q_g_z(x, z) * (~np.iscomplex(q_g_z(x, z))).astype(float),
+                         0, r_max, epsabs=1e-12, epsrel=1e-6)
+    return ans[0]
 
 
-def q_res_s(dz, func):
-    ans = integrate.quad(lambda x: x * func(x) * q_s_z(x, dz) * (~np.iscomplex(q_s_z(x, dz))).astype(float),
-                               0, r_max, epsabs=1e-12, epsrel=1e-6)
-    return 2 * np.pi * ans[0]
+def q_res_s(z, func):
+    ans = integrate.quad(lambda x: x * func(x) * q_s_z(x, z) * (~np.iscomplex(q_s_z(x, z))).astype(float),
+                         0, r_max, epsabs=1e-12, epsrel=1e-6)
+    return ans[0]
 
 
 # Calculation
 n = 200
 z = np.linspace(-2 * Rsp, 2 * Rsp, n)
 
-axial_g = [q_res_g(x, intensity_gaussian_tem00) for x in z]
-axial_s = [q_res_s(x, intensity_gaussian_tem00) for x in z]
+axial_g = gauss_peak() * [q_res_g(x, gauss) for x in z]
+axial_s = gauss_peak() * [q_res_s(x, gauss) for x in z]
 
 f_0 = n1 * P / constants.c  # net force
 
@@ -143,9 +51,9 @@ z = -z[::-1]
 
 # Graphics
 fig2 = plt.figure(2, figsize=(10, 6))
-plt.plot(z, f_0*axial_g, 'b-.', lw=1, label='$F_{g}$')
-plt.plot(z, f_0*axial_s, 'r--', lw=1, label='$F_{s}$')
-plt.plot(z, f_0*axial, 'k', lw=1, label='$F_{t}$')
+plt.plot(z, f_0*axial_g, '-.', lw=1, label='$F_{g}$')
+plt.plot(z, f_0*axial_s, '--', lw=1, label='$F_{s}$')
+plt.plot(z, f_0*axial, lw=1, label='$F_{t}$')
 plt.xlabel('r, m', fontsize=18)
 plt.ylabel('F, N', fontsize=18)
 plt.legend(fontsize=18)
diff --git a/Python/trap_forces_efficiencies.py b/Python/trap_forces_efficiencies.py
index 0d5d521..9a27e0e 100644
--- a/Python/trap_forces_efficiencies.py
+++ b/Python/trap_forces_efficiencies.py
@@ -1,62 +1,22 @@
 # These calculations are based on Ashkin's article "Forces of a single-beam
 # gradient laser trap on a dielectric sphere in the ray optics regime
 
-# all distances in mm
-
 import numpy as np
 import matplotlib.pyplot as plt
+import seaborn as sns
+from functions import *
 
 
-a = 1.0 * 1e-6  # radius of the bead
-n1 = 1.33  # index of refraction of the medium
-n = 1.8  # n2/n1
-n2 = n * n1  # index of refraction of the fused silica
-c0 = 3 * 1e8  # speed of light
-
-
-# Fresnel reflectivity
-def r_f(th, psi):
-    return (np.tan(th - np.arcsin(n1 / n2 * np.sin(th))) ** 2 /
-            np.tan(th + np.arcsin(n1 / n2 * np.sin(th))) ** 2) * np.cos(psi) ** 2 + \
-           (np.sin(th - np.arcsin(n1 / n2 * np.sin(th))) ** 2 /
-            np.sin(th + np.arcsin(n1 / n2 * np.sin(th))) ** 2) * np.sin(psi) ** 2
-
-
-# Fresnel transparency
-def t_f(th, psi):
-    return 1 - r_f(th, psi)
-
-
-# angle of refraction
-def r(th):
-    return np.arcsin(n1 / n2 * np.sin(th))
-
-
-# force factors
-def q_s(th, psi):
-    return 1 + r_f(th, psi) * np.cos(2 * th) - t_f(th, psi) ** 2 * \
-           (np.cos(2 * th - 2 * r(th)) + r_f(th, psi) * np.cos(2*th)) / \
-           (1 + r_f(th, psi) ** 2 + 2 * r_f(th, psi) * np.cos(2*r(th)))
-
-
-def q_g(th, psi):
-    return r_f(th, psi) * np.sin(2 * th) - t_f(th, psi) ** 2 * \
-           (np.sin(2 * th - 2 * r(th)) + r_f(th, psi) * np.sin(2 * th)) / \
-           (1 + r_f(th, psi) ** 2 + 2 * r_f(th, psi) * np.cos(2 * r(th)))
-
-
-def q_mag(th, psi):
-    return np.sqrt(q_s(th, psi) ** 2 + q_g(th, psi) ** 2)
-
+sns.set()
 
 t = np.linspace(0, np.pi / 2, 1000)
 t_deg = t * 180 / np.pi
 pol = np.pi / 4
 
 plt.figure(figsize=(13, 8))
-plt.plot(t_deg, q_s(t, pol), 'r--', label='Q_s')
-plt.plot(t_deg, -q_g(t, pol), 'b-.', label='Q_g')
-plt.plot(t_deg, q_mag(t, pol), 'k', label='Q_t')
+plt.plot(t_deg, q_s(t, pol), '--', label='Q_s')
+plt.plot(t_deg, -q_g(t, pol), '-.', label='Q_g')
+plt.plot(t_deg, q_mag(t, pol), label='Q_t')
 plt.grid()
 plt.xlabel(r'$\theta$, deg', fontsize=18)
 plt.ylabel('Q', fontsize=18)
diff --git a/Python/trap_forces_transverse.py b/Python/trap_forces_transverse.py
index 9d9cae2..30f5ad7 100644
--- a/Python/trap_forces_transverse.py
+++ b/Python/trap_forces_transverse.py
@@ -4,133 +4,36 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy import integrate
-from scipy import special
-from scipy import constants
-
-n1 = 1.3337  # index of refraction of the immersion medium
-n2 = 1.4607 # index of refraction of the fused silica at wavelength 523 nm
-n = n2 / n1 # n2/n1
-NA = 1.25  # numerical aperture
-th_max = np.arcsin(NA / n1)  # maximum angle of incidence
-f = 2.0e-3  # objective lens focus or WD
-r_max = f * np.tan(th_max)  # radius of a Gaussian beam (1:1 with input aperture condition)
-Rsp = 1.03e-6  # sphere radius
-P = 14e-3  # power of the laser
+import seaborn as sns
+from functions import *
 
 
-# Angles
-def r(th):
-    return np.arcsin(n1 / n2 * np.sin(th, dtype=np.cfloat))
-
-
-def phi(dr):
-    return np.arctan(dr / f, dtype=np.cfloat)
-
-
-def gamma(db, dr):
-    return np.arccos(np.cos(np.pi / 2 - phi(dr)) * np.cos(db), dtype=np.cfloat)
-
-
-def thi(db, dr, dy):
-    return np.arcsin(dy / Rsp * np.sin(gamma(db, dr)), dtype=np.cfloat)
-
-
-# Fresnel reflectivity
-def r_f(th, psi):
-    return (np.tan(th - r(th)) ** 2 / np.tan(th + r(th)) ** 2) * np.cos(psi) ** 2 + \
-           (np.sin(th - r(th)) ** 2 / np.sin(th + r(th)) ** 2) * np.sin(psi) ** 2
-
-
-# Fresnel transparency
-def t_f(th, psi):
-    return 1 - r_f(th, psi)
-
-
-# force factors
-def q_s(th, psi):
-    return 1 + r_f(th, psi) * np.cos(2 * th) - t_f(th, psi) ** 2 * \
-           (np.cos(2 * th - 2 * r(th)) + r_f(th, psi) * np.cos(2*th)) / \
-           (1 + r_f(th, psi) ** 2 + 2 * r_f(th, psi) * np.cos(2*r(th)))
-
-
-def q_g(th, psi):
-    return r_f(th, psi) * np.sin(2 * th) - t_f(th, psi) ** 2 * \
-           (np.sin(2 * th - 2 * r(th)) + r_f(th, psi) * np.sin(2 * th)) / \
-           (1 + r_f(th, psi) ** 2 + 2 * r_f(th, psi) * np.cos(2 * r(th)))
-
-
-def q_mag(th, psi):
-    return np.sqrt(q_s(th, psi) ** 2 + q_g(th, psi) ** 2)
-
-
-# Average factors (circular polarization)
-def q_s_avg(th):
-    return 0.5 * (q_s(th, 0) + q_s(th, np.pi/2))
-
-
-def q_g_avg(th):
-    return 0.5 * (q_g(th, 0) + q_g(th, np.pi/2))
-
-
-def q_mag_avg(th):
-    return np.sqrt(q_s_avg(th) ** 2 + q_g_avg(th) ** 2)
-
-
-def q_g_z(db, dr, dy):
-    return q_g_avg(thi(db, dr, dy)) * np.cos(phi(dr), dtype=np.cfloat)
-
-
-def q_s_z(db, dr, dy):
-    return q_s_avg(thi(db, dr, dy)) * np.sin(gamma(db, dr), dtype=np.cfloat)
-
-
-# Intensity profile
-a = 1.0
-w0 = a * r_max
-
-
-def intensity_uniform(dr):
-    return 1 / (np.pi * r_max ** 2)
-
-
-def intensity_gaussian_tem00(dr):
-    amp = (1 - np.exp(-2 * r_max ** 2 / w0 ** 2))
-    p_0 = np.pi * w0 ** 2 / 2
-    return 1 / (amp * p_0) * np.exp(-2 * dr ** 2 / w0 ** 2)
-
-
-def intensity_bessel(dr):
-    p_0 = np.exp(0.5) * w0 * 0.0025 / 4.81
-
-    def temp_f(w):
-        return w * special.jv(0, 2.405 / w0 * w) ** 2
-    amp = integrate.quad(temp_f, 0, r_max)
-    return 1 / (2 * np.pi * amp[0]) * special.jv(0, 2.405 / w0 * dr) ** 2
+sns.set()
 
 
 # Intensity profile graphics
 rho = np.linspace(-r_max, r_max, 500)
 fig1 = plt.figure(1, figsize=(10, 6))
-plt.plot(rho, intensity_gaussian_tem00(rho), 'k')
-plt.grid()
+I = gauss(rho)
+I0 = np.max(I)
+plt.plot(rho, I / I0)
+plt.fill_between(rho, I / I0, 0, alpha=0.3)
 plt.xlabel('r, m', fontsize=18)
 plt.ylabel('I(r)', fontsize=18)
-plt.show()
 
 
 # Integration
 def q_res_g(dy, func):
-    ans = integrate.dblquad(lambda dr, db, dy: dr * func(dr) * q_g_z(db, dr, dy) * (~np.iscomplex(q_g_z(db, dr, P))).astype(float),
+    ans = integrate.dblquad(lambda dr, db, dy: dr * func(dr) * q_g_y(db, dr, dy) * (~np.iscomplex(q_g_y(db, dr, dy))).astype(float),
                             0, 2 * np.pi, 0, r_max, args=(dy, ),
-                            epsabs=1e-4, epsrel=1e-6)
+                            epsabs=1e-8, epsrel=1e-6)
     return ans[0]
 
 
 def q_res_s(dy, func):
-    ans = integrate.dblquad(lambda dr, db: dr * func(dr) * q_s_z(db, dr, dy) * (~np.iscomplex(q_g_z(db, dr, P))).astype(float),
+    ans = integrate.dblquad(lambda dr, db: dr * func(dr) * q_s_y(db, dr, dy) * (~np.iscomplex(q_g_y(db, dr, dy))).astype(float),
                             0, 2 * np.pi, 0, r_max,
-                            epsabs=1e-4, epsrel=1e-6)
+                            epsabs=1e-8, epsrel=1e-6)
     return ans[0]
 
 
@@ -138,21 +41,17 @@ def q_res_s(dy, func):
 n = 150
 y = np.linspace(-2 * Rsp, 2 * Rsp, n)
 
-transverse_g = np.abs(np.array([q_res_g(x, intensity_gaussian_tem00) for x in y]))
-transverse_s = np.array([q_res_s(x, intensity_gaussian_tem00) for x in y])
-
-f_0 = n1 * P / constants.c  # net force
-
+transverse_g = gauss_peak() * F0 * np.abs(np.array([q_res_g(x, gauss) for x in y]))
+transverse_s = gauss_peak() * F0 * np.array([q_res_s(x, gauss) for x in y])
 transverse = transverse_g + transverse_s
 
 
 # Graphics
 fig2 = plt.figure(2, figsize=(10, 6))
-plt.plot(y, transverse_g, 'b-.', lw=1, label='$F_{g}$')
-plt.plot(y, transverse_s, 'r--', lw=1, label='$F_{s}$')
-plt.plot(y, transverse, 'k', lw=1, label='$F_{t}$')
+plt.plot(y, transverse_g, '-.', lw=1, label='$F_{g}$')
+plt.plot(y, transverse_s, '--', lw=1, label='$F_{s}$')
+plt.plot(y, transverse, lw=1, label='$F_{t}$')
 plt.xlabel('r, m', fontsize=18)
 plt.ylabel('F, m', fontsize=18)
 plt.legend()
-plt.grid()
 plt.show()
\ No newline at end of file