forked from ilent2/ott
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathslm_to_focalplane.m
153 lines (119 loc) · 4.27 KB
/
slm_to_focalplane.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
% Calculates how a SLM pattern would look at the focal plane of a microscope
%
% This file is part of the optical tweezers toolbox.
% See LICENSE.md for information about using/distributing this file.
% Add the toolbox path
addpath('../');
% Close open figures
close all;
% Turn off change warnings for examples
ott.change_warnings('off');
% Numerical aperture of microscope
NA = 1.02;
% Refractive index of medium
index_medium = 1.33;
% Polarisation of incident beam
polarisation = [ 0 1 ];
% Range [in number of wavelengths in medium] we want to visualise fields
% near focal plane
range = 3.0;
% Describes how well matched the beam is to the SLM, to underfill
% changes this value to less than 1.
fillfactor = 1.0;
%% Setup the incident beam and SLM pattern
% Define a grid representing pixels on the SLM
x = linspace(-1, 1, 512);
y = linspace(-1, 1, 512);
x_idx = 1:length(x);
y_idx = 1:length(y);
[xx, yy] = meshgrid(x, y);
[xx_idx, yy_idx] = meshgrid(x_idx, y_idx);
% Calculate the incident beam amplitude (Gaussian beam, LG00)
incident_mode = ott.utils.lgmode(0,0, ...
sqrt(xx.^2+yy.^2)/fillfactor, atan2(yy,xx));
% Generate the phase pattern
slm_phase = zeros(size(incident_mode));
% Set the entire pattern to checkerboard (zero amplitude)
slm_phase(logical(mod(xx_idx + yy_idx, 2))) = pi;
% Unmask the region we want our pattern on
xxyy = ott.utils.rotz(45) * [xx(:) yy(:) zeros(size(xx(:)))].';
slm_phase(sinc(sum([1 0 0] * xxyy, 1)*4) > 0) = pi/2;
%% Show the incident mode and the SLM pattern
figure(1);
imagesc(incident_mode);
title('Incident mode amplitude');
xlabel('X [pixels]'); ylabel('Y [pixels]');
axis('image');
figure(2);
imagesc(slm_phase);
title('SLM phase pattern');
xlabel('X [pixels]'); ylabel('Y [pixels]');
axis('image');
%% Calculate the beam shape coefficients for the paraxial beam
% For the visualisation we need to specify the range where we
% want to plot the fields, so we manually calculate Nmax
Nmax = ott.utils.ka2nmax(2*pi*index_medium*range/2.0);
% Combined the phase and incident beam and make a beam object
combined_beam = exp(1i*slm_phase) .* incident_mode;
beam = ott.BscPmParaxial(-NA, combined_beam, 'Nmax', Nmax, ...
'polarisation', polarisation, 'index_medium', index_medium);
%% Visualise the near fields (transverse to beam axis)
% Create a grid of point we want to view
nx = 80;
ny = 80;
xrange = linspace(-0.5, 0.5, nx)*range/index_medium;
yrange = linspace(-0.5, 0.5, ny)*range/index_medium;
[xx, yy] = meshgrid(xrange, yrange);
xyz = [xx(:) yy(:) zeros(size(xx(:)))].';
% Calculate the electric field
E = beam.emFieldXyz(xyz);
% Calculate the intensity of the E field
Ei=reshape(sqrt(sum(real(E).^2,1)),[nx,ny]);
% Calculate the phase at the focal plane
Ep = reshape(angle(E(3, :)),[nx,ny]);
% Calculate the radiance
I=reshape(sum(abs(E).^2,1),[nx,ny]);
figure(3);
subplot(1, 3, 1);
imagesc(xrange, yrange, Ei);
axis('image'); xlabel('X [\lambda_m]'); ylabel('Y [\lambda_m]');
title('E field intensity');
subplot(1, 3, 2);
imagesc(xrange, yrange, Ep);
axis('image'); xlabel('X [\lambda_m]'); ylabel('Y [\lambda_m]');
title('E field phase');
subplot(1, 3, 3);
imagesc(xrange, yrange, I);
axis('image'); xlabel('X [\lambda_m]'); ylabel('Y [\lambda_m]');
title('radiance');
set(gcf, 'Name','At focal plane','NumberTitle','off');
%% Visualise the near fields (parallel to beam axis)
% Create a grid of point we want to view
nx = 80;
ny = 80;
xrange = linspace(-0.5, 0.5, nx)*range/index_medium;
yrange = linspace(-0.5, 0.5, ny)*range/index_medium;
[xx, yy] = meshgrid(xrange, yrange);
xyz = [xx(:) zeros(size(xx(:))) yy(:)].';
% Calculate the electric field
E = beam.emFieldXyz(xyz);
% Calculate the intensity of the E field
Ei=reshape(sqrt(sum(real(E).^2,1)),[nx,ny]);
% Calculate the phase at the focal plane
Ep = reshape(angle(E(3, :)),[nx,ny]);
% Calculate the radiance
I=reshape(sum(abs(E).^2,1),[nx,ny]);
figure(4);
subplot(1, 3, 1);
imagesc(xrange, yrange, Ei);
axis('image'); xlabel('X [\lambda_m]'); ylabel('Z [\lambda_m]');
title('E field intensity');
subplot(1, 3, 2);
imagesc(xrange, yrange, Ep);
axis('image'); xlabel('X [\lambda_m]'); ylabel('Z [\lambda_m]');
title('E field phase');
subplot(1, 3, 3);
imagesc(xrange, yrange, I);
axis('image'); xlabel('X [\lambda_m]'); ylabel('Z [\lambda_m]');
title('radiance');
set(gcf, 'Name','Along beam axis','NumberTitle','off');