% This is Laguerre source plane intensities - 20.10.2006
function Laguerre_sou
clear;clc;clf;warning off MATLAB:divideByZero;close all
linw = [2 2 3.0 3.0 2 2 2 2 2] ; sitil = ['-k ';'--k';':k ';'-.k';'-*k';':dk';'-.k';':k '];
lamda = 1.55e-6; k = 2*pi/lamda;
narr = [1 2 0 2] ;marr = [0 0 1 1];%marr = [0 0 0 0];
sxset = -1.e-1:0.15e-2:1e-1;syset = sxset;
%sxset = sxset*0.8;syset = sxset; 
sxset = sxset*0.3;syset = sxset; 
%%%%%%%%%%%%%%%%%%%%%%%%%%      Bessel BEAM     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sx,sy] = meshgrid(sxset,syset);[Fis, s] = cart2pol(sx,sy);
%%%%% Source settings %%%%%%%%%%
alfas = 1e-2;alfa = 1/(k*alfas^2);alfacs = alfas/sqrt(2); 
alfasarr = [1e-2 1e-2 1e-2 1e-2 1e-2];alfasarr = alfasarr*1;maxgenarr = [];

alfasstr = [' \it\alpha_s = \rm\bf' num2str(alfas/1e-2,'%4.1f') ' cm '];
nmstr1 = [' \itn\rm\bf = ' num2str(narr(1),'%3.0f') ' , ' '\itm\rm\bf = ' num2str(marr(1),'%3.0f')];
nmstr2 = [' \itn\rm\bf = ' num2str(narr(2),'%3.0f') ' , ' '\itm\rm\bf = ' num2str(marr(2),'%3.0f')];
nmstr3 = [' \itn\rm\bf = ' num2str(narr(3),'%3.0f') ' , ' '\itm\rm\bf = ' num2str(marr(3),'%3.0f')];
nmstr4 = [' \itn\rm\bf = ' num2str(narr(4),'%3.0f') ' , ' '\itm\rm\bf = ' num2str(marr(4),'%3.0f')];
legc = alfasstr;%{alfasstr;' ';Lstr};
%%%%%%%%%%%%% Source Beam %%%%%%%%%%%%%%%%%%
I_sarr = [];
for is = 1:length(narr)
    n = narr(is);m = marr(is);
    alfas = alfasarr(is);alfa = 1/(k*alfas^2);
%%% In radial coordinates    
u = (-j)^m*(sqrt(2)*s/alfas).^m.*exp(-k*alfa*s.^2).*exp(j*m*Fis).*LaguerreH(m,n,2*s.^2/alfas^2);

I_s = u.*conj(u);maxgen = max(max(I_s));I_s = I_s / max(max(I_s));maxgenarr = [maxgenarr maxgen];
sxn = sx / 1e-2;syn = sy / 1e-2;I_sarr = cat(3,I_sarr,I_s);end
%%%%%%%%%%%%% Receiver Beam %%%%%%%%%%%%%%%%%%

%%%% Plot along s_x axis. 
figure(1)
I_s = I_sarr(fix((length(sxset) + 1)/2),:,1);plot(sxset/1e-2,I_s,'-k','LineWidth',2);hold on;
I_s1 = I_sarr(fix((length(sxset) + 1)/2),:,2);plot(sxset/1e-2,I_s1,'--k','LineWidth',3);hold on
I_s2 = I_sarr(fix((length(sxset) + 1)/2),:,3);plot(sxset/1e-2,I_s2,':k','LineWidth',3);hold on
I_s3 = I_sarr(fix((length(sxset) + 1)/2),:,4);plot(sxset/1e-2,I_s3,'-.k','LineWidth',3);hold off

set(gcf,'Color',[1 1 1]);set(gca,'FontSize',16)
xlabel('\its\rm\bf  axis in cm','FontSize',14,'FontWeight','bold');
zlabel('\fontname{Verdana}\fontsize{12}\it I_s_N ','FontSize',12,'FontWeight','bold');%grid on;
%save('C:\Documents and Settings\HTERemote\My Documents\MATLAB\Scin\Area_scin\BesselJ0_I0\plot_data\Fig22I','sxset','I_s','I_s1','I_s2','I_s3','-ASCII')
axis ([min(sxset/1e-2)*1.05 max(sxset/1e-2)*1.05 0 max(I_s)*1.05]);
text(max(sxset/1e-2)*0.9,max(I_s)*.7,legc,'FontSize',20,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');
text(max(sxset/1e-2)*0.1,max(I_s)*.9,nmstr1,'FontSize',20,'FontWeight','bold','BackgroundColor','white');
text(max(sxset/1e-2)*0.2,max(I_s)*.9,nmstr2,'FontSize',20,'FontWeight','bold','BackgroundColor','white');
text(max(sxset/1e-2)*0.3,max(I_s)*.9,nmstr3,'FontSize',20,'FontWeight','bold','BackgroundColor','white');
text(max(sxset/1e-2)*0.4,max(I_s)*.9,nmstr4,'FontSize',20,'FontWeight','bold','BackgroundColor','white');
 
%%%%%%%% Four Individual Intensity Plots %%%%%%%%%%%%%%%%%
figure(2)
for is = 1:length(narr)
I_s = I_sarr(:,:,is);subplot(2,2,is);meshc(sxn,syn,real(I_s));view([35 45]);colormap([0 0 0])
xlabel('\its_x\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',-14);
ylabel('\its_y\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',30)
zlabel('\fontname{Verdana}\fontsize{12}\itI_s_N','FontSize',12,'FontWeight','bold')
nstr = [' \itn \rm\bf= ' num2str(narr(is),'%3.0f') ' '];
mstr = [' \itm \rm\bf= ' num2str(marr(is),'%3.0f') ' '];
leg1 = ['Source intensity plot for Laguerre-Gaussian beam with ' nstr '\rm\bf' mstr];
%leg1 = ['Source intensity picture for a Bessel-Gaussian beam with '];
title(leg1,'FontSize',10)
set(gcf,'Renderer','Zbuffer');set(gcf,'Color','White');set(gca,'FontSize',11);
axis ([min(min(sxn))*1 max(max(sxn))*1 min(min(syn))*1 max(max(syn))*1 0 1]);
end
text(min(min(sxn))*0.9,min(min(syn))*0.1,max(max(I_s))*1.3,legc,'FontSize',10,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');
function y = funb(a,b)
y = factorial(a)/factorial(b)/factorial(a-b);