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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
|
% FOURIER OPTICS DEMONSTRATION
%
% Calculates the field intensity on a screen due to diffraction in the Fresnel
% approximation using the Fourier transform method. The field u(x,y) in the source plane
% is a highly concave sperical phase front (ROC=5 mm) incident on a pair of rectangular
% apertures. The field u'(x',y') in the "field plane" is calculated via Fourier transform.
%
% You can change the parameters below (e.g. roc, w, a, b, d) to investigate the effect of
% changing the slit positions and widths.
% ---------------------------------------------------------------------------------------------------------------------
%% The first part of the code is intended to be modified by the user to change the physical and
% computational parameters
% Physical Parameters
% -------------------
c = 3e8; % speed of light in meters per second
epsilon0 = 8.854e-12; % vacuum permittivity in Farads per meter (SI units)
lambda = 633e-9; % optical wavelength in meters
% Source plane
% ------------
xmax=0.002; ymax=xmax; % size of source plane = [-xmax,xmax]*[symax,symax] meters
Nx = 2^nextpow2(512); Ny = 2^nextpow2(Nx); % number of points in source plane grid = Nx*Ny
dx = 2*xmax/(Nx-1); dy=2*ymax/(Ny-1); % interpixel distances (sampling intervals) in the source plane (meters)
x = repmat( ((0:Nx-1)-floor(Nx/2)) *dx, Ny,1); % source plane points (in scaled distance) at which to evaluate the source field
y = repmat( ((0:Ny-1)-floor(Ny/2)).'*dy, 1,Nx); % source plane points (in scaled distance) at which to evaluate the source field
% ABCD Matrix Components
% ----------------------
% Here assumed to be for a free-space/lens/free-space system (f=Inf corresponds to no lens)
L1 = 1e-3;
L2 = 999e-3;
f = Inf;
%f = -10e-3; % Uncomment this line to see an example of Fresnel (near-field) diffraction
M = [[1 L2];[0 1]] * [[1 0];[-1/f 1]] * [[1 L1];[0 1]];
disp('M = ');
dispmat(M);
AA = M(1,1); BB = M(1,2); CC = M(2,1); DD = M(2,2);
% Aperture
% --------
% Field amplitude is non-zero at these values of x, y (i.e. where it passes through the aperture))
% The apertures are defined as logical matrixes that are used to index the source field
% distribution, i.e. Usource(~aperture)=0; UIsource(aperture)= <something nonzero>.
a = 100*1e-6;
wideslit = abs(x)<a & abs(y)<8*a;
a = 25*1e-6;
narrowslit = abs(x)<a & abs(y)<16*a;
a=20*1e-6; % slit width a and interslit distance d (center to center, in meters)
d=8*a; % f=Inf, roc=Inf; w=Inf; gives classic Fraunhofer fringes.
sixslits = (abs(x+3*d)<a | abs(x+2*d)<a | abs(x+d)<a | abs(x)<a | abs(x-d)<a | abs(x-2*d)<a | abs(x-3*d)<a) & abs(y)<16*a;
a = 200*1e-6; % length of side of equilateral triangle (in meters)
triangle = (y<sqrt(3)*x+a/2/sqrt(3)) & (y<-sqrt(3)*x+a/2/sqrt(3))& (y>-a/2/sqrt(3));
a=150e-6;
circle = x.^2+y.^2 < a^2;
aperture = narrowslit; % Choose one of singleslit, sixslits, circle, or triangle
% Source Field
% ------------
% Here, the field is assumed to be due to a Gaussian beam of width "w" and phasefront radius of curv. "roc" incident on the aperture.
roc = 0.5; % radius of curvature of phasefront at source plane, in meters
w = 500e-6; % Gaussian beam width (meters)
I0 = 1e6; % (cycle averaged) maximum intensity of the light in source plane in watts per meter^2
E0 = sqrt(2*I0/c/epsilon0); % field amplitude in the source plane in Newtons/Coulomb.
k=2*pi/lambda; %
r=sqrt(x.^2+y.^2); % radius from center of aperture
Upreap = E0*exp(-r.^2/w^2).*exp(1i*k*r.^2/2/roc)... % field amplitude of beam incident on the source plane aperture
.* exp(1i*0*pi/4*randn(size(x))); % (diverging beam if roc is positive). Possible incoherent phase (this line).
Usource = Upreap; Usource(~aperture)=0; % source plane field amplitude will be zero everywhere except in the aperture
Isource = epsilon0*c/2*abs(Usource).^2; % Intensity of the field in the source plane (W/m^2)
%% ==========================================================================================================================
% |+|+|+|+| THE COMPUTATION OCCURS BETWEEN THIS LINE AND THE ONE LIKE IT BELOW |+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|
% ==========================================================================================================================
% Change of variablesa and variable initialization
% ------------------------------------------------
% h, below is a scale factor to change from the physical units (meters) to new units in which all physical lengths are scaled
% by h=1/sqrt(L*lambda). In the new units, the Fresnel integral becomes a pure Fourier transform multiplied by a phase factor.
% We now scale all physical lengths to the new units before performing the fourier tranform.
h = sqrt(1/BB/lambda); % scaling factor
dX=h*dx; dY=h*dy; % source plane spatial sampling (pixel) intervals in the new units
X = h*x; % source plane points (in scaled distance) at which to evaluate the source field
Y = h*y; % source plane points (in scaled distance) at which to evaluate the source field
dF = 1/dX/Nx; dG=1/dY/Ny; % corresponding spatial sampling interval in field plane after fft2
F=repmat(([0:Nx-1]-floor(Nx/2)) *dF,Ny,1); % Field plane, x-domain (in scaled length)
G=repmat(([0:Ny-1]-floor(Ny/2)).'*dG,1,Nx); % Field plane, y-domain (in scaled length)
df=dF/h; dg=dG/h; % field plane sampling intervals (in meters)
f = F/h; g = G/h; % Field plane, x and y-domains (in meters)
% Perform 2D FFT on and scale correctly
% -------------------------------------
Ufield = -1i/lambda/BB*exp(1i*pi*DD*((F).^2+(G).^2))... % * * * * * * HERE IT IS ! * * * * * * * * *
.*fftshift( fft2( exp(1i*pi*AA*(X.^2+Y.^2)).*Usource )*dx*dy );
Ifield = epsilon0*c/2*abs(Ufield).^2; % get the intensity
% =========================================================================================================================
% |+|+|+|+| EVERYTHING BELOW THIS LINE IS JUST FLUFF |+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|+|
% =========================================================================================================================
%% Report statistics of the calculation and Check energy conservation
% ------------------------------------------------------------------
% The total power in the source plane is the same as in the field plane. As a
% (non-exhaustive) check that things are scaled correctly, check that the total diffracted
% power is approximately equal to the power in the source plane.
inpow = trapz(trapz(Isource))*dx*dx;
outpow = trapz(trapz(Ifield))*df*dg;
disp(' ');
disp('--------------------------------------------------------------------------');
disp(['Source plane pixel size: dx*dy = ',num2str(1e6*dx,'%0.2f'),' * ',num2str(1e6*dy,'%0.2f'),' microns^2.']);
disp(['Source plane size: x*y = ',num2str(1e6*2*xmax,'%0.1f'),' * ',num2str(1e6*2*ymax,'%0.1f'),' microns^2']);
disp(['Field plane pixel size: df*dg = ',num2str(1e3*df,'%0.2f'),' * ',num2str(1e3*dg,'%0.2f'),' mm^2.']);
disp(['Field plane size: f*g = ',num2str(1e3*2*(f(end)-f(1)),'%0.1f'),' * ',num2str(1e3*2*(g(end)-g(1)),'%0.1f'),' mm^2']);
disp(' ');
disp(['Power in the source plane: Pin = ',num2str(inpow*1000),' mW']);
disp(['Power in the field plane: Pout = ',num2str(outpow*1000),' mW']);
disp('---------------------------------------------------------------------------');
disp(' ');
% Make a red colormap to use to display the laser beam intesity
% -------------------------------------------------------------
cb=colormap('bone'); % use the values in column one of the 'bone' colormap
cmred=[cb(:,2)*1, cb(:,2)*0.2, cb(:,2)*0.1]; % as the columns of the new colormap but scale them separately
% Display source plane field amplitude (Fig. 1)
% ---------------------------------------------
figure(1); % open a figure window
ax1 = pcolor(x*1e6,y*1e6,Isource/1000); % plot the intensity in mW/mm^2 the source plane as a fn of x,y (in microns)
xlabel('x ({\mu}m)'); % label the axes
ylabel('y ({\mu}m)');
axis square % make the figure display both axes to the same scale (no stretching)
set(ax1,'linestyle','none');
caxis([min([max(Isource(aperture))*0.9,...
min(Isource(aperture))]),max(Isource(aperture))]/1200);
% colormap('jet');
colormap(cmred);
shading interp;
aperturemat=zeros(Nx,Ny); % construct a matrix representing the aperture so that we can superimpose
aperturemat(aperture)=1e-15; % a contour plot of the aperture onto the intensity plot of the source field
hold on;
[tmp cont1]=contour(x*1e6,...
y*1e6,aperturemat,[0.5e-15,0.5e-15]); % Shows contours of the edge of the aperture
set(cont1,'color',[1 1 1]*0.5);
cbar1=colorbar;
ylabel(cbar1,'Intensity (mW/mm^2)');
hold off;
title(['Source Plane Intensity']);
% Display field plane field amplitude (Fig. 2)
% ---------------------------------------------
figure(2);
ax2 = pcolor(f*1e6,g*1e6,(Ifield/1000)); % plot the intensity in the field plane as a fn of physical field plane locations (in mm)
view(2); % If "surf" is substituted for "pcolor", this plot can be rotated from the top-view orientation
shading interp;
xlabel('x ({\mu}m)'); % label the axes
ylabel('y ({\mu}m)');
axis square
axis tight % make the figure display both axes to the same scale (no stretching)
set(ax2,'linestyle','none');
caxis(([min(min(Ifield))/12000 (max(max(Ifield)))/1250]));
title(['Diffracted Intensity in the Field Plane']);
colormap('jet');
%colormap(cmred);
cbar2=colorbar;
ylabel(cbar2,'Intensity (mW/mm^2)');
% hold on;
% contour(f*1e3,g*1e3,abs(Ifield)./max(max(abs(Ifield))),1/exp(2));
% hold off;
|
