% Fig_6_3to6
% Essential Semiconductor Laser Device Physics
% Ntot Fock state photon probability at output of
% lossless 50:50 beam splitter and indistinguishable photons
% line 31 is Ntot 
clear all
clf;
t=cputime;
%plotting parameters + fontsizes
FS      = 18;	%label fontsize 14
FSN     = 16;	%number fontsize 12
LW      = 1.5;	%linewidth

% Change default axes fonts.
set(0,'DefaultAxesFontName', 'Times New Roman')
set(0,'DefaultAxesFontSize', FSN)

% Change default text fonts.
set(0,'DefaultTextFontname', 'Times New Roman')
set(0,'DefaultTextFontSize', FSN)
%*********************************
% Constants and input parameters *
%*********************************
hbar    =1.0545715e-34;		% Reduced Planck's constant (J s)
eye     =complex(0.0,1.0);	% Square Root of -1
echarge =1.6021764e-19;		% Electron charge (C)
clight  =2.99792458e8;      % Seed of light (m s-1)
lambda0 =1500.0e-9;         % Resonant wavelength (m)
E0      =clight*hbar*2*pi/(echarge*lambda0); %Resonant energy (eV)

Ntot=8;            %total number of photons 2, 8, 32, 64
P=zeros(Ntot+2);    %preallocate and add extra row and column for 3D plot
                    %color square is pixel above vertex 

for ii=1:1:Ntot+1
    n1=double(ii)-1.0;  %number of photons input port 1
    n2=Ntot-n1;         %number of photons input port 2
    const=((-1.0)^n1)*((0.5)^(Ntot/2.0));
for    jj=1:Ntot+1
    N1=double(jj)-1.0;  %number of photons output port 1
    N2=Ntot-N1;         %number of photons output port 2
 
    Ptmp=0.0;
%********************************************************************** 
    for kk0=1:N1+1
        
fCC1=1.0;fCC2=1.0;fCC3=1.0;fCC4=1.0;%initialize flag

kk=double(kk0)-1.0;%k-counter

n1kk=n1-kk;
kk1=kk;
if (kk1 < 0.0); kk1 = 0.0; fCC1 = 0.0; %end%CC1 = 0;
elseif (n1kk < 0.0); n1kk = 0.0; fCC1 = 0.0; %end
elseif (kk > n1); fCC1 = 0.0; end%CC1 = 0; 
CC1= factorial(n1)/(factorial(kk1)*factorial(n1kk));
CC1=CC1*fCC1;
 
N1kk=N1-kk;
kk3=kk;
if (kk3) < 0.0; kk3 = 0.0; fCC3 = 0.0; %end%CC3 = 0; 
elseif (N1kk) < 0.0; N1kk = 0.0; fCC3 = 0.0; %end%CC3 = 0; 
elseif (kk > N1); fCC3=0.0; end%CC3 = 0; 
CC3= factorial(N1)/(factorial(kk3)*factorial(N1kk));
CC3=CC3*fCC3;

n2N1kk=n2-(N1-kk);
N1kk=N1-kk;
if (N1kk) < 0.0; N1kk = 0.0; fCC2 = 0.0; %end%CC2=0;
elseif (n2N1kk) < 0.0; n2N1kk = 0.0; fCC2 = 0.0; %end%CC2=0; 
elseif (N1kk > n2) ; fCC2 = 0.0; end%CC2=0; 
CC2= factorial(n2)/(factorial(N1kk)*factorial(n2N1kk)); 
CC2=CC2*fCC2;

N2n1kk=N2-(n1-kk);
n1kk=n1-kk;
if (n1kk) < 0.0; n1kk = 0.0; fCC4 = 0.0; %end%CC4 = 0; 
elseif (N2n1kk) < 0.0; N2n1kk = 0.0; fCC4 = 0.0; %end%CC4 = 0; 
elseif (n1kk > N2) ; fCC4=0.0; end%CC2=0; 
CC4= factorial(N2)/(factorial(n1kk)*factorial(N2n1kk));
CC4=CC4*fCC4;
 
 Ptmp=Ptmp+(const*((-1)^kk)*sqrt(CC1*CC2*CC3*CC4));
 
    end
    P(ii,jj)=(abs(Ptmp))^2;
end
end

%  P %print probability matrix 
%*********************************************************************
% Plot data
%*********************************************************************
xtickstep=1;
if Ntot>32; xtickstep=8;
elseif Ntot>16;xtickstep=4
elseif Ntot>8;xtickstep=2
end
    
ttl1=['Lossless 50:50 beam splitter: \it{n}\rm{_{tot}}= ',...
    num2str(Ntot)];
xp=[0:1:Ntot+1];%add extra pixel
xmin=xp(1);xmax=xp(end);
Pmax=1.1*max(P(1,:));
%********************************************************************
figure(1)
surf(xp,xp,P)
shading interp;  %do color interpolation
%                 az = 0;          %set up the angle of view
%                 el = 90;         %looking from directly above
%                 view(az, el);    %set the view
                colormap(jet); %winter cool autumn jet gray copper bone
xlabel('Input, (\it{n}\rm{_1}, \it{n}\rm{_2}= \it{n}\rm{_{tot}}- \it{n}\rm{_1})','Fontsize',FS,'Fontname','Times New Roman');
ylabel('Output, (\it{n}\rm{_3}, \it{n}\rm{_4}= \it{n}\rm{_{tot}}- \it{n}\rm{_3})','Fontsize',FS,'Fontname','Times New Roman');
zlabel('Probability','Fontsize',FS,'Fontname','Times New Roman');
title(ttl1,'Fontsize',FS,'Fontname','Times New Roman'); 
axis([ xmin xmax xmin xmax 0.0 Pmax]);
set(gca,'XTick',0:xtickstep:xmax-1);
set(gca,'YTick',0:xtickstep:xmax-1);

ttl2=['\it{n}\rm{_{tot}}= \it{n}\rm{_1}= ',...
    num2str(Ntot),' into port 1 (blue), \it{n}\rm{_1}= ',...
    num2str(ceil(Ntot/2)),' into port 1 (red)'];
%********************************************************************
figure(2)
plot(xp,P(1,:),'b','LineWidth',LW);
hold on
plot(xp,P(ceil((Ntot+1)/2),:),'r','LineWidth',LW);
xlabel('Output, (\it{n}\rm{_3}, \it{n}\rm{_4}= \it{n}\rm{_{tot}}- \it{n}\rm{_3})','Fontsize',FS,'Fontname','Times New Roman');
ylabel('Probability, \it{P}\rm{(}\it{n}\rm{_3}, \it{n}\rm{_4}= \it{n}\rm{_{tot}}- \it{n}\rm{_3})','Fontsize',FS,'Fontname','Times New Roman');
title(ttl2,'Fontsize',FS,'Fontname','Times New Roman'); 
axis([ xmin xmax-1 0 Pmax]);
set(gca,'XTick',0:xtickstep:xmax-1);
hold off
%********************************************************************
if Ntot>=2
    ttl3=['Lossless 50:50 beam splitter: N_{tot}=',...
    num2str(Ntot),', n_1=',...
    num2str(ceil(Ntot/2)),' into port 1'];
figure(3)
plot(xp,P(ceil((Ntot+1)/2),:),'LineWidth',LW)
xlabel('Output, (N_1, N_2=N_{tot}-N_1)');
ylabel('Probability, P(N_1, N_2=N_{tot}-N_1)');
title(ttl3); 
axis([ xmin xmax-1 0 Pmax]);
end
if Ntot>=4
ttl4=['Lossless 50:50 beam splitter: n_{tot}=',...
    num2str(Ntot),', n_1=',...
    num2str(ceil(Ntot/4)),' into port 1'];
%********************************************************************
figure(4)
plot(xp,P(ceil((Ntot+1)/4),:),'LineWidth',LW)
xlabel('Output, (n_3, n_4= n_{tot}-n_3)','Fontsize',FS,'Fontname','Times New Roman');
ylabel('Probability, P(n_3, n_4= n_{tot}-n_3)','Fontsize',FS,'Fontname','Times New Roman');
title(ttl4,'Fontsize',FS,'Fontname','Times New Roman'); 
axis([ xmin xmax-1 0 Pmax]);
set(gca,'XTick',0:xtickstep:xmax-1);
end

%********************************************************************
figure(5)%2D flat color plot
pcolor(xp,xp,P);
colormap(winter); %winter cool autumn jet gray copper bone
axis square;
%colorbar;
xlabel('Input, (\it{n}\rm{_1}, \it{n}\rm{_2}= \it{n}\rm{_{tot}}- \it{n}\rm{_1})','Fontsize',FS,'Fontname','Times New Roman');
ylabel('Output, (\it{n}\rm{_3}, \it{n}\rm{_4}= \it{n}\rm{_{tot}}- \it{n}\rm{_3})','Fontsize',FS,'Fontname','Times New Roman');
zlabel('Probability','Fontsize',FS,'Fontname','Times New Roman');
title(ttl1,'Fontsize',FS,'Fontname','Times New Roman'); 
axis([ xmin xmax xmin xmax ]);
set(gca,'XTick',0.5:xtickstep:xmax);
set(gca,'YTick',0.5:xtickstep:xmax);
set(gca,'XTickLabel',0:xtickstep:xmax-1);
set(gca,'YTickLabel',0:xtickstep:xmax-1);

%print probabilities to screen
P(1:Ntot+1,1:Ntot+1)
e=cputime-t