% Fig_6_8
% essential Semiconductor Laser Device Physics
% Photon number resolved Fock state photon at FP resonator
% Uses lossless 50:50 beam splitter to build model of FP resonator
% and gives probability of detecting n_phot indistinguishable photons
% using photon number resolving detector
% See Dowling PRA 80, 043822 (2009)

clear all
clf;
tcpu=cputime;
%plotting parameters + fontsizes
FS      = 14;	%label fontsize
FSN     = 12;	%number fontsize
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 *
%*********************************

% Universal constants
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)

ss=char('lightBlue','black','red','green','blue','magenta');
%plot curves in different colors
s=char('c','k','r','g','b','m');
%plot characters in different colors
% red, green, yellow, magenta, blue, black, white
% gray, darkGreen, orange, and lightBlue

npoints=900;    %number of points in plot
% resonator parameters
nLC=5.0;        %resonator length in units of lambda0
LC=nLC*lambda0; %resonator length
%t=1/sqrt(2.0); % On resonance, quarter wave relector dielectric has 
                % t^2=1/(1+((1-nr^2)/2nr)^2) and t^2=0.5 when nr=1+sqrt(2)

r2=0.7;         %modulus of reflection amplitude squared
r=sqrt(r2);     %modulus of reflection amplitude

nphot=4;                   %number of photons at input port 1
kphot=zeros(1,nphot);      %pre-allocate n-photon array

kmin=0.01*pi/lambda0;% 0.001
kmax=kmin*35.0; %50
pmax=1.1;                  %maximum on probability plot
%kstep=(kmax-kmin)/npoints;
T=zeros(1,npoints);        % initialize Transmission amplitude array
kk=linspace(0,kmax+kmin,npoints); % set up wave number k array
P=zeros(nphot,npoints);    %pre-allocate proability of detecting n-photons
PE=zeros(nphot,npoints);   %pre-allocate energy-weighted photon probability

econst=exp(-2.0*eye*sqrt(1-r2));
for i=1:npoints    
phi=kk(i)*LC;
%econst=1.0;
T(i)=((1-r2)*econst)/(((r2*econst*exp(2*eye*phi))-1.0));
end
T2=T.*conj(T);

for ii=1:nphot  %loop for number of transmitted n photons <= nphot
    const=factorial(nphot)/(factorial(ii)*factorial(nphot-ii));
   for i=1:npoints    
P(ii,i)=const*(T2(i)^ii)*((1.0-T2(i))^(nphot-ii));%probability
   end
PE(ii,:)=P(ii,:)*double(ii)/nphot;  %weight by number of photons
                                    %for energy spectrum  
end
Ptot=sum(PE);   %energy spectrum is sum of probabilities for all discrete 
                %photon number outputs weighted by number of photons

%*********************************************************
%Plot figures
%*********************************************************
ttl1=['Classical (blue), Quantum (red), \it{r}\rm = ',...
    num2str(r,'%5.3f'),', \it{L}\rm{_C} = ',num2str(nLC),...
    '\lambda_0, \lambda_0 = ',num2str(lambda0*1e6),' \mum'];
%*********************************************************
figure(1)
plot(kk*LC,T2,'b');              %Classical spectrum (blue)
hold on
plot(kk*LC,Ptot,'r:');           %Quantum spectrum (red)
xlabel('Phase, \phi(\it{kL}\rm{_C})',...
    'Fontsize',FS,'Fontname','Times New Roman');
ylabel('Transmission energy spectrum, |\it{T}\rm{(\phi)|^2}',...
    'Fontsize',FS,'Fontname','Times New Roman');
axis([ 0 kmax*LC 0 pmax ]);
title(ttl1,'Fontsize',FS,'Fontname','Times New Roman'); 
hold off
%*********************************************************
ttl2=['\it{P}\rm(\it{n}\rm{, }',...
    '\it{N_{\rm{tot}}}\rm{), }\it{N_{\rm{tot}}}\rm = ',...
    num2str(nphot),', \it{r}\rm = ',...
    num2str(r),', \it{L}\rm{_C} = ',num2str(nLC),...
    '\lambda_0 = ',num2str(LC*1e6),...
    ' \mum, \lambda_0 = ',num2str(lambda0*1e6),' \mum'];
%*********************************************************
figure(2)
for i=1:nphot
    cindex=mod(i,6)+1 ;          %color index for curves
    plot(kk*LC,P(i,:),s(cindex));%loop through number of kphotons <= nphot
                                 %and each curve in different color
ttlp=['\color{',ss(cindex,:),'}{n}\color{black} = ',num2str(i)];
text(4.1,(1.08-i*0.05),ttlp);    %color plot lables to match curves
    hold on
end
hold off
grid on
xlabel('Phase, \phi(\it{kL}\rm{_C})',...
    'Fontsize',FS,'Fontname','Times New Roman');
ylabel('Photon # resolving trans. prob., \it{P}\rm(\it{n}\rm{, }\it{N_{\rm{tot}}}\rm{)}',...
    'Fontsize',FS,'Fontname','Times New Roman');
axis([ 0 kmax*LC 0 pmax ]);
title(ttl2,'Fontsize',FS,'Fontname','Times New Roman'); 
%*********************************************************
% plot Figure 6.8
%*********************************************************
figure(3)
for i=1:nphot
    cindex=mod(i,6)+1 ;          %color index for curves
    plot(kk*LC,P(i,:),s(cindex),'LineWidth',LW);%loop through number of kphotons <= nphot
                                 %and each curve in different color
ttlp=['\color{',ss(cindex,:),'}{n}\color{black} = ',num2str(i)];
text(4.1,(1.08-i*0.05),ttlp);    %color plot lables to match curves
    hold on
end
plot(kk*LC,Ptot,'m--','LineWidth',LW)
hold off
grid on
xlabel('Phase, \phi(\it{kL}\rm{_C})',...
    'Fontsize',FS,'Fontname','Times New Roman');
ylabel('Photon # resolving trans. prob., \it{P}\rm(\it{n}\rm{, }\it{N_{\rm{tot}}}\rm{)}',...
    'Fontsize',FS,'Fontname','Times New Roman');
axis([ 0 kmax*LC 0 pmax ]);
title(ttl2,'Fontsize',FS,'Fontname','Times New Roman'); 
%*********************************************************
%log10 plot
%*********************************************************
pmax=max(max(-log10(P)));
figure(4)
for i=1:nphot
    cindex=mod(i,6)+1 ;          %color index for curves
    plot(kk*LC,-log10(P(i,:)),s(cindex),'LineWidth',LW);%loop through number of kphotons <= nphot
                                 %and each curve in different color
ttlp=['\color{',ss(cindex,:),'}{n}\color{black} = ',num2str(i)];
text(4.1,(pmax-i*0.7),ttlp);    %color plot lables to match curves
    hold on
end
plot(kk*LC,-log10(Ptot),'m--','LineWidth',LW)
hold off
grid on
xlabel('Phase, \phi(\it{kL}\rm{_C})',...
    'Fontsize',FS,'Fontname','Times New Roman');
ylabel('Photon # resolving trans. prob., -log_{10}\it{P}\rm(\it{n}\rm{, }\it{N_{\rm{tot}}}\rm{)}',...
    'Fontsize',FS,'Fontname','Times New Roman');
axis([ 0 kmax*LC 0 pmax ]);
title(ttl2,'Fontsize',FS,'Fontname','Times New Roman'); 

e=cputime-tcpu