% This is receiver plane intensities of Gaussian beams - 17.04.2012
clear;clc;clf;warning off MATLAB:divideByZero;close all
lamda = 1.55e-6; k = 2*pi/lamda;Ac = 1;zarr = [0 500 1e3 5e3];
sxst = -1.e-1:0.2e-2:1e-1;syst = sxst;
rxst = -1.e-1:0.2e-2:1e-1;ryst = rxst;
sxset = sxst*0.2;syset = syst*0.2;
%%%%% Source settings %%%%%%%%%%
Fs = 1e15;alfas = 5e-2;alfa = 1./(k*alfas^2) + j./(2*Fs);
rxset1 = rxst*0.2;ryset1 = rxset1;rxset2 = rxst*0.4;ryset2 = rxset2;rxset3 = rxst*1;ryset3 = rxset3;
rxmat = [sxset;rxset1; rxset2; rxset3];rymat = [syset;ryset1; ryset2; ryset3];
%%%%%%%%%%%%%%%%%%%%%%%%%%      Source Gaussian BEAM     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[sx,sy]=meshgrid(sxset,syset);[Fis, s] = cart2pol(sx,sy);

%%%%%%%%%%%%% Source Beam %%%%%%%%%%%%%%%%%%
[sx,sy]=meshgrid(sxset,syset);[Fis, s] = cart2pol(sx,sy);
Us = Ac*exp(-k*alfa*s.^2); % In cylindrical coordinates
I_s = Us.*conj(Us);maxgen = max(max(I_s));I_s = I_s / max(max(I_s));
I_rarr = I_s; 

for ib = 2:length(zarr)
z = zarr(ib);
rxset = rxmat(ib,:);ryset = rymat(ib,:);[rx,ry]=meshgrid(rxset,ryset);[Fir, r] = cart2pol(rx,ry);
%%%% Expression in terms of alpha and electric field
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);
%%%% Intensity expression in terms of alpha_s and F_s
I_r = Ac^2*k^2*alfas^4*Fs^2/(k^2*alfas^4*Fs^2 - 2*k^2*alfas^4*Fs*z + 4*Fs^2*z^2 + k^2*alfas^4*z^2)* ....
    exp(-r.^2*2*k^2*alfas^2*Fs^2/(k^2*alfas^4*Fs^2 - 2*k^2*alfas^4*Fs*z + 4*Fs^2*z^2 + k^2*alfas^4*z^2));
I_r = I_r / maxgen;
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)
if ib == 1;xlabel('\its_x\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',-12);end
if ib == 1;ylabel('\its_y\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',25);end
zlabel('\fontname{Verdana}\fontsize{12}\itI_r_N','FontSize',12,'FontWeight','bold')

zstr = [' \itz \rm\bf= ' num2str(zarr(ib)/1e0,'%3.0f') ' m '];
leg1 = ['Receiver intensity plot for Gaussian beam with ' zstr];if ib == 1;leg1 = ['Source intensity plot for Gaussian beam with ' zstr];end
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 max(max(I_r))]);
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');
if ib == 1;xlabel('\its_x\rm\bf axis in cm','FontSize',10,'FontWeight','bold');end
if ib == 1;ylabel('\its_y\rm\bf axis in cm','FontSize',10,'FontWeight','bold');end
zstr = [' \itz \rm\bf= ' num2str(zarr(ib)/1e0,'%3.0f') ' m '];
leg1 = ['Receiver intensity plot for Gaussian beam with ' zstr];if ib == 1;leg1 = ['Source intensity plot for Gaussian beam with ' zstr];end
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
figure(3)
%%%%%%%% Four Individual Crosssection Plots %%%%%%%%%%%%%%%%%
for ib = 1:length(zarr)
I_r = I_rarr(:,:,ib);subplot(2,2,ib);
rxset = rxmat(ib,:);ryset = rymat(ib,:);
R = sqrt(rxset.^2+ryset.^2);RR = [-R(1:round(numel(R)/2)) R(round(numel(R)/2)+1:end)];RRn = RR / 1e-2;
plot(RRn,diag(I_r),'-k','LineWidth',2);
set(gcf,'Renderer','Zbuffer');set(gcf,'Color','White');set(gca,'FontSize',14);grid on;
axis ([min(RRn)*1 max(RRn)*1 0 max(max(I_r))*1.05]);
xlabel('\itr\rm\bf axis in cm','FontSize',10,'FontWeight','bold');
ylabel('\fontname{Verdana}\fontsize{12}\itI_r_N','FontSize',12,'FontWeight','bold')
if ib == 1;xlabel('\its\rm\bf axis in cm','FontSize',10,'FontWeight','bold');end
if ib == 1;ylabel('\fontname{Verdana}\fontsize{12}\itI_s_N','FontSize',12,'FontWeight','bold');end
end