%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mutual Intensity function formed by inputing a LED light beam into an optical system. 
% The LED light is incoherent and quasi-monochromatic (i.e., temporally coherent).
% First, the van Cittert-Zernike Theorem is invoked to estimate the mutual intensity function 
% after an initial  free space propagation.
% Paraxial approximation is used for the optical system.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all ;
clear; clc

%% Source plane parameters 
% Define a hard rectangular aperture
rect=@(s) (abs(s)<=0.5) ; 
L= 10 ;          % Lateral dimension (10 mm) on the source plane 
N=1024;       % #samples (even) 
dx=L/N;        % sampling interval= 10 um 

xr=-L/2:dx:L/2-dx ;    % source coordinates 
x1=xr;
x2=x1 ;

% With xr, prepare a mesh of spatial frequency
% ur = -N/(2*L):1/L:(N/2-1)/L ;
% [u1,v1]=meshgrid(ur,ur);

% Parameters of optical source
wl=0.65*10^(-3);        % wacelength (mm)
k=2.0*pi/wl;                 % wavenumber
aa=1.0 ;                           % Source diameter (mm)

% lens data
R=50.0 ;                       % Radius of curvature, mm
d0=9.0 ;                       % Max thickness
n=1.52 ;                       % Refractive Index
Dia=50. ;                      % Lens diameter                  

%% Phase shifts induced by a planar-spherical lens
% with a central thickness of d0.
d1=d0-R*(1-sqrt(1-(x1.^2)/R^2));
d2=d0-R*(1-sqrt(1-(x2.^2)/R^2));
dd1=rect(x1/Dia).*(n*d1+(d0-d1));
dd2=rect(x2/Dia).*(n*d2+(d0-d2));

t1=exp(j*k*dd1);
t2=exp(-j*k*dd2);

%% Prepare mutual intensity right before lens1
% by invoking van Citter-Zernike Theorem
z1=93.;    
bb=wl*z1/aa ;  
p=1.0+1j;
J11=zeros(N,N,'like',p);
J22=zeros(N,N,'like',p);

for idx=1:N
  sd12(:)=x2(idx)^2-x1(:).^2 ;
  x12(:)=x2(idx)-x1(:);
  J11(:,idx)=(aa/(pi*z1^2))*exp(-j*pi/(wl*z1)*sd12(:)).*sinc(x12(:)./bb);
end

% Note: The radius of coherent spots developed on the first lens
% r_c=lambda*z/r_s=0.65*z1/0.5=120 um
% dx=L/N1= 10 um--> The spot is sufficiently sampled (by 12 times).
%
% Plot of mutual intensity function right before lens1
[X,Y]=meshgrid(xr,xr) ;
JJ=abs(J11)./max(max(abs(J11))) ;        
font_size=14;
line_width=2;
Marker_size=10;
figure(1), mesh(X,Y,JJ,'Markersize',Marker_size,'linewidth',line_width)
view(0,90);
axis equal;
colorbar;
set(gca,'fontsize',font_size,'linewidth',line_width)
set(gcf,'Color','w')
xlabel('x1(mm)','fontsize',font_size)
ylabel('x2(mm)','fontsize',font_size)
title('Coherence Factor right before Lens 1')
xlim([-2 2]); ylim([-2 2])

%% Mutual intensity right behind lens1
for idx=1:N
  J11(:,idx)=t2(idx)*J11(:,idx).*t1(:) ;
end
% Free space propagation to lens 2: for Mutual intensity
z2=30. ;
% 1) Quadrature phase multiplication
for idx=1:N
    sd12 = x1.^2-x2(idx)^2 ;
    J11(:,idx)=exp(j*pi/(wl*z2)*sd12(:)).*J11(:,idx) ;
end
% 2) Scaled Fourier Transform
J22=fft2(J11);
J22=fftshift(J22);
J22=fliplr(J22);
% 3) Quadrature phase multiplication
for idx=1:N
    sd12 = x1.^2-x2(idx)^2 ;
    J11(:,idx)=exp(j*pi/(wl*z2)*sd12(:)).*J22(:,idx) ;
end
% Mutual intensity right after lens2
for idx=1:N
  J11(:,idx)=t2(idx)*J11(:,idx).*t1(:) ;
end

% Plot of mutual intensity function right behind Lens 2
%JJ=abs(J11)./max(max(abs(J11))) ;        
%font_size=14;
%line_width=2;
%Marker_size=10;
%figure(2),mesh(X,Y,JJ,'Markersize',Marker_size,'linewidth',line_width)
%view(0,90);
%axis equal;
%colorbar;
%set(gca,'fontsize',font_size,'linewidth',line_width)
%set(gcf,'Color','w')
%xlabel('x1(mm)','fontsize',font_size)
%ylabel('x2(mm)','fontsize',font_size)
%title('Coherence Factor right behind Lens 2')
%xlim([-5 5]); ylim([-5 5])

%% Mutual intensity propagating to a slit 
% (or a pinhole in a configuration of axial assymetry)
z3=z1;
% 1) Quadrature phase multiplication
for idx=1:N
    sd12 = x1.^2-x2(idx)^2 ;
    J11(:,idx)=exp(j*pi/(wl*z3)*sd12(:)).*J11(:,idx) ;
end
% 2) Scaled Fourier Transform
J22=fft2(J11);
J22=fftshift(J22);
J22=fliplr(J22);
% 3) Quadrature phase multiplication
for idx=1:N
    sd12=x1.^2-x2(idx)^2 ;
    J11(:,idx)=exp(j*pi/(wl*z3)*sd12(:)).*J22(:,idx);
end

%% After the slit (or pinhole in a configuration of axial assymetry)
% Pinhole function with a diameter of 0.2 mm
P=0.2; 
PF1=rect(x1/P);
PF2=rect(x2/P); 
%
for idx=1:N
  J11(:,idx)=PF2(idx)*J11(:,idx).*PF1(:);
end

% Plot the mutual intensity function right behind the slit
JJ=abs(J11)./max(max(abs(J11))) ;        
font_size=14;
line_width=2;
Marker_size=10;
figure(2),mesh(X,Y,JJ,'Markersize',Marker_size,'linewidth',line_width)
view(0,90);
axis equal;
colorbar;
set(gca,'fontsize',font_size,'linewidth',line_width)
set(gcf,'Color','w')
xlabel('x1(mm)','fontsize',font_size)
ylabel('x2(mm)','fontsize',font_size)
title('Coherence Factor right behind the slit')
xlim([-1 1]); ylim([-1 1])

%% Free propagating the Mutual intensity to Lens 3
z4=z1;
% 1) Quadrature phase multiplication
for idx=1:N
    sd12 = x1.^2-x2(idx)^2 ;
    J11(:,idx)=exp(j*pi/(wl*z4)*sd12(:)).*J11(:,idx) ;
end
% 2) Scaled Fourier Transform
J22=fft2(J11);
%J22=fftshift(J22);
J22=fliplr(J22);
% 3) Quadrature phase multiplication
for idx=1:N
    sd12 = x1.^2-x2(idx)^2 ;
    J11(:,idx)=exp(j*pi/(wl*z4)*sd12(:)).*J22(:,idx) ;
end
% Mutual intensity right after Lens 3
for idx=1:N
  J11(:,idx)=t2(idx)*J11(:,idx).*t1(:) ;
end

%% Figures Preparation 
% Plot coherence factor on x1-x2 meshgrid
JJ=abs(J11)./max(max(abs(J11))) ;        
font_size=14;
line_width=2;
Marker_size=10;
figure(3),mesh(X,Y,JJ,'Markersize',Marker_size,'linewidth',line_width)
view(0,90);
axis equal;
colorbar;
set(gca,'fontsize',font_size,'linewidth',line_width)
set(gcf,'Color','w')
xlabel('x1(mm)','fontsize',font_size)
ylabel('x2(mm)','fontsize',font_size)
title('Coherence Factor right after Collimated Lens')
xlim([-2 2]); ylim([-2 2])

% Plot coherent factor profile along the x1-axis 
% (Fix x2 on the optical axis)
N12=N/2;
Jcd=JJ(:,N12);
font_size=14;
line_width=2;
Marker_size=10;
figure(4),plot(x1,Jcd,'Markersize',Marker_size,'linewidth',line_width)
set(gca,'fontsize',font_size,'linewidth',line_width)
set(gcf,'Color','w')
xlabel('x1(mm)','fontsize',font_size)
ylabel('Jc','fontsize',font_size)
title('Coherence Factor Profile along x1-axis')
xlim([-2 2])