Fixed multipliers, refactored ray efficiencies

Matlab/Forces_Ashkin_Ray_Efficiencies.m

@@ -1,48 +0,0 @@
-close all
-clear
-format compact
-clc
-
-% 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
-a = 1.0e-6; % radius of the bead
-n1 = 1.0; % index of rafraction of the medium
-n = 1.4607; % n2/n1
-n2 = n*n1; % index of refraction of the fused silica
-c0 = 3e8; % speed of light
-
-%reflectivity
-R = @(th,psi) (tan(th-asin(n1/n2*sin(th))).^2./...
-    tan(th+asin(n1/n2*sin(th))).^2).*cos(psi).^2+...
-    (sin(th-asin(n1/n2*sin(th))).^2./...
-    sin(th+asin(n1/n2*sin(th))).^2).*sin(psi).^2;
-
-%transparency
-T = @(th,psi) 1-R(th,psi);
-
-r = @(th) asin(n1/n2*sin(th));
-
-% Factors
-Qs = @(th, psi) 1 + R(th, psi) .* cos(2*th) - T(th,psi).^2 .* (cos(2*th -...
-    2*r(th)) + R(th, psi) .* cos(2*th)) ./ (1 + R(th,psi).^2 +...
-    2*R(th,psi) .* cos(2*r(th)));
-
-Qg = @(th, psi) R(th, psi) .* sin(2*th) - T(th,psi).^2 .* (sin(2*th -...
-    2*r(th)) + R(th, psi) .* sin(2*th)) ./ (1 + R(th,psi).^2 +...
-    2*R(th,psi) .* cos(2*r(th)));
-
-Qmag = @(th, psi) sqrt(Qs(th, psi).^2 + Qg(th, psi).^2);
-
-t = linspace(0, pi/2, 1000);
-t_deg = t*180/pi;
-pol = pi/4;
-
-figure
-plot(t_deg, Qs(t, pol),'r--', t_deg, -Qg(t, pol),'b-.', t_deg, Qmag(t, pol),'k');
-grid
-xlabel('\theta, deg')
-ylabel('Q')
-legend('Q_s','Q_g','Q_t','location','northwest')
-sdf('my')

Matlab/axial.asv

@@ -1,70 +0,0 @@
-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

Matlab/axial.m

@@ -9,17 +9,6 @@ format compact
 
 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

Matlab/ray_efficiencies.m

@@ -0,0 +1,23 @@
+close all
+clear
+format compact
+clc
+
+% 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
+
+load_constants
+
+t = linspace(0, pi/2, 1000);
+t_deg = t * 180/pi;
+psi = pi / 4;
+
+figure
+plot(t_deg, qs_factor(t, psi, n1, n2), ...
+    t_deg, -qg_factor(t, psi, n1, n2), ...
+    t_deg, qmag_factor(t, psi, n1, n2));
+grid
+title(['\psi = ', num2str(psi * 180/pi), '^\circ'])
+xlabel('\theta, ^\circ')
+ylabel('Q')
+legend('Q_s','Q_g','Q_t','location','northwest')

Matlab/transverse.m

@@ -21,13 +21,14 @@ xlabel('r, m')
 ylabel('I(r)')
 
 % Integration
-Qres_g = @(y) integral2(@(beta, r) r .* gauss(r, r_max, w0) .* ...
+G0 = gauss_peak(r_max, w0);
+Qres_g = @(y) G0 * 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);
+    'Method', 'iterated', 'AbsTol', 1e-8, 'RelTol', 1e-6);
 
-Qres_s = @(y) integral2(@(beta, r) r .* gauss(r, r_max, w0) .* ...
+Qres_s = @(y) G0 * 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);
+    'Method', 'iterated', 'AbsTol', 1e-8, 'RelTol', 1e-6);
 
 % Calulation
 N = 150;
@@ -43,7 +44,7 @@ for ii = 1:N
 end
 close(wb);
 
-Transverse = abs(Transverse_g) + Transverse_s;
+Transverse = Transverse_g + Transverse_s;
 
 %Graphics
 figure