%
% fo.m Fourier Optics shell.

% This is a skeleton script to do Fourier Optics problems.  In it, we
% generate an object electric field, take the Fourier transform to
% simulate the field in the pupil plane, truncate it to simulate
% the aperture stop, and then inverse Fourier transform it to get
% the image back.  You can generate your own field, and do whatever
% you want in the pupil plane.  For example, you could introduce a
% quartic phase error to simulate spherical aberration, or use a
% strange shape for the aperture. 

rev=1.02; % Allows presetting some variables; NA, lambda, field
disp(['fo, Fourier Optics Demo Shell ',num2str(rev)]);
nxsamp=1024;nysamp=1024; % Number of samples 
xpix=0.2e-6;ypix=0.2e-6; % Pixel sizes in meters
if(~exist('lambda','var'));lambda=633e-9;end;  % light wavelength, meters
%% Set up x and y axes with origin in the center
xaxis=fftaxisshift(([0:nxsamp-1])*xpix); 
yaxis=fftaxisshift(([0:nysamp-1])*ypix);
%% Set up spatial frequency axes also centered. Note:  fftaxis
% generates a "frequency' axis from a "time" axis or in this case a
% spatial frequency axis from a transverse distance axis.
% fftaxisshift shifts the axis to match what fftshift does to the
% Fourier transform.
fxaxis=fftaxisshift(fftaxis(xaxis));
fyaxis=fftaxisshift(fftaxis(yaxis));
%% Set up direction cosine axes.  u and v will be limited by the
% aperture to u^2+v^2<NA^2
uaxis=fxaxis*lambda;
vaxis=fyaxis*lambda;

%% Make 2D arrays of x,y,fx,fy,u,v
[x,y]=meshgrid(xaxis,yaxis);
[fx,fy]=meshgrid(fxaxis,fyaxis);
[u,v]=meshgrid(uaxis,vaxis);
% Compute the square of the radius; Useful for field patterns that
% have circular symmetry.
rsquared=x.^2+y.^2;

%% Generate a field that is 1 inside some maximum radius rmax
rmax=5e-6;
if(~exist('field','var'));
  field=zeros(size(x));
field(find(rsquared<rmax.^2))=1;
end;

%% Take the Fourier transform.
Field=fftshift(fft2(ifftshift(field)));

% Image the original field (if it's not real, this needs to be
% modified to image real and imaginary parts or amplitude and phase
% separately. 
figure;imagesc(xaxis*1e6,yaxis*1e6,field);axis image;colorbar;
xlabel('x, \mum');ylabel('y, \mum');
title('field');

% Image against frequency units.  Show in dB
% Note the FT units are not useful in the sense of being electric
% fields but as long as we end up taking the IFT back to the object
% or image plane it doesn't matter.
figure;imagesc(fxaxis/1e3,fyaxis/1e3,20*log10(abs(Field)));
moose=caxis;caxis(moose(2)+[-50,0]);
axis image;colorbar;
xlabel('f_x, Cycles/mm');ylabel('f_y, Cycles/mm');
title('Magnitude of FT, dB');

% Image amplitude in dB against direction cosines
figure;imagesc(uaxis,vaxis,20*log10(abs(Field)));
moose=caxis;caxis(moose(2)+[-50,0]); % Limit the scale to 50 dB
axis image;colorbar;
xlabel('u');ylabel('v');
title('Magnitude of FT in Pupil, dB');

% Same for the phase
figure;imagesc(uaxis,vaxis,angle(Field)*180/pi);axis image;
caxis(180*[-1,1]);
colormap('hsv');colorbar;
xlabel('u');ylabel('v');
title('Phase in degrees');

%% Block the light that doesn't go through the NA

if(~exist('NA','var'));NA=0.95;end;
Aperture=double((u.^2+v.^2)<NA^2);
Field2=Field.*Aperture;
figure;imagesc(uaxis,vaxis,20*log10(abs(Field2)));
moose=caxis;caxis(moose(2)+[-50,0]); % Limit the scale to 50 dB
axis image;colorbar;
xlabel('u');ylabel('v');
title('Magnitude of FT after Aperture, dB');

%% Inverse Fourier Transform to get back to the field.
field2=fftshift(ifft2(ifftshift(Field2)));

figure;imagesc(xaxis*1e6,yaxis*1e6,abs(field2).^2);axis image;colorbar;
xlabel('x, \mum');ylabel('y, \mum');
title('Image');
