% This is source plane intensities of sinusoidal and hyperbolic Gaussian beams - 25.10.2011
%function SinoHyp_ECE635
clear;clc;clf;warning off MATLAB:divideByZero;close all
lamda = 1.55e-6; k = 2*pi/lamda;
sxst = -1.e-1:0.2e-2:1e-1;syst = sxst;rxst = sxst;ryst = sxst;

sxset1 = sxst*0.2;syset1 = sxset1;sxset2 = sxst*0.4;syset2 = sxset2;sxset3 = sxst*0.2;syset3 = sxset3;sxset4 = sxst*0.3;syset4 = sxset4;
sxmat = [sxset1; sxset2; sxset3; sxset4];symat = [syset1; syset2; syset3; syset4];
rxset1 = rxst*0.4;ryset1 = rxset1;rxset2 = rxst*0.6;ryset2 = rxset2;rxset3 = rxst*1.2;ryset3 = rxset3;rxset4 = rxst*4;ryset4 = rxset4;
rxmat = [rxset1; rxset2; rxset3; rxset4];rymat = [ryset1; ryset2; ryset3; ryset4];
%%%%%%%%%%%%%%%%%%%%%%%%%%      Gaussian BEAM     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ALarr = [1 1;j -j;1 1;1 -1];Dmat = [0 25 50 200 200 250];
alfasarr = [1e-2 1e-2;1e-2 1e-2;1e-2 1e-2;1e-2 1e-2];
Dsarr = [j*Dmat(4) -j*Dmat(4);j*Dmat(4) -j*Dmat(4);j*Dmat(4) -j*Dmat(4);
         j*Dmat(4) -j*Dmat(4)];
Fsarr = [5e15 5e15;5e15 5e15;5e15 5e15;5e15 5e15];
alfaarr = 1./(k*alfasarr.^2) + j./(2*Fsarr);zarr = [0.2e3 0.5e3 1e3 3e3];
%%%% Legends %%%%%%
legtitle = [' cos Gaussian beam         ';
            ' sine Gaussian beam        ';
            ' cosh Gaussian beam        ';
            ' sinh Gaussian             '];
%%%%%%%%%%%%% Source Beam %%%%%%%%%%%%%%%%%%
I_sarr = [];
ib = 1; % Select one beam
sxset = sxmat(ib,:);syset = symat(ib,:);[sx,sy]=meshgrid(sxset,syset);[Fis, s] = cart2pol(sx,sy);us = zeros(size(s));

for iL = 1:length(ALarr(1,:))
AL = ALarr(ib,iL);alfa = alfaarr(ib,iL);Ds = Dsarr(ib,iL);
us = us + AL*exp(-k*alfa*s.^2 + (sin(Fis) + cos(Fis))*Ds.*s); % In cylindrical coordinates
end
I_s = us.*conj(us);maxgen = max(max(I_s));I_s = I_s / max(max(I_s));
I_sarr = cat(3,I_sarr,I_s); 

%%%%%%%%%%%%% Receiver Beam %%%%%%%%%%%%%%%%%%

I_rarr = [];
for iz = 1:length(zarr)     
z = zarr(iz);
rxset = rxmat(iz,:);ryset = rymat(iz,:);[rx,ry]= meshgrid(rxset,ryset);[Fir, r] = cart2pol(rx,ry);ur = zeros(size(r));
for iL = 1:length(ALarr(1,:))
AL = ALarr(ib,iL);alfa = alfaarr(ib,iL);Ds = Dsarr(ib,iL);
ur = ur + AL/(1 + 2*j*alfa*z)*exp(-(k*alfa*r.^2 - (sin(Fir) + cos(Fir))*Ds.*r - j*Ds^2*z/k)/(1 + 2*j*alfa*z)); % In cylindrical coordinates
end
I_r = ur.*conj(ur);I_r = I_r / maxgen;
I_rarr = cat(3,I_rarr,I_r); 
end

figure(1)
%%%%%%%% Four Individual 3D Intensity Plots %%%%%%%%%%%%%%%%%
for iz = 1:length(zarr)
I_r = I_rarr(:,:,iz);subplot(2,2,iz);
rxset = rxmat(iz,:);ryset = rymat(iz,:);[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('\its_x\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',-12);
ylabel('\its_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 '];
Dsstr = [' \itD_s \rm\bf= ' num2str(Dsarr(ib,1),'%3.0f') ' m^-^1 '];
leg1 = [legtitle(ib,:) Dsstr];
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))*1.05]);
end
%figure(2)
%%%%%%%% Four Individual Contour Plots %%%%%%%%%%%%%%%%%
% for ib = 1:length(alfasarr)
% I_s = I_sarr(:,:,ib);subplot(2,2,ib);
% sxset = sxmat(ib,:);syset = symat(ib,:);[sx,sy]=meshgrid(sxset,syset);sxn = sx / 1e-2;syn = sy / 1e-2;
% [C,h] = contour(sxn,syn,I_s,'-k');axis square;fax = 1;
% axis ([min(min(sxn))*1 max(max(sxn))*1 min(min(syn))*1 max(max(syn))*1]);
% line([min(min(sxn))*fax max(max(sxn))*fax],[0 0],'LineWidth',1.8,'color','black','LineStyle','--');
% line([0 0],[min(min(syn))*fax max(max(syn))*fax],'LineWidth',1.8,'color','black','LineStyle','--');
% xlabel('\its_x\rm\bf axis in cm','FontSize',10,'FontWeight','bold');
% ylabel('\its_y\rm\bf axis in cm','FontSize',10,'FontWeight','bold');
% alfasstr = [' \it\alpha_s \rm\bf= ' num2str(alfasarr(ib)/1e-2,'%3.1f') ' cm '];
% Dsstr = [' \itD_s \rm\bf= ' num2str(Dsarr(ib,1),'%3.0f') ' m^-^1 '];
% leg1 = [legtitle(ib,:) Dsstr];
% 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