% 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;