% This is source plane intensities of different beams - 25.10.2011
clear;clc;clf;warning off MATLAB:divideByZero;close all
lamda = 1.55e-6; k = 2*pi/lamda;Ac = 1;zarr = [1 100 1e3 5e3];
rxst = -1.e-1:0.2e-2:1e-1;ryst = rxst;

rxset1 = rxst*0.2;ryset1 = rxset1;rxset2 = rxst*0.2;ryset2 = rxset2;rxset3 = rxst*0.8;ryset3 = rxset3;rxset4 = rxst*4;ryset4 = rxset4;
rxmat = [rxset1; rxset2; rxset3; rxset4];rymat = [ryset1; ryset2; ryset3; ryset4];
%%%%%%%%%%%%%%%%%%%%%%%%%%      Gaussian BEAM     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[sx,sy]=meshgrid(sxset,syset);[Fis, s] = cart2pol(sx,sy);
%%%%% Source settings %%%%%%%%%%
Fsarr = [1e15 1e15 1e15 1e15];%s = sqrt(sx.^2 + sy.^2);
%sxn = sx / 1e-2;syn = sy / 1e-2;
alfasarr = [0.1e-2 1e-2 2e-2 5e-2];alfaarr = 1./(k*alfasarr.^2) + j./(2*Fsarr);
%%%%%%%%%%%%% Source Beam %%%%%%%%%%%%%%%%%%
I_rarr = [];
for ib = 1:length(zarr)
z = zarr(ib);alfa = alfaarr(2);
rxset = rxmat(ib,:);ryset = rymat(ib,:);[rx,ry]=meshgrid(rxset,ryset);[Fir, r] = cart2pol(rx,ry);
ur = Ac/(1 + 2*j*alfa*z)*exp(-k*alfa*r.^2/(1 + 2*j*alfa*z)); % In cylindrical coordinates
I_r = ur.*conj(ur);maxgen = max(max(I_r));I_r = I_r / max(max(I_r));
I_rarr = cat(3,I_rarr,I_r); 
end
figure(1)
%%%%%%%% Four Individual 3D Intensity Plots %%%%%%%%%%%%%%%%%
for ib = 1:length(zarr)
I_r = I_rarr(:,:,ib);subplot(2,2,ib);
rxset = rxmat(ib,:);ryset = rymat(ib,:);[rx,ry]=meshgrid(rxset,ryset);rxn = rx / 1e-2;ryn = ry / 1e-2;
meshc(rxn,ryn,real(I_r));view([35 45]);colormap([0 0 0])
xlabel('\itr_x\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',-12);
ylabel('\itr_y\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',25)
zlabel('\fontname{Verdana}\fontsize{12}\itI_r_N','FontSize',12,'FontWeight','bold')
alfasstr = [' \it\alpha_s \rm\bf= ' num2str(alfasarr(ib)/1e-2,'%3.1f') ' cm '];
zstr = [' \itz \rm\bf= ' num2str(zarr(ib)/1e0,'%3.0f') ' m '];
leg1 = ['Receiver intensity plot for Gaussian beam with ' zstr];
title(leg1,'FontSize',10)
set(gcf,'Renderer','Zbuffer');set(gcf,'Color','White');set(gca,'FontSize',11);
axis ([min(min(rxn))*1 max(max(rxn))*1 min(min(ryn))*1 max(max(ryn))*1 0 1]);
end
figure(2)
%%%%%%%% Four Individual Contour Plots %%%%%%%%%%%%%%%%%
for ib = 1:length(zarr)
I_r = I_rarr(:,:,ib);subplot(2,2,ib);
rxset = rxmat(ib,:);ryset = rymat(ib,:);[rx,ry]=meshgrid(rxset,ryset);rxn = rx / 1e-2;ryn = ry / 1e-2;
[C,h] = contour(rxn,ryn,I_r,'-k');axis square;fax = 1;
axis ([min(min(rxn))*1 max(max(rxn))*1 min(min(ryn))*1 max(max(ryn))*1]);

line([min(min(rxn))*fax max(max(rxn))*fax],[0 0],'LineWidth',1.8,'color','black','LineStyle','--');
line([0 0],[min(min(ryn))*fax max(max(ryn))*fax],'LineWidth',1.8,'color','black','LineStyle','--');
xlabel('\itr_x\rm\bf axis in cm','FontSize',10,'FontWeight','bold');
ylabel('\itr_y\rm\bf axis in cm','FontSize',10,'FontWeight','bold');
alfasstr = [' \it\alpha_s \rm\bf= ' num2str(alfasarr(ib)/1e-2,'%3.1f') ' cm '];
zstr = [' \itz \rm\bf= ' num2str(zarr(ib)/1e0,'%3.0f') ' m '];
leg1 = ['Receiver intensity plot for Gaussian beam with ' zstr];
title(leg1,'FontSize',10)
%set(gcf,'Renderer','Zbuffer');set(gcf,'Color','White');set(gca,'FontSize',11);
set(h,'LineWidth',1.2);H = clabel(C,h);set(H,'FontSize',6,'LineWidth',3)
set(gca,'FontSize',11);set(gcf,'Color',[1 1 1])
end