Started work on the 1.4.0 examples

examples/axial_equilibrium.m

@@ -1,50 +1,47 @@
-% Example file for finding spring constants
+% ott.axial_equilibrium can be used to calculate the trapping position
+% and spring constant along the z axis.  If the beam is rotated, the
+% method can also be used for other axies too.
 %
 % This file is part of the optical tweezers toolbox.
 % See LICENSE.md for information about using/distributing this file.
 
-import ott.*
-import ott.utils.*
-
 % Make warnings less obtrusive
 ott_warning('once');
 change_warnings('off');
 
 % Specify refractive indices
 n_medium = 1.34;
 n_particle = 1.59;
-n_relative = n_particle/n_medium;
 
-% If you want to give all measurements in wavelengths in the surrounding
-% medium, then:
-wavelength = 1;
-% wavelength = wavelength0 / n_medium;
-% else you can give it in any units you want. Only k times lengths matters
-k = 2*pi/wavelength;
+% Specify the wavelength in freespace [m]
+wavelength = 1064.0e-9;
+
+% Specify the particle radius (sphere)
+radius = 1.0*wavelength/n_medium;
 
-radius = 1.5;
-Nmax = ka2nmax(k*radius);
+%% Calculate the beam
 
-if Nmax < 12
-    Nmax = 12;
-end
+beam = ott.BscPmGauss('angle_deg', 50, ...
+    'polarisation', [ 1 0 ], 'power', 1.0);
 
-% a Gaussian beam: w0 = 2/(k*tan(theta))
-beam_angle = 50; % Convergence half-angle of 50 degrees
+%% Calculate the particle T-matrix
 
-% Polarisation. [ 1 0 ] is plane-polarised along the x-axis, [ 0 1 ] is
-% y-polarised, and [ 1 -i ] and [ 1 i ] are circularly polarised.
-polarisation = [ 1 0 ];
+T = ott.Tmatrix.simple('sphere', radius, ...
+    'n_medium', n_medium, ...
+    'n_particle', n_particle, ...
+    'wavelength0', wavelength);
 
-[a,b] = bsc_pointmatch_farfield(Nmax,1,[ 0 0 beam_angle 1 polarisation 90 ]);
+%% Find the equilibrium and trap stiffness in the x and x directions
 
-% If you're going to do a range of particles, then the T-matrix has to
-% calculated inside the loop.
+% Find the equilibrium in the z-direction
+[z,kz] = ott.axial_equilibrium(T, beam)
 
-% To search for a refractive index, I recommend a bisection search. For an
-% example of bisection search, see find_axial_equilibrium.m
+% Translate the beam to the z-axis equilibrium
+beam = beam.translateZ(z);
 
-T = tmatrix_mie(Nmax,k,k*n_relative,radius);
+% Rotate the beam about the y axis (so the beam is aligned with the x axis)
+beam = beam.rotateY(pi/2.0);
 
-[z,k] = axial_equilibrium(T,a,b)
+% Calculate the equilibrium in the x-direction
+[x,kx] = ott.axial_equilibrium(T, beam);
 

examples/beam_visualisation.m

@@ -1,8 +1,140 @@
-%Example farfield plotting code:
+% Demonstrates generation of near and far field visualisations of the beams
+%
+% Shows the radiance at the focal plane, transverse to the focal
+% plane and in the far field.
+%
+% Multiple beam configurations are included, use the beam_type
+% variable to vary the beam that is being visualised.
+%
+% This file is part of the optical tweezers toolbox.
+% See LICENSE.md for information about using/distributing this file.
 
-%build beam:
-[a1,b1]=bsc_pointmatch_farfield(20,1,[ 1 1 lg_mode_w0([1,1],65) 1 1 -1i 90 0 0 0 ]);
-[n,m]=combined_index(find(abs(a1)|abs(b1)));
+%% Setup the beam
+
+% Select a beam type to visualise (gaussian, lg, hg, ig, addition, scattered)
+beam_type = 'gaussian';
+
+switch beam_type
+  case 'gaussian'
+
+    % Create a simple Gaussian beam with circular polarisation
+    beam = ott.BscPmGauss('NA', 1.02, 'polarisation', [ 1 1i ]);
+
+  case 'lg'
+
+    % Create a simple LG03 beam with circular polarisation
+    beam = ott.BscPmGauss('lg', [ 0 3 ], ...
+        'NA', 1.02, 'polarisation', [ 1 1i ]);
+
+  case 'hg'
+
+    % Create a HG23 beam with circular polarisation
+    %
+    % The beam waist calculation doesn't yet handle HG beams so
+    % we estimate the beam waist using a magic formula
+    beam_angle = asin(NA/n_medium);
+    w0=1.0/pi/tan(abs(theta/180*pi));
+    beam = ott.BscPmGauss('hg', [ 2 3 ], ...
+        'polarisation', polarisation, 'w0', w0);
+
+  case 'ig'
+
+    % Create a IG beam, not sure what would look good yet
+    %
+    % The beam waist calculation doesn't yet handle IG beams so
+    % we estimate the beam waist using a magic formula
+    beam_angle = asin(NA/n_medium);
+    w0=1.0/pi/tan(abs(theta/180*pi));
+    beam = ott.BscPmGauss('ig', [ 2 3 1 2 ], ...
+        'polarisation', polarisation, 'w0', w0);
+
+  case 'addition'
+
+    % Separation between beams
+    dx = 1.0;
+
+    % Create a simple Gaussian beam with circular polarisation
+    beam = ott.BscPmGauss('NA', 1.02, 'polarisation', [ 1 1i ]);
+    beam.Nmax = beam.Nmax + ott.utils.ka2nmax(dx);
+
+    % Displace the beams by +/- dx/2
+    beam1 = beam.translateXyz([dx/2, 0, 0]);
+    beam2 = beam.translateXyz([-dx/2, 0, 0]);
+
+    % Join the beams to create our beam (adding a phase shift to the 2nd beam)
+    beam = beam1 + beam2 * exp(2 * pi * 1i * 1/2);
+
+  case 'scattered'
+
+    % Create a simple Gaussian beam with circular polarisation
+    ibeam = ott.BscPmGauss('NA', 1.02, 'polarisation', [ 1 1i ]);
+
+    % Create a spherical particle to scatter the beam
+    tmatrix = ott.Tmatrix.simple('sphere', 1.0, 'wavelength0', 1.0, ...
+        'n_medium', 1.0, 'n_particle', 1.2);
+
+    % Scatter the beam to create the final beam
+    beam = tmatrix * ibeam;
+
+  otherwise
+    error('Unsupported beam type specified');
+end
+
+%% Generate plot of the intensity at the focal plane
+
+% Create the grid of points we want to view
+nx = 40;
+ny = 40;
+xrange = linspace(-8, 8, nx);
+yrange = linspace(-8, 8, ny);
+[xx, yy] = meshgrid(xrange, yrange);
+xyz = [xx(:) yy(:) zeros(size(xx(:)))].';
+
+% Calculate the E and H fields
+[E, H] = beam.emFieldXyz(xyz);
+
+% Calculate the intensity of the E field
+Ei=reshape(sqrt(sum(real(E).^2,1)),[nx,ny]);
+
+% Calculate the radiance
+I=reshape(sum(abs(E).^2,1),[nx,ny]);
+
+figure(1);
+subplot(1, 2, 1);
+imagesc(xrange, yrange, Ei);
+title('E field intensity (focal plane)');
+subplot(1, 2, 2);
+imagesc(xrange, yrange, I);
+title('radiance (focal plane)');
+
+%% Generate plot of the intensity at the focal plane (same as above)
+
+% Create the grid of points we want to view
+nx = 40;
+ny = 40;
+xrange = linspace(-8, 8, nx);
+yrange = linspace(-8, 8, ny);
+[xx, yy] = meshgrid(xrange, yrange);
+xyz = [xx(:) zeros(size(xx(:))) yy(:)].';
+
+% Calculate the E and H fields
+[E, H] = beam.emFieldXyz(xyz);
+
+% Calculate the intensity of the E field
+Ei=reshape(sqrt(sum(real(E).^2,1)),[nx,ny]);
+
+% Calculate the radiance
+I=reshape(sum(abs(E).^2,1),[nx,ny]);
+
+figure(2);
+subplot(1, 2, 1);
+imagesc(xrange, yrange, Ei);
+title('E field intensity (cross-section)');
+subplot(1, 2, 2);
+imagesc(xrange, yrange, I);
+title('radiance (cross-section)');
+
+%% Generate a figure showing the farfield
 
 %build grid:
 nt=80;
@@ -12,26 +144,26 @@
 [~,theta,phi]=xyz2rtp(x,y,z);
 
 %find far-field in theta, phi:
-[E,H]=farfield(n,m,a1,b1,0*a1,0*b1,theta(:),phi(:));
-
-
-%% find radiant flux:
-I=reshape(sum(abs(E).^2,2),[nt+1,nt+1]);
-
-%% find circular polarisation: 
-% horizontal and vertical require a coordinate transformation:
-% (theta,phi)->(x,y)
-R=reshape(abs(E(:,2)+1i*E(:,3)).^2,[nt+1,nt+1]); %can't remember which is left and right
-L=reshape(abs(E(:,2)-1i*E(:,3)).^2,[nt+1,nt+1]);
-
-%%
-figure(1)
-surf(x,y,z,R,'facecolor','interp','edgecolor','none')
-axis equal
-view([0,270])
-title('right')
-figure(2)
-surf(x,y,z,L,'facecolor','interp','edgecolor','none')
-axis equal
-view([0,270])
-title('left')
+[E,H]=beam.farfield(theta(:),phi(:));
+
+% Calculate the intensity of the E field
+Ei=reshape(sqrt(sum(real(E).^2,1)),[nt+1,nt+1]);
+
+% Calculate the average phase of the E field
+% This should be similar to the pattern we put on a SLM
+Ep=reshape(sum(arg(E),1)/2.0),[nt+1,nt+1]);
+
+% Calculate the radiance
+I=reshape(sum(abs(E).^2,1),[nt+1,nt+1]);
+
+figure(3);
+subplot(1, 2, 1);
+surf(x,y,z,Ei,'facecolor','interp','edgecolor','none')
+title('E field intensity (farfield)');
+subplot(1, 2, 2);
+surf(x,y,z,Ep,'facecolor','interp','edgecolor','none')
+title('E field phase (farfield)');
+subplot(1, 2, 3);
+surf(x,y,z,I,'facecolor','interp','edgecolor','none')
+title('radiance (farfield)');
+