% Example of partial temporal coherence
% partial temporal coherence revealed via Young's double-pinhole experiment
% with relatively delayed light beams with finite spectral bandwidth 
close all; clear all; clc

% Aperture function definition
circ=@(ri) (abs(ri)<=1.) ;

%% Some global parameters
lambda0=600e-9;     %central wavelength (m) 
c=3e8;                     %speed of light 
k0=2*pi/lambda0;    %center wavenumber 
nu0=c/lambda0;      %central frequency 

%% Spectral parameters of Gaussian lineshape 
N=51;                 % number of spectral components (odd) 
delnu=5e9;          % FWHM (Hz) spectral width of the light source 
b=delnu/(2*sqrt(log(2))); % FWHM scaling 
dnu=4*delnu/N;      % spectral resolution
%
%% source plane parameters 
L1=50e-3;           % side dimension of source plane (50mm)
M=128;               % # samples (even) 
dx1=L1/M;           % sample interval 
x1=-L1/2:dx1:L1/2-dx1; % source coords 
x1=fftshift(x1);     % shift x coord 
[X1,Y1]=meshgrid(x1,x1); 
%
%% beam parameters 
w=0.5e-3;             % beam waist radius = 0.5 mm 
dels=5e-3;          % transverse separation 
f=0.5;               % focal dist for Fraunhofer = 10 cm
lf=lambda0*f;   
ilf = 1/lf;

% Time delay mesh
Delay = 0:1e-2:8e-2; 
deld=2e-2;        % delay distance = 2 cm
ND = length(Delay);

%% Loop through spectral components 
stack = zeros(M,M,1,ND); 
for ii = 1:ND
  I2=zeros(M); 
  deld=Delay(ii) ; % delay distance
  for n=1:N 
    % spectral density function 
      nu=(n-(N+1)/2)*dnu+nu0; 
      S=1/(sqrt(pi)*b)*exp(-(nu-nu0)^2/b^2); 
      k=2*pi*nu/c; 
    % source 
      u=circ(sqrt((X1-dels/2).^2+Y1.^2)/w)... 
       +circ(sqrt((X1+dels/2).^2+Y1.^2)/w)... 
       *exp(j*k*deld); 
    %
    % Fraunhofer pattern 
      u2=ilf*(fft2(u))*dx1^2; 
    % weighted irradiance and sum 
      I2=I2+S*(abs(u2).^2)*dnu; 
  end 
%  
  I2=ifftshift(I2); % normalize/center irradiance 
  stack(:,:,1,ii) = I2;
end

%% Plot out the result on obs coords
x2=(-1/(2*dx1):1/L1:1/(2*dx1)-1/L1)*lf; 
y2=x2; 
%  
figure(1)          % irradiance image 
colormap('jet'); 
for ii = 1:ND
   subplot(3,3,ii)
   imagesc(x2,y2,stack(:,:,1,ii)); 
   xlabel('x (m)'); ylabel('y (m)');  
   axis square; axis xy; 
   title(['Delay: ' num2str(Delay(ii))]);
end
%  
figure(2)              % irradiance profile 
for ii = 1:ND
   subplot(3,3,ii)
   plot(x2,stack(M/2+1,:,1,ii)); 
   xlabel('x (m)'); ylabel('Irradiance');  
   title(['Delay: ' num2str(Delay(ii))]);
end