% This is receiver plane beam sizes of Gaussian beams - 23.04.2012
clear;clc;clf;warning off MATLAB:divideByZero;close all
lamda = 1.51e-6; k = 2*pi/lamda;zarr = 0:0.1e2:5e3;
% lamda = 0.5e-6; k = 2*pi/lamda;zarr = 0:0.1e2:6e3;
Fs = 100;alfas = 1.2e-2;

%%%% Beam size expression in terms of alpha_s and F_s
alfar = sqrt((k^2*alfas^4*Fs^2 - 2*k^2*alfas^4*Fs*zarr + 4*Fs^2*zarr.^2 + k^2*alfas^4*zarr.^2)/(k^2*alfas^2*Fs^2));
alfaf = 2*Fs/(k*alfas)
zB = k^2*alfas^4*Fs/(4*Fs^2 + k^2*alfas^4) %%% Distance of beam waist from source plane
alfaB = sqrt(4*alfas^2*Fs^2/(4*Fs^2 + k^2*alfas^4))  %%% Beam waist
zR1 = (k^2*alfas^4*Fs - 2*k*alfas^2*Fs^2)/(4*Fs^2 + k^2*alfas^4)  %%%% Rayleigh range for convergent beam
zR2 = (k^2*alfas^4*Fs + 2*k*alfas^2*Fs^2)/(4*Fs^2 + k^2*alfas^4)  %%%% Rayleigh range for convergent beam
zR = 0.5*k*alfas^2
figure(1)
plot(zarr/1e3,alfar/1e-2,'-k','LineWidth',2);
set(gcf,'Renderer','Zbuffer');set(gcf,'Color','White');set(gca,'FontSize',14);grid on;
axis ([min(zarr/1e3)*1 max(zarr/1e3)*1 0 max(alfar/1e-2)*1.05]);
ylabel('\it\alpha_r\rm\bf in cm - Beam size','FontSize',16,'FontWeight','bold');
xlabel('\itz\rm\bf in km - Propagation distance','FontSize',16,'FontWeight','bold')
figure(2)
% %%%%%%%% Illustration for beam diffraction %%%%%%%%%%%%%%%%%
alfar1 = -alfar;
zeroc = zeros(1,length(alfar));
plot(zarr/1e3,alfar/1e-2,'-k','LineWidth',3);hold on;plot(zarr/1e3,alfar1/1e-2,'-k','LineWidth',3);hold on
plot(zarr/1e3,zeroc,'--k','LineWidth',2);
set(gcf,'Renderer','Zbuffer');set(gcf,'Color','White');set(gca,'FontSize',16);
ylabel('\its \rm\bf or \itr \rm\bf axis in cm - Beam size borders','FontSize',16,'FontWeight','bold');
xlabel('\itz\rm\bf in km - Propagation distance','FontSize',16,'FontWeight','bold')
axis ([min(zarr/1e3)*1 max(zarr/1e3)*1 -max(alfar/1e-2)*1.05 max(alfar/1e-2)*1.05]);grid on