From 95fa3c4f3fa58e7501a5dc33f355c38f37f97d9a Mon Sep 17 00:00:00 2001 From: BrokenVoodooDoll Date: Wed, 22 Mar 2023 15:07:19 +0400 Subject: [PATCH] Fixed python calculations, added gitignore --- .gitignore | 5 +++ Matlab/bessel.m | 4 +- Matlab/test.m | 16 +++++++ Matlab/{axial.m => trap_forces_axial.m} | 12 +++--- ...iciencies.m => trap_forces_efficiencies.m} | 0 ...{transverse.m => trap_forces_transverse.m} | 0 Matlab/uniform.m | 4 +- Python/functions.py | 42 ++++++++++++------- Python/trap_forces_axial.py | 25 +++++------ Python/trap_forces_transverse.py | 13 +++--- README.md | 2 +- 11 files changed, 77 insertions(+), 46 deletions(-) create mode 100644 .gitignore create mode 100644 Matlab/test.m rename Matlab/{axial.m => trap_forces_axial.m} (68%) rename Matlab/{ray_efficiencies.m => trap_forces_efficiencies.m} (100%) rename Matlab/{transverse.m => trap_forces_transverse.m} (100%) diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f1d9e68 --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +# Python +__pycache__/ +*.pyc +# Matlab +*.asv \ No newline at end of file diff --git a/Matlab/bessel.m b/Matlab/bessel.m index 281e296..78cf547 100644 --- a/Matlab/bessel.m +++ b/Matlab/bessel.m @@ -1,6 +1,6 @@ % Bessel beam function b = bessel(r, w0, r_max) - ring_radius = 2.405; % radisu of the first ring of the besselj_0 + ring_radius = 2.405; % radius 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; + b = A * besselj(0, ring_radius/w0 * r).^2; end \ No newline at end of file diff --git a/Matlab/test.m b/Matlab/test.m new file mode 100644 index 0000000..a100e02 --- /dev/null +++ b/Matlab/test.m @@ -0,0 +1,16 @@ +load_constants + +fprintf('q_s = %e\n', qs_factor(pi / 4, pi / 4, n1, n2)) +fprintf('q_g = %e\n', qg_factor(pi / 4, pi / 4, n1, n2)) +fprintf('q_mag = %e\n', qmag_factor(pi / 4, pi / 4, n1, n2)) +fprintf('qs_avg = %e\n', qs_avg_factor(pi / 4, n1, n2)) +fprintf('qg_avg = %e\n', qg_avg_factor(pi / 4, n1, n2)) +fprintf('q_mag_avg = %e\n', qmag_avg_factor(pi / 4, n1, n2)) +fprintf('phi_i = %e\n', phi_i(Rsp, f)) +fprintf('gamma = %e\n', gamma_angle(pi / 4, Rsp, f)) +fprintf('th_i_z = %e\n', th_i_z(Rsp, Rsp, Rsp, f)) +fprintf('th_i_y = %e\n', th_i_y(pi / 4, Rsp, Rsp, Rsp, f)) +fprintf('qg_z = %e\n', qg_z_factor(Rsp, Rsp, n1, n2, Rsp, f)) +fprintf('qs_z = %e\n', qs_z_factor(Rsp, Rsp, n1, n2, Rsp, f)) +fprintf('qg_y = %e\n', qg_y_factor(pi / 4, Rsp, Rsp, n1, n2, Rsp, f)) +fprintf('qs_y = %e\n', qs_y_factor(pi / 4, Rsp, Rsp, n1, n2, Rsp, f)) \ No newline at end of file diff --git a/Matlab/axial.m b/Matlab/trap_forces_axial.m similarity index 68% rename from Matlab/axial.m rename to Matlab/trap_forces_axial.m index 41af23a..7c44dad 100644 --- a/Matlab/axial.m +++ b/Matlab/trap_forces_axial.m @@ -20,13 +20,13 @@ xlabel('r, ΠΌ') ylabel('I(r)') % Integration -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_g = @(z) 2*pi / (pi * w0 ^ 2) * integral(@(r) r .* gauss(r, w0, r_max) .* ... + iscomplex(qg_z_factor(r, z, n1, n2, Rsp, f)), 0, r_max, ... + 'AbsTol', 1e-12, 'RelTol', 1e-6); -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); +Qres_s = @(z) 2*pi / (pi * w0 ^ 2) * integral(@(r) r .* gauss(r, w0, r_max) .* ... + iscomplex(qs_z_factor(r, z, n1, n2, Rsp, f)), 0, r_max, ... + 'AbsTol', 1e-12, 'RelTol', 1e-6); % Calulation N = 200; diff --git a/Matlab/ray_efficiencies.m b/Matlab/trap_forces_efficiencies.m similarity index 100% rename from Matlab/ray_efficiencies.m rename to Matlab/trap_forces_efficiencies.m diff --git a/Matlab/transverse.m b/Matlab/trap_forces_transverse.m similarity index 100% rename from Matlab/transverse.m rename to Matlab/trap_forces_transverse.m diff --git a/Matlab/uniform.m b/Matlab/uniform.m index 7a10840..a4bca3e 100644 --- a/Matlab/uniform.m +++ b/Matlab/uniform.m @@ -1,4 +1,4 @@ % Beam with uniform distribution -function u = uniform(r, w0, r_max, P) - u = P / (pi*r_max^2); +function u = uniform(r, w0, r_max) + u = ones(size(r)); end \ No newline at end of file diff --git a/Python/functions.py b/Python/functions.py index c3d0d68..7bb567f 100644 --- a/Python/functions.py +++ b/Python/functions.py @@ -8,7 +8,7 @@ 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 +P = 51.7e-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 @@ -106,23 +106,37 @@ def q_s_y(beta, r, y): # 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) + # the fraction of power that falls on the pupil of the micro lens + A = (1 - np.exp(-2*r_max ** 2 / w0 ** 2)) + return 2 * A * 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] + ring_radius = 2.405; # radius of the first ring of the besselj_0 + # the fraction of power that falls on the pupil of the micro lens + A = ring_radius / (w0 * r_max * np.exp(0.5)) * 2 * np.pi * \ + integrate.quad(lambda r: r * special.jv(0, ring_radius / w0 * r) ** 2, + 0, r_max, epsabs=1e-12, epsrel=1e-6)[0] + return 2 * A* special.jv(0, ring_radius / w0 * r) ** 2 def uniform(r): - return P / (np.pi * r_max ** 2) * np.ones(r.shape) + return np.ones(r.shape) + +def test(): + print("q_s = ", q_s(np.pi / 4, np.pi / 4)) + print("q_g = ", q_g(np.pi / 4, np.pi / 4)) + print("q_mag = ", q_mag(np.pi / 4, np.pi / 4)) + print("q_s_avg = ", q_s_avg(np.pi / 4)) + print("q_g_avg = ", q_g_avg(np.pi / 4)) + print("q_mag_avg = ", q_mag_avg(np.pi / 4)) + print("phi_i = ", phi_i(Rsp)) + print("gamma = ", gamma(np.pi / 4, Rsp)) + print("th_i_z = ", th_i_z(Rsp, Rsp)) + print("th_i_y = ", th_i_y(np.pi / 4, Rsp, Rsp)) + print("q_g_z = ", q_g_z(Rsp, Rsp)) + print("q_s_z = ", q_s_z(Rsp, Rsp)) + print("q_g_y = ", q_g_y(np.pi / 4, Rsp, Rsp)) + print("q_s_y = ", q_s_y(np.pi / 4, Rsp, Rsp)) + \ No newline at end of file diff --git a/Python/trap_forces_axial.py b/Python/trap_forces_axial.py index d341d0c..c5c32d5 100644 --- a/Python/trap_forces_axial.py +++ b/Python/trap_forces_axial.py @@ -24,36 +24,33 @@ plt.ylabel('I(r)', fontsize=18) # Integration 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), + ans = integrate.quad(lambda dr: dr * func(dr) * q_g_z(dr, z) * (~np.iscomplex(q_g_z(dr, z))).astype(float), 0, r_max, epsabs=1e-12, epsrel=1e-6) - return ans[0] - + return 2 * np.pi * ans[0] / (np.pi * w0 ** 2) 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), + ans = integrate.quad(lambda dr: dr * func(dr) * q_s_z(dr, z) * (~np.iscomplex(q_s_z(dr, z))).astype(float), 0, r_max, epsabs=1e-12, epsrel=1e-6) - return ans[0] + return 2 * np.pi * ans[0] / (np.pi * w0 ** 2) # Calculation n = 200 z = np.linspace(-2 * Rsp, 2 * Rsp, n) -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] +axial_g = [q_res_g(dz, gauss) for dz in z] +axial_s = [q_res_s(dz, gauss) for dz in z] -f_0 = n1 * P / constants.c # net force - -axial_g = np.array(axial_g[::-1]) -axial_s = np.array(axial_s[::-1]) +axial_g = F0 * np.array(axial_g[::-1]) +axial_s = F0 * np.array(axial_s[::-1]) axial = axial_g + axial_s z = -z[::-1] # Graphics fig2 = plt.figure(2, figsize=(10, 6)) -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.plot(z, axial_g, '-.', lw=1, label='$F_{g}$') +plt.plot(z, axial_s, '--', lw=1, label='$F_{s}$') +plt.plot(z, 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_transverse.py b/Python/trap_forces_transverse.py index 30f5ad7..b0bed50 100644 --- a/Python/trap_forces_transverse.py +++ b/Python/trap_forces_transverse.py @@ -26,23 +26,22 @@ plt.ylabel('I(r)', fontsize=18) def q_res_g(dy, func): 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-8, epsrel=1e-6) - return ans[0] + epsabs=1e-6, epsrel=1e-6) + return ans[0] / (np.pi * w0 ** 2) def q_res_s(dy, func): 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-8, epsrel=1e-6) - return ans[0] - + epsabs=1e-6, epsrel=1e-6) + return ans[0] / (np.pi * w0 ** 2) # Calculation n = 150 y = np.linspace(-2 * Rsp, 2 * Rsp, n) -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_g = F0 * np.abs(np.array([q_res_g(dy, gauss) for dy in y])) +transverse_s = F0 * np.array([q_res_s(dy, gauss) for dy in y]) transverse = transverse_g + transverse_s diff --git a/README.md b/README.md index cb184bf..d7d4f3a 100644 --- a/README.md +++ b/README.md @@ -1 +1 @@ -# optical-trap-forces \ No newline at end of file +# Optical Trap Forces \ No newline at end of file