% This is Modified Bessel and Bessel beam source plane intensities - 28.12.2006
% Note that summation is implemented
function Bessel_sou
clear;clc;clf;warning off MATLAB:divideByZero;close all
lamda = 1.55e-6; k = 2*pi/lamda;
narr = [0 1 2 3] ;
sxset = -1.e-1:0.3e-2:1e-1;syset = sxset;
sxset = sxset*0.4;syset = sxset; 
%%%%%%%%%%%%%%%%%%%%%%%%%%      Bessel BEAM     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% Source settings %%%%%%%%%%
Fs = 1e15;
alfas = 1e-2;alfa = 1/(k*alfas^2) + j/(2*Fs);
aBarr = [1/alfas 1/alfas 1/alfas 1/alfas];
sxsetarr = [1.0*sxset;1.5*sxset;2*sxset;4*sxset];sysetarr = sxsetarr; 
sxsetarr = [1.0*sxset;1.0*sxset;1*sxset;1*sxset];sysetarr = sxsetarr;
%%%%%%%%%%%%% Source Beam %%%%%%%%%%%%%%%%%%

I_sarr = [];
for in = 1:length(narr)
    n = narr(in);aB = aBarr(in);
    sxset = sxsetarr(in,:);syset = sysetarr(in,:);
[sx,sy]=meshgrid(sxset,syset);[Fis, s] = cart2pol(sx,sy);    
u_s = besseli(n,aB*s).*exp(-k*alfa*s.^2).*exp(-j*n*Fis); 
I_s = u_s.*conj(u_s);maxgen = max(max(I_s));I_s = I_s / max(max(I_s));
I_sarr = cat(3,I_sarr,I_s); 
end

%%%%%%%% Four Individual Intensity Plots %%%%%%%%%%%%%%%%%
for is = 1:length(narr)
sxset = sxsetarr(is,:);syset = sysetarr(is,:);[sx,sy]=meshgrid(sxset,syset); 
sxn = sx / 1e-2;syn = sy / 1e-2;
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',35)
zlabel('\fontname{Verdana}\fontsize{12}\itI_s_N','FontSize',12,'FontWeight','bold')
nstr = [' \itn \rm\bf= ' num2str(narr(is),'%3.0f') ' '];
aBstr = [' \ita_B = \rm\bf' num2str(aBarr(is),'%3.1f') ' m^-^1 '];

leg1 = ['Source intensity plot for a modified Bessel Gaussian beam with' nstr aBstr];

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
alfasstr = [' \it\alpha_s = \rm\bf' num2str(alfas/1e-2,'%3.1f') ' cm '];

leg1 = alfasstr;
text(min(min(sxn))*0.3,min(min(syn))*0.3,max(max(I_s))*(+1.3),leg1,'FontSize',10,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');