% This is Bessel source plane intensities - 20.10.2006
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] ;%narr = [4 4 4 4];
sxset = -1.e-1:0.2e-2:1e-1;syset = sxset;
%sxset = sxset*0.8;syset = sxset; 
sxset = sxset*1;syset = sxset; 
%%%%%%%%%%%%%%%%%%%%%%%%%%      Bessel BEAM     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sx,sy] = meshgrid(sxset,syset);[Fis, s] = cart2pol(sx,sy);
%%%%% Source settings %%%%%%%%%%
Fs = 1e15;
alfas = 1e-2;alfasarr = [1e-2 1e-2 1e-2 1e-2];alfasarr = alfasarr*1;
alfaarr = 1./(k*alfasarr.^2) + j/(2*Fs);
aBarr = [1/alfasarr(1) 1/alfasarr(1) 1/alfasarr(1) 1/alfasarr(1)];
%aBarr = [3.76/alfa0arr(1) 3.8/alfa0arr(1) 3.74/alfa0arr(1) 3.77/alfa0arr(1)];
%%%%%%%%%%%%% Source Beam %%%%%%%%%%%%%%%%%%
I_sarr = [];
for is = 1:length(narr)
alfa = alfaarr(is);
%    n = narr(is+0);aB = 0.01*aBarr(4); % For Fig.1 
    n = narr(is);aB = aBarr(is); % For Fig.2
u = besselj(n,aB*s).*exp(-j*n*Fis); % New definition
%u = besselj(n,aB*sqrt(sx.^2 + sy.^2)).*exp(-k*alfa*sx.^2 - k*alfa*sy.^2).*(sx - j*sy).^n./(sx.^2 + sy.^2).^(0.5*n); % In cartesian coordinates
I_s = u.*conj(u);maxgen = max(max(I_s));I_s = I_s / max(max(I_s));
%I_sxy = u_xy.*conj(u_xy);I_sxy = (I_sxy) / max(max(I_sxy));
sxn = sx / 1e-2;syn = sy / 1e-2;I_sarr = cat(3,I_sarr,I_s); 
end
figure(1)

%%%%%%%% Four Individual Intensity Plots %%%%%%%%%%%%%%%%%
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') ' '];
aBstr = [' \ita_B = \rm\bf' num2str(aBarr(is),'%3.1f') ' m^-^1 '];
%leg1 = ['Source intensity plot for Bessel-Gaussian beam with ' nstr];
leg1 = ['Source intensity plot for a 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 max(max(I_s))]);
end
alfasstr = [' \it\alpha_s = \rm\bf' num2str(alfas/1e-2,'%3.1f') ' cm '];
aBstr = [' \it\fontname{century gothic}a\fontname{helvetica} = \rm\bf' num2str(aBarr(3),'%3.1f') ' m^-^1 '];
aBstr = [' \ita_B = \rm\bf' num2str(aBarr(4),'%3.1f') ' m^-^1 '];
nstr = [' \itn \rm\bf= ' num2str(narr(2),'%3.0f') ' '];
leg1 = [aBstr ',' alfasstr];leg1 = alfasstr;
text(min(min(sxn))*0.9,min(min(syn))*0.1,max(max(I_s))*(+1.3),leg1,'FontSize',10,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');