From 74231472c916632ce4d1b0ae4b7b3532e139a894 Mon Sep 17 00:00:00 2001 From: BrokenVoodooDoll Date: Mon, 20 Mar 2023 22:50:56 +0400 Subject: [PATCH] Refactoring --- Matlab/Forces_Ashkin_Transverse.m | 106 ------------------------------ Matlab/angles.m | 21 ------ Matlab/axial.asv | 70 ++++++++++++++++++++ Matlab/axial.m | 33 ++++++---- Matlab/bessel.m | 4 ++ Matlab/bessel_peak.m | 3 + Matlab/coefficients.m | 9 --- Matlab/factors.m | 37 ----------- Matlab/gamma_angle.m | 3 + Matlab/gauss.m | 4 ++ Matlab/gauss_peak.m | 4 ++ Matlab/intensity.m | 7 -- Matlab/load_constants.m | 15 +++-- Matlab/phi_i.m | 3 + Matlab/qg_avg_factor.m | 3 + Matlab/qg_factor.m | 8 +++ Matlab/qg_y_factor.m | 3 + Matlab/qg_z_factor.m | 3 + Matlab/qmag_avg_factor.m | 3 + Matlab/qmag_factor.m | 3 + Matlab/qs_avg_factor.m | 3 + Matlab/qs_factor.m | 9 +++ Matlab/qs_y_factor.m | 3 + Matlab/qs_z_factor.m | 3 + Matlab/reflectivity.m | 6 ++ Matlab/th_i_y.m | 3 + Matlab/th_i_z.m | 3 + Matlab/th_r.m | 3 + Matlab/theta_angle.m | 3 + Matlab/transmittance.m | 3 + Matlab/transverse.m | 56 ++++++++++++++++ Matlab/uniform.m | 4 ++ 32 files changed, 240 insertions(+), 201 deletions(-) delete mode 100644 Matlab/Forces_Ashkin_Transverse.m delete mode 100644 Matlab/angles.m create mode 100644 Matlab/axial.asv create mode 100644 Matlab/bessel.m create mode 100644 Matlab/bessel_peak.m delete mode 100644 Matlab/coefficients.m delete mode 100644 Matlab/factors.m create mode 100644 Matlab/gamma_angle.m create mode 100644 Matlab/gauss.m create mode 100644 Matlab/gauss_peak.m delete mode 100644 Matlab/intensity.m create mode 100644 Matlab/phi_i.m create mode 100644 Matlab/qg_avg_factor.m create mode 100644 Matlab/qg_factor.m create mode 100644 Matlab/qg_y_factor.m create mode 100644 Matlab/qg_z_factor.m create mode 100644 Matlab/qmag_avg_factor.m create mode 100644 Matlab/qmag_factor.m create mode 100644 Matlab/qs_avg_factor.m create mode 100644 Matlab/qs_factor.m create mode 100644 Matlab/qs_y_factor.m create mode 100644 Matlab/qs_z_factor.m create mode 100644 Matlab/reflectivity.m create mode 100644 Matlab/th_i_y.m create mode 100644 Matlab/th_i_z.m create mode 100644 Matlab/th_r.m create mode 100644 Matlab/theta_angle.m create mode 100644 Matlab/transmittance.m create mode 100644 Matlab/transverse.m create mode 100644 Matlab/uniform.m diff --git a/Matlab/Forces_Ashkin_Transverse.m b/Matlab/Forces_Ashkin_Transverse.m deleted file mode 100644 index d5ac37b..0000000 --- a/Matlab/Forces_Ashkin_Transverse.m +++ /dev/null @@ -1,106 +0,0 @@ -%all distances in m -close all -clear -clc -format compact - -% 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". -% There are transverse forces only - -n1 = 1.33; % index of refraction of the immersion medium -n2 = 1.6; % index of refraction of the fused silica at wavelength 523 nm -n = n2/n1; % n2/n1 -c0 = 3e8; % speed of light -NA = 1.25; % numerical aperture -th_max = asin(NA/n1); % maximum angle of incidence -f = 100.0e-3; % objective lens focus or WD -r_max = f*tan(th_max); % radius of a Gaussian beam (1:1 with input aperture condition) -Rsp = 1.0e-6; % sphere radius -P = 20.0e-3; % power of the laser - -thr = @(th) asin(n1/n2*sin(th)); % refraction angle - -%reflectivity -R = @(th,psi) (tan(th-thr(th)).^2./tan(th+thr(th)).^2).*cos(psi).^2+... - (sin(th-thr(th)).^2./sin(th+thr(th)).^2).*sin(psi).^2; - -%transparency -T = @(th,psi) 1-R(th,psi); - -% Factors -Qs = @(th, psi) 1 + R(th, psi) .* cos(2*th) - T(th,psi).^2 .* (cos(2*th -... - 2*thr(th)) + R(th, psi) .* cos(2*th)) ./ (1 + R(th,psi).^2 +... - 2*R(th,psi) .* cos(2*thr(th))); - -Qg = @(th, psi) R(th, psi) .* sin(2*th) - T(th,psi).^2 .* (sin(2*th -... - 2*thr(th)) + R(th, psi) .* sin(2*th)) ./ (1 + R(th,psi).^2 +... - 2*R(th,psi) .* cos(2*thr(th))); - -% Average factors (circular polarization -Qs_avg = @(th) 0.5*(Qs(th, 0) + Qs(th, pi/2)); -Qg_avg = @(th) 0.5*(Qg(th, 0) + Qg(th, pi/2)); - -% Angles -phi = @(r) atan(r/f); -gamma = @(beta,r) acos(cos(pi/2-phi(r)).*cos(beta)); -thi = @(beta,r,y) asin(y/Rsp.*sin(gamma(beta,r))); - -Qgy = @(beta,r,y) Qg_avg(thi(beta,r,y)).*cos(phi(r)); -Qsy = @(beta,r,y) Qs_avg(thi(beta,r,y)).*sin(gamma(beta,r)); - -% Intensity profile -a = 1.0; -w0 = a*r_max; - -%I = @(r) P/(pi*r_max^2); % uniform distribution - -A = (1-exp(-2*r_max.^2/w0^2)); -I0 = 2*P/(pi*w0^2)*A; -I = @(r) I0*exp(-2*r.^2/w0^2); % Gaussian TEM00 beam - -%P0 = exp(0.5)*w0*0.0025/4.81; -%A = 2*pi*integral(@(r) r.*besselj(0,2.405/w0*r).^2,0,r_max)/P0; -%I = @(r) 1/(A*P0)*besselj(0,2.405/w0*r).^2; % Bessel beam - -% Intensity profile graphics -rho = linspace(-r_max, r_max, 500); -figure -plot(rho, I(rho)/max(I(rho)),'k') -grid -xlabel('r, m') -ylabel('I(r)') -sdf('my') - -% Integration -Qres_g = @(y) 1/(pi*r_max^2)*integral2(@(beta,r) r.*I(r).*... - iscomplex(Qgy(beta,r,y)),0,2*pi,0,r_max,... - 'Method','iterated','AbsTol',1e-6,'RelTol',1e-6); - -Qres_s = @(y) 1/(pi*r_max^2)*integral2(@(beta,r) r.*I(r).*... - iscomplex(Qsy(beta,r,y)),0,2*pi,0,r_max,... - 'Method','iterated','AbsTol',1e-6,'RelTol',1e-6); - -% Calulation -N = 150; -y = linspace(-2*Rsp,2*Rsp,N); -Transverse_g = zeros(1,N); -Transverse_s = zeros(1,N); - -for ii = 1:N - Transverse_g(ii) = abs(Qres_g(y(ii))); - Transverse_s(ii) = Qres_s(y(ii)); -end - -Transverse = abs(Transverse_g) + Transverse_s; - -F0 = n1*P/c0; % net force; - -%Graphics -figure -plot(y,F0*Transverse_g,'r--',y,F0*Transverse_s,'b-.',y,F0*Transverse,'k') -legend('F_{g}','F_{s}','F_{t}') -xlabel('r, μ') -ylabel('F, Ν') -grid -sdf('my') diff --git a/Matlab/angles.m b/Matlab/angles.m deleted file mode 100644 index 798852a..0000000 --- a/Matlab/angles.m +++ /dev/null @@ -1,21 +0,0 @@ -load_constants - -function angle = th_r(theta) - angle = asin(n1 / n2 * sin(theta)); -end - -function angle = ph_i(r) - angle = atan(r / f); -end - -function angle = th_i(r, z) - angle = asin(z / Rsp .* sin(ph_i(r))); -end - -function angle = gamma(beta, r) - angle = acos(cos(pi/2 - ph_i(r)) .* cos(beta)); -end - -function angle = theta(beta, r, y) - angle = asin(y / Rsp .* sin(gamma(beta, r))); -end \ No newline at end of file diff --git a/Matlab/axial.asv b/Matlab/axial.asv new file mode 100644 index 0000000..e661f62 --- /dev/null +++ b/Matlab/axial.asv @@ -0,0 +1,70 @@ +clear +close all +clc +format compact + +% 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". +% There are axial forces only + +load_constants + +% Factors plots + +theta = linspace(0, pi/2, 500); +figure +plot(theta, qs_factor(theta, pi/4, n1, n2), ... + theta, -qg_factor(theta, pi/4, n1, n2), ... + theta, qmag_factor(theta, pi/4, n1, n2)) +grid +xlabel('$\theta$, $^{\circ}$', 'Interpreter', 'latex') +ylabel('Q') + +% Intensity profile plots +rho = linspace(-r_max, r_max, 500); +figure +I = gauss(rho, r_max, w0, P); +I0 = max(I); +plot(rho, I/I0, 'k') +grid +xlabel('r, ΠΌ') +ylabel('I(r)') + +% Integration +G0 = gauss_peak(r_max, w0, P); +Qres_g = @(z) 2 * pi * G0 * integral2(@(beta, r) r .* gauss(r, r_max, w0, P) .* ... + iscomplex(qg_z_factor(r, z, n1, n2, Rsp, f)), 0, 2*pi, 0, r_max, ... + 'Method', 'iterated', 'AbsTol', 1e-6, 'RelTol', 1e-6); + +Qres_s = @(z) 1 / (pi * r_max^2) * integral2(@(beta, r) r .* gauss(r, r_max, w0, P) .* ... + iscomplex(qs_z_factor(r, z, n1, n2, Rsp, f)), 0, 2*pi, 0, r_max, ... + 'Method', 'iterated', 'AbsTol', 1e-6, 'RelTol', 1e-6); + +% Calulation +N = 200; +z = linspace(-2*Rsp, 2*Rsp, N); +Axial_g = zeros(1, N); +Axial_s = zeros(1, N); + +wb = waitbar(0, 'Calculating...'); +for ii = 1:N + Axial_g(ii) = Qres_g(z(ii)); + Axial_s(ii) = Qres_s(z(ii)); + waitbar(ii / N, wb, 'Calculating...'); +end +close(wb); + +Axial_g = fliplr(Axial_g); +Axial_s = fliplr(Axial_s); +Axial = Axial_g + Axial_s; +z = -fliplr(z); + +% Plots +figure +plot(z, F0 * Axial_g, 'b-.', ... + z, F0 * Axial_s, 'r--', ... + z, F0 * Axial, 'k') +legend('F_{g}','F_{s}','F_{t}') +xlabel('r, ΠΌ') +ylabel('F, Н') +grid \ No newline at end of file diff --git a/Matlab/axial.m b/Matlab/axial.m index b6cbffe..deca9d6 100644 --- a/Matlab/axial.m +++ b/Matlab/axial.m @@ -8,13 +8,14 @@ format compact % There are axial forces only load_constants -import factors.* % Factors plots theta = linspace(0, pi/2, 500); figure -plot(theta, Qs(theta, pi/4), theta, -Qg(theta, pi/4), theta, Qmag(theta, pi/4)) +plot(theta, qs_factor(theta, pi/4, n1, n2), ... + theta, -qg_factor(theta, pi/4, n1, n2), ... + theta, qmag_factor(theta, pi/4, n1, n2)) grid xlabel('$\theta$, $^{\circ}$', 'Interpreter', 'latex') ylabel('Q') @@ -22,26 +23,28 @@ ylabel('Q') % Intensity profile plots rho = linspace(-r_max, r_max, 500); figure -plot(rho, I(rho)/max(I(rho)),'k') +I = gauss(rho, r_max, w0); +I0 = max(I); +plot(rho, I/I0, 'k') grid xlabel('r, ΠΌ') ylabel('I(r)') % Integration -G0 = 2/(A*pi^2*r_max^2*w0^2); -Qres_g = @(z) G0 * integral2(@(beta,r) r.*I(r).*... - iscomplex(Qgz(r,z)),0,2*pi,0,r_max,... - 'Method','iterated','AbsTol',1e-6,'RelTol',1e-6); +G0 = gauss_peak(r_max, w0); +Qres_g = @(z) 2 * pi * G0 * integral2(@(beta, r) r .* gauss(r, r_max, w0) .* ... + 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) G0 * integral2(@(beta,r) r.*I(r).*... - iscomplex(Qsz(r,z)),0,2*pi,0,r_max,... - 'Method','iterated','AbsTol',1e-6,'RelTol',1e-6); +Qres_s = @(z) 2 * pi * G0 * integral2(@(beta, r) r .* gauss(r, r_max, w0) .* ... + iscomplex(qs_z_factor(r, z, n1, n2, Rsp, f)), 0, 2*pi, 0, r_max, ... + 'Method', 'iterated', 'AbsTol', 1e-12, 'RelTol', 1e-6); % Calulation N = 200; -z = linspace(-2*Rsp,2*Rsp,N); -Axial_g = zeros(1,N); -Axial_s = zeros(1,N); +z = linspace(-2*Rsp, 2*Rsp, N); +Axial_g = zeros(1, N); +Axial_s = zeros(1, N); wb = waitbar(0, 'Calculating...'); for ii = 1:N @@ -58,7 +61,9 @@ z = -fliplr(z); % Plots figure -plot(z,F0*Axial_g,'b-.',z,F0*Axial_s,'r--',z,F0*Axial,'k') +plot(z, F0 * Axial_g, 'b-.', ... + z, F0 * Axial_s, 'r--', ... + z, F0 * Axial, 'k') legend('F_{g}','F_{s}','F_{t}') xlabel('r, ΠΌ') ylabel('F, Н') diff --git a/Matlab/bessel.m b/Matlab/bessel.m new file mode 100644 index 0000000..2561504 --- /dev/null +++ b/Matlab/bessel.m @@ -0,0 +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; +end \ No newline at end of file diff --git a/Matlab/bessel_peak.m b/Matlab/bessel_peak.m new file mode 100644 index 0000000..66c1834 --- /dev/null +++ b/Matlab/bessel_peak.m @@ -0,0 +1,3 @@ +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/coefficients.m b/Matlab/coefficients.m deleted file mode 100644 index a6808e6..0000000 --- a/Matlab/coefficients.m +++ /dev/null @@ -1,9 +0,0 @@ -function reflectivity = R(th, psi) - reflectivity = (tan(th - th_r(th)).^2 ./ tan(th + th_r(th)).^2) .* ... - cos(psi).^2 + ... - (sin(th - th_r(th)).^2 ./ sin(th + th_r(th)).^2) .* sin(psi).^2; -end - -function transparency = T(th, psi) - transparency = 1 - R(th, psi); -end \ No newline at end of file diff --git a/Matlab/factors.m b/Matlab/factors.m deleted file mode 100644 index 4f4ee0a..0000000 --- a/Matlab/factors.m +++ /dev/null @@ -1,37 +0,0 @@ -% Factors -function factor = Qs(th, psi) - factor = 1 + R(th, psi) .* cos(2*th) - T(th, psi).^2.*... - (cos(2*th - 2*th_r(th)) + R(th, psi) .* cos(2*th)) ./ ... - (1 + R(th, psi).^2 + 2*R(th, psi) .* cos(2*th_r(th))); -end - -function factor = Qg(th, psi) - factor = R(th, psi) .* sin(2*th) - T(th, psi).^2 .* ... - (sin(2*th - 2*th_r(th)) + R(th, psi) .* sin(2*th)) ./ ... - (1 + R(th, psi).^2 + 2*R(th, psi) .* cos(2*th_r(th))); -end - -function factor = Qmag(th, psi) - factor = sqrt(Qs(th, psi).^2 + Qg(th, psi).^2); -end - -% Average factors -function factor = Qs_avg(th) - factor = 0.5*(Qs(th, 0) + Qs(th, pi/2)); -end - -function factor = Qg_avg(th) - factor = 0.5*(Qg(th, 0) + Qg(th, pi/2)); -end - -function factor = Qmag_avg(th) - factor = sqrt(Qs_avg(th).^2 + Qg_avg(th).^2); -end - -function factor = Qgz(r, z) - factor = -Qg_avg(th_i(r, z)) .* sin(ph_i(r)); -end - -function factor = Qsz(r, z) - factor = Qs_avg(th_i(r, z)) .* cos(ph_i(r)); -end \ No newline at end of file diff --git a/Matlab/gamma_angle.m b/Matlab/gamma_angle.m new file mode 100644 index 0000000..2337ed3 --- /dev/null +++ b/Matlab/gamma_angle.m @@ -0,0 +1,3 @@ +function angle = gamma_angle(beta, r, focus) + angle = acos(cos(pi/2 - phi_i(r, focus)) .* cos(beta)); +end \ No newline at end of file diff --git a/Matlab/gauss.m b/Matlab/gauss.m new file mode 100644 index 0000000..4e15277 --- /dev/null +++ b/Matlab/gauss.m @@ -0,0 +1,4 @@ +% Gaussian TEM00 beam +function g = gauss(r, r_max, w0) + g = 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 new file mode 100644 index 0000000..4216cd5 --- /dev/null +++ b/Matlab/gauss_peak.m @@ -0,0 +1,4 @@ +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/intensity.m b/Matlab/intensity.m deleted file mode 100644 index 691b684..0000000 --- a/Matlab/intensity.m +++ /dev/null @@ -1,7 +0,0 @@ -load_constants - -function g = gauss(r) - A = (1-exp(-2*r_max.^2 / w0^2)); - I0 = 2*P / (pi * w0^2 * A); - g = I0 * exp(-2 * r.^2 / w0^2); % Gaussian TEM00 beam -end \ No newline at end of file diff --git a/Matlab/load_constants.m b/Matlab/load_constants.m index 49effdd..295920c 100644 --- a/Matlab/load_constants.m +++ b/Matlab/load_constants.m @@ -1,13 +1,14 @@ -Rsp = 1.0e-6; % sphere radius [m] -n1 = 1.33; % index of refraction of the immersion medium -n2 = 1.6; % index of refraction of the fused silica at wavelength 523 nm +Rsp = 1.03e-6; % sphere radius [m] +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; % relative refraction index c0 = 3e8; % speed of light [m/s] NA = 1.25; % numerical apertur e -f = 100.0e-3; % objective lens focus or WD [m] -P = 20.0e-3; % power of the laser [W] +f = 2.0e-3; % objective lens focus or WD [m] +P = 51.7e-3; % power of the laser [W] F0 = n1*P/c0; % resulting force [N] th_max = asin(NA/n1); % maximum incidence angle + +ratio = 1.0; % the ratio of the beam radius to the aperture radius r_max = f * tan(th_max); % radius of a Gaussian beam (1:1 with input aperture condition) -aperture = 1.0; -w0 = aperture * r_max; % Gaussian beam waist radius [m] \ No newline at end of file +w0 = ratio * r_max; % Gaussian beam waist radius [m] \ No newline at end of file diff --git a/Matlab/phi_i.m b/Matlab/phi_i.m new file mode 100644 index 0000000..68b9a36 --- /dev/null +++ b/Matlab/phi_i.m @@ -0,0 +1,3 @@ +function angle = phi_i(r, focus) + angle = atan(r / focus); +end \ No newline at end of file diff --git a/Matlab/qg_avg_factor.m b/Matlab/qg_avg_factor.m new file mode 100644 index 0000000..90ad0eb --- /dev/null +++ b/Matlab/qg_avg_factor.m @@ -0,0 +1,3 @@ +function factor = qg_avg_factor(th, n1, n2) + factor = 0.5*(qg_factor(th, 0, n1, n2) + qg_factor(th, pi/2, n1, n2)); +end \ No newline at end of file diff --git a/Matlab/qg_factor.m b/Matlab/qg_factor.m new file mode 100644 index 0000000..879ed8b --- /dev/null +++ b/Matlab/qg_factor.m @@ -0,0 +1,8 @@ +function factor = qg_factor(th, psi, n1, n2) + R = reflectivity(th, psi, n1, n2); + T = transmittance(th, psi, n1, n2); + theta_refl = th_r(th, n1, n2); + factor = R .* sin(2*th) - T.^2 .* ... + (sin(2*th - 2*theta_refl) + R .* sin(2*th)) ./ ... + (1 + R.^2 + 2*R .* cos(2*theta_refl)); +end \ No newline at end of file diff --git a/Matlab/qg_y_factor.m b/Matlab/qg_y_factor.m new file mode 100644 index 0000000..09effdc --- /dev/null +++ b/Matlab/qg_y_factor.m @@ -0,0 +1,3 @@ +function factor = qg_y_factor(beta, r, y, n1, n2, Rsp, focus) + factor = qg_avg_factor(th_i_y(beta, r, y, Rsp, focus), n1, n2).*sin(gamma_angle(beta, r, focus)); +end \ No newline at end of file diff --git a/Matlab/qg_z_factor.m b/Matlab/qg_z_factor.m new file mode 100644 index 0000000..4ca9c7c --- /dev/null +++ b/Matlab/qg_z_factor.m @@ -0,0 +1,3 @@ +function factor = qg_z_factor(r, z, n1, n2, Rsp, focus) + factor = -qg_avg_factor(th_i_z(r, z, Rsp, focus), n1, n2) .* sin(phi_i(r, focus)); +end \ No newline at end of file diff --git a/Matlab/qmag_avg_factor.m b/Matlab/qmag_avg_factor.m new file mode 100644 index 0000000..4752a05 --- /dev/null +++ b/Matlab/qmag_avg_factor.m @@ -0,0 +1,3 @@ +function factor = qmag_avg_factor(th, n1, n2) + factor = sqrt(qs_avg_factor(th, n1, n2).^2 + qg_avg_factor(th, n1, n2).^2); +end \ No newline at end of file diff --git a/Matlab/qmag_factor.m b/Matlab/qmag_factor.m new file mode 100644 index 0000000..d2ee7b1 --- /dev/null +++ b/Matlab/qmag_factor.m @@ -0,0 +1,3 @@ +function factor = qmag_factor(th, psi, n1, n2) + factor = sqrt(qs_factor(th, psi, n1, n2).^2 + qg_factor(th, psi, n1, n2).^2); +end \ No newline at end of file diff --git a/Matlab/qs_avg_factor.m b/Matlab/qs_avg_factor.m new file mode 100644 index 0000000..67ce068 --- /dev/null +++ b/Matlab/qs_avg_factor.m @@ -0,0 +1,3 @@ +function factor = qs_avg_factor(th, n1, n2) + factor = 0.5*(qs_factor(th, 0, n1, n2) + qs_factor(th, pi/2, n1, n2)); +end \ No newline at end of file diff --git a/Matlab/qs_factor.m b/Matlab/qs_factor.m new file mode 100644 index 0000000..a7ab290 --- /dev/null +++ b/Matlab/qs_factor.m @@ -0,0 +1,9 @@ +function factor = qs_factor(th, psi, n1, n2) + R = reflectivity(th, psi, n1, n2); + T = transmittance(th, psi, n1, n2); + theta_refl = th_r(th, n1, n2); + + factor = 1 + R .* cos(2*th) - T.^2.*... + (cos(2*th - 2*theta_refl) + R .* cos(2*th)) ./ ... + (1 + R.^2 + 2*R .* cos(2*theta_refl)); +end \ No newline at end of file diff --git a/Matlab/qs_y_factor.m b/Matlab/qs_y_factor.m new file mode 100644 index 0000000..8315bfc --- /dev/null +++ b/Matlab/qs_y_factor.m @@ -0,0 +1,3 @@ +function factor = qs_y_factor(beta, r, y, n1, n2, Rsp, focus) + factor = qs_avg_factor(th_i_y(beta, r, y, Rsp, focus), n1, n2).*cos(phi_i(r, focus)); +end \ No newline at end of file diff --git a/Matlab/qs_z_factor.m b/Matlab/qs_z_factor.m new file mode 100644 index 0000000..d04adeb --- /dev/null +++ b/Matlab/qs_z_factor.m @@ -0,0 +1,3 @@ +function factor = qs_z_factor(r, z, n1, n2, Rsp, focus) + factor = qs_avg_factor(th_i_z(r, z, Rsp, focus), n1, n2) .* cos(phi_i(r, focus)); +end \ No newline at end of file diff --git a/Matlab/reflectivity.m b/Matlab/reflectivity.m new file mode 100644 index 0000000..10d5978 --- /dev/null +++ b/Matlab/reflectivity.m @@ -0,0 +1,6 @@ +function R = reflectivity(th, psi, n1, n2) + theta_refl = th_r(th, n1, n2); + R = (tan(th - theta_refl).^2 ./ tan(th + theta_refl).^2) .* ... + cos(psi).^2 + ... + (sin(th - theta_refl).^2 ./ sin(th + theta_refl).^2) .* sin(psi).^2; +end \ No newline at end of file diff --git a/Matlab/th_i_y.m b/Matlab/th_i_y.m new file mode 100644 index 0000000..438cdcd --- /dev/null +++ b/Matlab/th_i_y.m @@ -0,0 +1,3 @@ +function angle = th_i_y(beta, r, y, Rsp, focus) + angle = asin(y / Rsp .* sin(gamma_angle(beta, r, focus))); +end \ No newline at end of file diff --git a/Matlab/th_i_z.m b/Matlab/th_i_z.m new file mode 100644 index 0000000..11d29b8 --- /dev/null +++ b/Matlab/th_i_z.m @@ -0,0 +1,3 @@ +function angle = th_i_z(r, z, Rsp, focus) + angle = asin(z / Rsp .* sin(phi_i(r, focus))); +end \ No newline at end of file diff --git a/Matlab/th_r.m b/Matlab/th_r.m new file mode 100644 index 0000000..e7aa0ee --- /dev/null +++ b/Matlab/th_r.m @@ -0,0 +1,3 @@ +function angle = th_r(theta, n1, n2) + angle = asin(n1 / n2 * sin(theta)); +end \ No newline at end of file diff --git a/Matlab/theta_angle.m b/Matlab/theta_angle.m new file mode 100644 index 0000000..acf2b75 --- /dev/null +++ b/Matlab/theta_angle.m @@ -0,0 +1,3 @@ +function angle = theta_angle(beta, r, y, Rsp, focus) + angle = asin(y / Rsp .* sin(gamma(beta, r, focus))); +end \ No newline at end of file diff --git a/Matlab/transmittance.m b/Matlab/transmittance.m new file mode 100644 index 0000000..1b1558e --- /dev/null +++ b/Matlab/transmittance.m @@ -0,0 +1,3 @@ +function T = transmittance(th, psi, n1, n2) + T = 1 - reflectivity(th, psi, n1, n2); +end \ No newline at end of file diff --git a/Matlab/transverse.m b/Matlab/transverse.m new file mode 100644 index 0000000..bc278d0 --- /dev/null +++ b/Matlab/transverse.m @@ -0,0 +1,56 @@ +%all distances in m +close all +clear +clc +format compact + +% 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". +% There are transverse forces only + +load_constants + +% Intensity profile graphics +rho = linspace(-r_max, r_max, 500); +I = gauss(rho, r_max, w0); +I0 = max(I); +figure +plot(rho, I/I0, 'k') +grid +xlabel('r, m') +ylabel('I(r)') + +% Integration +Qres_g = @(y) integral2(@(beta, r) r .* gauss(r, r_max, w0) .* ... + iscomplex(qg_y_factor(beta, r, y, n1, n2, Rsp, f)), 0, 2*pi, 0, r_max, ... + 'Method', 'iterated', 'AbsTol', 1e-6, 'RelTol', 1e-6); + +Qres_s = @(y) integral2(@(beta, r) r .* gauss(r, r_max, w0) .* ... + iscomplex(qs_y_factor(beta, r, y, n1, n2, Rsp, f)), 0, 2*pi, 0, r_max, ... + 'Method', 'iterated', 'AbsTol', 1e-6, 'RelTol', 1e-6); + +% Calulation +N = 150; +y = linspace(-2*Rsp, 2*Rsp, N); +Transverse_g = zeros(1, N); +Transverse_s = zeros(1, N); + +wb = waitbar(0, 'Calculating...'); +for ii = 1:N + Transverse_g(ii) = abs(Qres_g(y(ii))); + Transverse_s(ii) = Qres_s(y(ii)); + waitbar(ii / N, wb, 'Calculating...'); +end +close(wb); + +Transverse = abs(Transverse_g) + Transverse_s; + +%Graphics +figure +plot(y, F0 * Transverse_g,'r--', ... + y, F0 * Transverse_s, 'b-.', ... + y, F0*Transverse, 'k') +legend('F_{g}','F_{s}','F_{t}') +xlabel('r, μ') +ylabel('F, Ν') +grid diff --git a/Matlab/uniform.m b/Matlab/uniform.m new file mode 100644 index 0000000..7a10840 --- /dev/null +++ b/Matlab/uniform.m @@ -0,0 +1,4 @@ +% Beam with uniform distribution +function u = uniform(r, w0, r_max, P) + u = P / (pi*r_max^2); +end \ No newline at end of file