% fourier optics.m Homework for EECE4646
%  C. DiMarzio, Northeastern Univ. Spring 2017
%  Fourier Optics in Microscopy

% Parameters 

% laser
lambda=630e-9;k=2*pi/lambda;
% sphere with diameter diam and index contrast n-nbkg
%diam=80e-6;n=1.38;nbkg=1.36;
diam=80e-6;n=1.365;nbkg=1.36;
%xc=50e-6;yc=20e-6; % Center of sphere
xc=0e-6;yc=20e-6; % Center of sphere
% Pupil
NA=0.75;outer=1;  % NA is NA of objective.  outer is transmission
                  % of pupil within the NA

% set up a field of view.
fftsize=1024; % actual is twice this
pixel=lambda/4; % Pixel size

xax=[-fftsize+1:fftsize]*pixel;
yax=xax;
[x,y]=meshgrid(xax,yax);
rsq=(x-xc).^2+(y-yc).^2;

ppl=2*sqrt((diam/2)^2-rsq).*(rsq<(diam/2).^2);
figure;imagesc(xax,yax,ppl);
colorbar;axis image;
title('Physical Path Length');
xlabel('x');ylabel('y');

opl=ppl*(n-nbkg);
figure;imagesc(xax,yax,opl);colorbar;
colorbar;axis image;
title('Optical Path Length');
xlabel('x');ylabel('y');


E=exp(1j*k*opl);

figure;imagesc(xax,yax,imag(E));
moose=caxis;caxis([-1,1]*max(moose));
colorbar;axis image;
title('Imag(E)');
xlabel('x');ylabel('y');

% set up the pupil

% Field before pupil
Epupil=fft2(E);

fxax=fftaxisshift(fftaxis(xax));
fyax=fftaxisshift(fftaxis(yax));
uax=fxax*lambda;vax=fyax*lambda;
[u,v]=meshgrid(uax,vax);
% Pupil transmission function
Tpupil=double((u.^2+v.^2)<NA^2)*outer;

figure;imagesc(uax,vax,...
	       fftshift(20*log10(abs(Epupil))));
colorbar;axis image;
title('Before Pupil, dB');
xlabel('u');ylabel('v');
hold on;
circle(0,0,NA,'k-'); % Plot the pupil 
hold off;

% field after pupil
Epupil2=Epupil.*ifftshift(Tpupil);

% Image field
Eimage=ifft2(Epupil2);

figure;imagesc(xax,yax,abs(Eimage).^2);colorbar;axis image;
moose=caxis;caxis([0,1]*max(moose));
title('Image Irradiance');
xlabel('x');ylabel('y');

figure;imagesc(xax,yax,angle(Eimage)*180/pi);
colormap(hsv);colorbar;axis image;
moose=caxis;caxis([-1,1]*180);
title('Image Phase, Deg.');
xlabel('x');ylabel('y');

