% 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.0e-6; k = 2*pi/lamda;
sxst = -1.e-1:0.2e-2:1e-1;syst = sxst;rxst = sxst;ryst = sxst;sig2 = 1e0;
sigstr = [' \it\sigma_s \rm\bf= ' num2str(sqrt(sig2)/1e-2,'%3.1f') ' cm '];
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.8;ryset2 = rxset2;rxset3 = rxst*1.5;ryset3 = rxset3;rxset4 = rxst*8;ryset4 = rxset4;
rxmat = [rxset1; rxset2; rxset3; rxset4];rymat = [ryset1; ryset2; ryset3; ryset4];
%%%%%%%%%%%%%%%%%%%%%%%%%%      Gaussian BEAM     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ALarr = [1 0;j -j;1 1;1 -1];Dmat = [0 25 50 200 200 250];
alfasxarr = [1e-2 1e-2;1e-2 1e-2;1e-2 1e-2;1e-2 0.8e-2];alfasyarr = [1e-2 1e-2;1e-2 1e-2;1e-2 1e-2;1e-2 0.8e-2];
alfasxarr = alfasxarr/sqrt(2);alfasyarr = alfasyarr/sqrt(2);
Dxarr = [j*Dmat(1) -j*Dmat(1);j*Dmat(4) -j*Dmat(4);1*Dmat(4) -1*Dmat(4);
          1*Dmat(1) -1*Dmat(1)];
Dyarr = [j*Dmat(4) -j*Dmat(4);j*Dmat(4) -j*Dmat(4);1*Dmat(4) -1*Dmat(4);
          1*Dmat(1) -1*Dmat(1)];     
Fsarr = [5e15 5e15;5e15 5e15;5e15 5e15;5e15 5e15];lamstr = [' \it\lambda = \rm\bf' num2str(lamda/1e-6) ' \mum'];
alfaxarr = 1./(k*alfasxarr.^2) + j./(2*Fsarr);alfayarr = 1./(k*alfasyarr.^2) + j./(2*Fsarr);Larr = [0.05e3 0.5e3 1e3 5e3];
alfasstr = [' \it\alpha_s \rm\bf= ' num2str(alfasxarr(1)/1e-2,'%3.1f') ' cm '];secstr = {sigstr;' ';lamstr};
%%%% Legends %%%%%%
legtitle = [' cos Gaussian beam         ';
            ' sine Gaussian beam        ';
            ' cosh Gaussian beam        ';
            ' sinh Gaussian beam        '];
%%%%%%%%%%%%% Source Beam %%%%%%%%%%%%%%%%%%
I_sarr = [];
ib = 4; % Select one beam
sxset = sxmat(ib,:);syset = symat(ib,:);[sx,sy]=meshgrid(sxset,syset);s1x = sx;s1y = sy;s2x = sx;s2y = sy;
I_s = zeros(size(sx));

for iL1 = 1:length(ALarr(1,:));for iL2 = 1:length(ALarr(1,:))
AL1 = ALarr(ib,iL1);AL2 = conj(ALarr(ib,iL2));
alfax1 = alfaxarr(ib,iL1);alfax2 = conj(alfaxarr(ib,iL2));alfay1 = alfayarr(ib,iL1);alfay2 = conj(alfayarr(ib,iL2));
Dx1 = Dxarr(ib,iL1);Dy1 = Dxarr(ib,iL1);Dx2 = conj(Dxarr(ib,iL2));Dy2 = conj(Dyarr(ib,iL2));
I_s = I_s + AL1*AL2*exp(-0.5*k*alfax1*s1x.^2 + Dx1*s1x -0.5*k*alfay1*s1y.^2 + Dy1*s1y).* ....
                    exp(-0.5*k*alfax2*s2x.^2 + Dx2*s2x -0.5*k*alfay2*s2y.^2 + Dy2*s2y); 
end;end
I_s = exp(-(s1x.^2 + s2x.^2 + s1y.^2 + s2y.^2 - 2*s1x.*s2x - 2*s1y.*s2y)/(2*sig2)).*I_s;maxgen = max(max(I_s));
%I_s = I_s / max(max(I_s));

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

I_rarr = [];
for iL = 1:length(Larr)     
L = Larr(iL);
rxset = rxmat(iL,:);ryset = rymat(iL,:);[rx,ry]= meshgrid(rxset,ryset);r1x = rx;r1y = ry;r2x = rx;r2y = ry;
I_r = zeros(size(rx));
for iL1 = 1:length(ALarr(1,:));for iL2 = 1:length(ALarr(1,:))
AL1 = ALarr(ib,iL1);AL2 = conj(ALarr(ib,iL2));
alfax1 = alfaxarr(ib,iL1);alfax2 = conj(alfaxarr(ib,iL2));alfay1 = alfayarr(ib,iL1);alfay2 = conj(alfayarr(ib,iL2));
Dx1 = Dxarr(ib,iL1);Dy1 = Dyarr(ib,iL1);Dx2 = conj(Dxarr(ib,iL2));Dy2 = conj(Dyarr(ib,iL2));
I_r = I_r + AL1*AL2/(k*alfax1*alfax2*sig2*L^2 + (alfax1 + alfax2)*L^2 + j*k*(alfax1 - alfax2)*sig2*L + k*sig2)^0.5/ ....
                    (k*alfay1*alfay2*sig2*L^2 + (alfay1 + alfay2)*L^2 + j*k*(alfay1 - alfay2)*sig2*L + k*sig2)^0.5* ....
                    exp(0.5/(L*(k*alfax1*sig2*L + L -j*k*sig2))*(sig2*(-j*k*r1x + Dx1*L).^2 + ....
(-j*(r1x - r2x)*k*L + j*(alfax1*L - j)*k^2*r2x*sig2 + Dx1*L^2 + (k*alfax1*sig2*L + L -j*k*sig2)*Dx2*L).^2/ ....
(k*(k*alfax1*alfax2*sig2*L^2 + (alfax1 + alfax2)*L^2 + j*k*(alfax1 - alfax2)*sig2*L + k*sig2)))).* ....
                    exp(0.5/(L*(k*alfay1*sig2*L + L -j*k*sig2))*(sig2*(-j*k*r1y + Dy1*L).^2 + ....
(-j*(r1y - r2y)*k*L + j*(alfay1*L - j)*k^2*r2y*sig2 + Dy1*L^2 + (k*alfax1*sig2*L + L -j*k*sig2)*Dy2*L).^2/ ....
(k*(k*alfay1*alfay2*sig2*L^2 + (alfay1 + alfay2)*L^2 + j*k*(alfay1 - alfay2)*sig2*L + k*sig2))));

 end;end
I_r = k*sig2*exp(j*k/(2*L)*(r1x.^2 - r2x.^2 + r1y.^2 - r2y.^2)).*I_r / maxgen;
I_rarr = cat(3,I_rarr,I_r);
end

figure(1)
%%%%%%%% Four Individual 3D Intensity Plots %%%%%%%%%%%%%%%%%
for iL = 1:length(Larr)
I_r = I_rarr(:,:,iL);subplot(2,2,iL);
rxset = rxmat(iL,:);ryset = rymat(iL,:);[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(alfasxarr(ib)/1e-2,'%3.1f') ' cm '];
Lstr = [' \itL \rm\bf= ' num2str(Larr(iL),'%3.0f') ' m '];
leg1 = [legtitle(ib,:) Lstr];
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]);
text(min(min(rxn))*1.5,min(min(ryn))*0.5,max(max(max(I_rarr)))*1,secstr,'FontSize',12,'FontWeight','bold','BackgroundColor','white');
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