From 8958bf5204d9354ccff578ce64918faa7c2c0518 Mon Sep 17 00:00:00 2001 From: BrokenVoodooDoll Date: Sun, 19 Mar 2023 13:45:34 +0400 Subject: [PATCH] Split up files --- Matlab/Forces_Ashkin_Axial.m | 112 ----------------------------------- Matlab/angles.m | 21 +++++++ Matlab/axial.m | 65 ++++++++++++++++++++ Matlab/coefficients.m | 9 +++ Matlab/factors.m | 37 ++++++++++++ Matlab/intensity.m | 7 +++ Matlab/iscomplex.m | 10 +--- Matlab/load_constants.m | 13 ++++ 8 files changed, 153 insertions(+), 121 deletions(-) delete mode 100644 Matlab/Forces_Ashkin_Axial.m create mode 100644 Matlab/angles.m create mode 100644 Matlab/axial.m create mode 100644 Matlab/coefficients.m create mode 100644 Matlab/factors.m create mode 100644 Matlab/intensity.m create mode 100644 Matlab/load_constants.m diff --git a/Matlab/Forces_Ashkin_Axial.m b/Matlab/Forces_Ashkin_Axial.m deleted file mode 100644 index 3e1bf21..0000000 --- a/Matlab/Forces_Ashkin_Axial.m +++ /dev/null @@ -1,112 +0,0 @@ -%all distances in m -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 - -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))); - -Qmag = @(th, psi) sqrt(Qs(th, psi).^2 + Qg(th, psi).^2); - -% 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)); -Qmag_avg = @(th) sqrt(Qs_avg(th).^2 + Qg_avg(th).^2); - -% Angles -phi = @(r) atan(r/f); -thi = @(r,z) asin(z/Rsp.*sin(phi(r))); - -Qgz = @(r,z) -Qg_avg(thi(r,z)).*sin(phi(r)); -Qsz = @(r,z) Qs_avg(thi(r,z)).*cos(phi(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 = P*2/(pi*w0^2)/A; -I = @(r) I0*exp(-2*r.^2/w0^2); % Gaussian TEM00 beam - -%A = 2*pi*integral(@(r) r.*besselj(0,2.405/w0*r).^2,0,r_max)/P0; -%w0_bb = 0.5*r_max; -%I0 = P*2/(pi*w0^2); -%I = @(r) I0*besselj(0,2.405/w0_bb*r).^2.*exp(-2*r.^2/w0^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, μ') -ylabel('I(r)') -sdf('my') - -% Integration -Qres_g = @(z) 1/(pi*r_max^2)*2/(A*pi*w0^2)*integral2(@(beta,r) r.*I(r).*... - iscomplex(Qgz(r,z)),0,2*pi,0,r_max,... - 'Method','iterated','AbsTol',1e-12,'RelTol',1e-6); - -Qres_s = @(z) 1/(pi*r_max^2)*2/(A*pi*w0^2)*integral2(@(beta,r) r.*I(r).*... - iscomplex(Qsz(r,z)),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); - -for ii = 1:N - Axial_g(ii) = Qres_g(z(ii)); - Axial_s(ii) = Qres_s(z(ii)); -end - -F0 = n1*P/c0; % net force - -Axial_g = fliplr(Axial_g); -Axial_s = fliplr(Axial_s); -Axial = Axial_g + Axial_s; -z = -fliplr(z); - -%Graphics -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 -sdf('my') \ No newline at end of file diff --git a/Matlab/angles.m b/Matlab/angles.m new file mode 100644 index 0000000..798852a --- /dev/null +++ b/Matlab/angles.m @@ -0,0 +1,21 @@ +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.m b/Matlab/axial.m new file mode 100644 index 0000000..b6cbffe --- /dev/null +++ b/Matlab/axial.m @@ -0,0 +1,65 @@ +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 +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)) +grid +xlabel('$\theta$, $^{\circ}$', 'Interpreter', 'latex') +ylabel('Q') + +% Intensity profile plots +rho = linspace(-r_max, r_max, 500); +figure +plot(rho, I(rho)/max(I(rho)),'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); + +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); + +% 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/coefficients.m b/Matlab/coefficients.m new file mode 100644 index 0000000..a6808e6 --- /dev/null +++ b/Matlab/coefficients.m @@ -0,0 +1,9 @@ +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 new file mode 100644 index 0000000..4f4ee0a --- /dev/null +++ b/Matlab/factors.m @@ -0,0 +1,37 @@ +% 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/intensity.m b/Matlab/intensity.m new file mode 100644 index 0000000..691b684 --- /dev/null +++ b/Matlab/intensity.m @@ -0,0 +1,7 @@ +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/iscomplex.m b/Matlab/iscomplex.m index 152359f..6f1f862 100644 --- a/Matlab/iscomplex.m +++ b/Matlab/iscomplex.m @@ -1,11 +1,3 @@ function S = iscomplex(A) -[n,m] = size(A); -S = zeros(n,m); - -for ii = 1:n - for jj = 1:m - S(ii,jj) = A(ii,jj)*double(isreal(A(ii,jj))); - end -end - + S = A .* double(imag(A) == 0); end \ No newline at end of file diff --git a/Matlab/load_constants.m b/Matlab/load_constants.m new file mode 100644 index 0000000..49effdd --- /dev/null +++ b/Matlab/load_constants.m @@ -0,0 +1,13 @@ +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 +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] +F0 = n1*P/c0; % resulting force [N] +th_max = asin(NA/n1); % maximum incidence angle +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