diff options
Diffstat (limited to 'examples/fourier_optics_demonstration.m')
| -rw-r--r-- | examples/fourier_optics_demonstration.m | 198 |
1 files changed, 198 insertions, 0 deletions
diff --git a/examples/fourier_optics_demonstration.m b/examples/fourier_optics_demonstration.m new file mode 100644 index 0000000..36abeb2 --- /dev/null +++ b/examples/fourier_optics_demonstration.m @@ -0,0 +1,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;
\ No newline at end of file |