% This is Generalized Hermite hyperbolic and sinosoidal Gaussian beam for ABCD system
% based on four summation (dervation dated 07.07.2005). - 10.07.2005, 1.09.2005
% in the last one, focusing parameters Fx and Fy are incorporated.
% Note that both Errfunc with summation is used. 
% Note that in cat I_s is no longer the first page
% This programme takes into account one optical element with x - y asymmetry.
% Notation is changed accordingly - 9.09.2005
% Revisited on 6.12.2011, works OK, addition of Errfuncw is incomplete.
% With the present settings, produces Fig. 3 of Optik paper
function HergenABCDF1 
clear;clc;clf reset;warning off MATLAB:divideByZero;close all
lamda = 1.55e-6;k = 2*pi/lamda;%lamda = 0.633e-6;k = 2*pi/lamda;% For Andrews example
alfasx = 1e-2;alfasy = 1*alfasx;Fsx = 5e2;Fsy = 5e2;% For source size of 1 cm
alfaax = 4e-2;alfaay = alfaax;FLx = 10e-2;FLy = 10e-2;Ld1 = 0.55e3;Ld2 = 12e-2;% For cos beam
alfaax = 5e-2;alfaay = alfaax;FLx = 5e-2;FLy = 5e-2;Ld1 = 1.2e3;Ld2 = 6e-2;% For cosh beam
%alfaax = 1e-2;alfaay = 10e-2;FLx = 20e-2;FLy = 20e-2;Ld1 = 0.55e3;Ld2 = 25e-2;% For TEM beam
alfaLx = 1/(k*alfaax^2) + j/FLx;alfaLy = 1/(k*alfaay^2) + j/FLy;
na = [0 1 2 3] ;ma = [0 1 2 3];Vxarr = [0 4.7 25 50 150 250];Vyarr = Vxarr;
tx1 = -0.65*alfasx;tx2 = +0.65*alfasx;ty1 = -0.65*alfasy;ty2 = +0.65*alfasy;% For cos beam
tx1 = -10*alfasx;tx2 = +10*alfasx;ty1 = -10*alfasy;ty2 = +10*alfasy;% For cosh beam
%tx1 = -2*alfasx;tx2 = +2*alfasx;ty1 = -0.1*alfasy;ty2 = +0.1*alfasy;% For cross beam
tx1n = tx1 / 1e-2;tx2n = tx2 / 1e-2;ty1n = ty1 / 1e-2;ty2n = ty2 / 1e-2;
sxsetn = -1.e-1:0.3e-2:1e-1;sysetn = sxsetn;
%sxsetn = [-10e-2 -5e-2 5e-2 10e-2];sysetn = sxsetn;
sxset = sxsetn*1.0;syset = sysetn*1.0;
sxset = sxsetn*0.2;syset = sysetn*0.2;% This setting is for cos beam 
sxset = sxsetn*0.5;syset = sysetn*0.5; % This setting is for cosh beam
%sxset = sxsetn*0.35;syset = sysetn*0.35;% This setting is for TEM beam 
%sxset = sxset*0.1;syset = syset*0.1;% This setting is for source size of 0.2 cm 
%%%%%%%%%%%%%%%%%%%%%%%  Generalized Hermite sinosiodal / hyperbolic beam  %%%%%%%%%%%%%%%%%%%%%%%%%%
% -------- Intensity Plot of General Hermite Gaussian Beam for ABCD System ------------
[sx,sy]=meshgrid(sxset,syset);
%R = sqrt(sxset.^2+syset.^2);RS = [-R(1:round(numel(R)/2)) R(round(numel(R)/2)+1:end)];
n = na(1);m = ma(1);Anm = 1;% For cos beam
n = na(2);m = ma(1);Anm = 1;% For cosh beam
%%%%%%%%%%% Setting Beam Type %%%%%%%%%%%%%%%%%%%%
L1 = 0.5; L2 = 0.5;Vx = Vxarr(6);Vy = Vx;Yx = -Vx;Yy = -Vy; % cos Beam
L1 = 0.5; L2 = 0.5;Vx = -j*Vxarr(5);Vy = Vx;Yx = -Vx;Yy = -Vy; %cosh Beam
%L1 = -0.5*j; L2 = 0.5*j;Vx = Vxarr(5);Vy = Vx;Yx = -Vx;Yy = -Vy; %sine Beam
L1 = 0.5; L2 = -0.5;Vx = -j*Vxarr(5);Vy = Vx;Yx = -Vx;Yy = -Vy; %Sinh Beam
%L1 = 0.5; L2 = 0.5;Vx = Vxarr(1);Vy = Vx;Yx = -Vx;Yy = -Vy; % For TEM or cross beam
lamstr = [' \it\lambda = \rm\bf' num2str(lamda/1e-6) ' \mum'];
ax =1/alfasx;ay = 1.00/alfasy;bx =0*ax; by = 0*ay;%ax =ax/sqrt(2);ay = ay/sqrt(2);
nstr = [' \itn \rm\bf= ' num2str(n,'%3.0f') ' '];mstr = [' \itm \rm\bf= ' num2str(m,'%3.0f') ' '];
Vxstr = [' \itV_x = \rm\bf' num2str(abs(Vx)) ' m^-^1 '];Vystr = [' \itV_y = \rm\bf' num2str(abs(Vy)) ' m^-^1 '];
alfastrsx = [' \it\alpha_s_x = \rm\bf' num2str(alfasx/1e-2,'%3.1f') ' cm '];
alfastrsy = [' \it\alpha_s_y = \rm\bf' num2str(alfasy/1e-2,'%3.1f') ' cm '];
alfaaxstr = [' \it\alpha_a_x = \rm\bf' num2str(alfaax/1e-2,'%3.1f') ' cm '];
alfaaystr = [' \it\alpha_a_y = \rm\bf' num2str(alfaay/1e-2,'%3.1f') ' cm '];
Fsxstr = [' \itF_s_x = \rm\bf' num2str(Fsx,'%3.0f') ' m '];
Fsystr = [' \itF_s_y = \rm\bf' num2str(Fsy,'%3.0f') ' m '];
FLxstr = [' \itF_l_x = \rm\bf' num2str(FLx/1e-2,'%3.0f') ' cm '];
FLystr = [' \itF_l_y = \rm\bf' num2str(FLy/1e-2,'%3.0f') ' cm '];
axstr = [' \it\fontname{century gothic}a\fontname{helvetica}_x = \rm\bf' num2str(ax,'%3.1f') ' m^-^1 '];
aystr = [' \it\fontname{century gothic}a\fontname{helvetica}_x = \rm\bf' num2str(ay,'%3.1f') ' m^-^1 '];
bxstr = [' \itb_x = \rm\bf' num2str(bx,'%3.1f')];
bystr = [' \itb_x = \rm\bf' num2str(by,'%3.1f')];
Ld1str = [' \itL\rm\bf_1 = ' num2str(Ld1/1e3) ' km'];
Ld2str = [' \itL\rm\bf_2 = ' num2str(Ld2/1e-2) ' cm'];
%%%%%%%%%%%%% Beam at Source Plane %%%%%%%%%%%%%%%%%%

leg1 = {[nstr ',' mstr ',' alfastrsx ',' alfastrsy ',' Vxstr ',' Vystr];[axstr ',' aystr ',' bxstr ',' bystr]};
u_s = L1*Anm*(H(n,ax*sx+bx).*H(m,ay*sy+by)).*exp(-(sx.^2/(2*alfasx^2)+sy.^2/(2*alfasy^2))).*exp(j*Vx*sx+j*Vy*sy) + ....
    L2*Anm*(H(n,ax*sx+bx).*H(m,ay*sy+by)).*exp(-(sx.^2/(2*alfasx^2)+sy.^2/(2*alfasy^2))).*exp(j*Yx*sx+j*Yy*sy);

I_s = u_s.*conj(u_s);maxgen = max(max(I_s));I_s = I_s / maxgen;
sxn = sx / 1e-2;syn = sy / 1e-2;

%%%%%%%%%%%%%%%%% Intensity Profile at Receiver Plane %%%%%%%%%%%
% Beam after Propogation
% Arrangement of the elements in ABCD matrix
Axarr = [1 1 1 + j*alfaLx*Ld2 1];Ayarr = [1 1 1 + j*alfaLy*Ld2 1];
Bxarr = [Ld1 Ld1 Ld1 + Ld2*(1 + j*alfaLx*Ld1) 0.1];Byarr = [Ld1 Ld1 Ld1 + Ld2*(1 + j*alfaLy*Ld1) 0.1];
Cxarr = [0 j*alfaLx j*alfaLx 0];Cyarr = [0 j*alfaLy j*alfaLy 0];
Dxarr = [1 1 + j*alfaLx*Ld1 1 + j*alfaLx*Ld1 1];Dyarr = [1 1 + j*alfaLy*Ld1 1 + j*alfaLy*Ld1 1];
pxsetarr = [7*sxset;3.5*sxset;0.8*sxset;0.6*sxset];pysetarr = pxsetarr;% For cos beam
%pxsetarr = [8*sxset;3.5*sxset;1*sxset;0.8*sxset];pysetarr = pxsetarr;% For TEM beam
%%% For cosh beam %%%%%%
Axarr = [1 1 1 + j*alfaLx*Ld2];Ayarr = [1 1 1 + j*alfaLy*Ld2];
Bxarr = [Ld1 Ld1 Ld1 + Ld2*(1 + j*alfaLx*Ld1)];Byarr = [Ld1 Ld1 Ld1 + Ld2*(1 + j*alfaLy*Ld1)];
Cxarr = [0 j*alfaLx j*alfaLx 0];Cyarr = [0 j*alfaLy j*alfaLy 0];
Dxarr = [1 1 + j*alfaLx*Ld1 1 + j*alfaLx*Ld1];Dyarr = [1 1 + j*alfaLy*Ld1 1 + j*alfaLy*Ld1];
pxsetarr = [2.2*sxset;1.8*sxset;0.35*sxset];pysetarr = pxsetarr;% For cosh beam
I_arr = [];%I_arr = cat(3,I_arr,I_s);
for imat = 1:length(Axarr)
    Ax = Axarr(imat);Ay = Ayarr(imat);Bx = Bxarr(imat);By = Byarr(imat);Dx = Dxarr(imat);Dy = Dyarr(imat);
    QAx = Bx*Fsx + j*k*alfasx^2*(Bx - Ax*Fsx);QAy = By*Fsy + j*k*alfasy^2*(By - Ay*Fsy);% According to new convention
    pxset = pxsetarr(imat,:);pyset = pysetarr(imat,:);
    [px,py]=meshgrid(pxset,pyset);pxn = px / 1e-2;pyn = py / 1e-2;
%%%%%%%%%%%%%%%% Summation over x %%%%%%%%%%%%%%%%%%%%%%%
VYarr = [Vx Yx];sumxa = [];
for iVY = 1:length(VYarr)
    VY = VYarr(iVY);sumx = 0; 
    QVY = Bx*VY - k*px;% According to new convention
    [row, col] = size(px);
gx1 = [];gx2 = [];
for irow = 1:row
    gx1r = [];gx2r = [];
    for icol = 1:col
        pxs = px(irow,icol);
        QVYs = Bx*VY - k*pxs;% According to new convention
        qoverp = j*Fsx*alfasx^2*QVYs/QAx;
        if qoverp <= tx1 & qoverp < tx2 & tx1 < tx2 ; gx1s = 1;gx2s = 1;end
        if tx1 < qoverp & tx1 < tx2 & qoverp < tx2 ; gx1s = -1;gx2s = 1;end
        if tx1 < tx2 & tx1 < qoverp & tx2 <= qoverp ; gx1s = -1;gx2s = -1;end
        gx1r = [gx1r gx1s];gx2r = [gx2r gx2s];
    end
    gx1 = [gx1;gx1r];gx2 = [gx2;gx2r];
end 
for Lx1 = 0:fix(n/2)
for Lnx1 = 0:n-2*Lx1
for kx1 = 0:Lnx1
    sume = -sqrt(pi)*(Errfuncs(sqrt(0.5/(Bx*Fsx*alfasx^2*QAx))*(tx1*QAx-j*Fsx*alfasx^2*QVY)) - ....
                      Errfuncs(sqrt(0.5/(Bx*Fsx*alfasx^2*QAx))*(tx2*QAx-j*Fsx*alfasx^2*QVY)));
    if kx1 > 0
    sume = 0;
    for kx2 = 0:fix((kx1-1)/2)    
    termexp = gx1.^(kx1-1).*exp(-0.5/(Bx*Fsx*alfasx^2*QAx)*(tx1*QAx-j*Fsx*alfasx^2*QVY).^2).* ....
                                (0.5/(Bx*Fsx*alfasx^2*QAx)*(tx1*QAx-j*Fsx*alfasx^2*QVY).^2).^((kx1-1)/2-kx2) - ....
              gx2.^(kx1-1).*exp(-0.5/(Bx*Fsx*alfasx^2*QAx)*(tx2*QAx-j*Fsx*alfasx^2*QVY).^2).* ....
                                (0.5/(Bx*Fsx*alfasx^2*QAx)*(tx2*QAx-j*Fsx*alfasx^2*QVY).^2).^((kx1-1)/2-kx2);   
    sume = sume + funE(kx1-1)/funE(kx1-1-2*kx2)/2^kx2*termexp;
    end
    termrem = sqrt(pi)*funE(kx1-1)*fix(abs(rem(kx1-1,2)))*(Errfuncs(sqrt(0.5/(Bx*Fsx*alfasx^2*QAx))*(tx1*QAx-j*Fsx*alfasx^2*QVY)) - ....
                                                           Errfuncs(sqrt(0.5/(Bx*Fsx*alfasx^2*QAx))*(tx2*QAx-j*Fsx*alfasx^2*QVY)));
    sume = sume - termrem;
    end         
    
TLx1 = funT(Lx1);facLx1 = funb(n,2*Lx1);facLnx1 = funb(n-2*Lx1,Lnx1);
sumx = sumx + (-1)^Lx1*2^(n-Lx1)*TLx1*facLx1*facLnx1*ax^Lnx1*bx^(n-2*Lx1-Lnx1)* ....
              (j*Fsx*alfasx^2*QVY/QAx).^Lnx1*0.5*factorial(Lnx1)/sqrt(0.5*QAx/(Bx*Fsx*alfasx^2)).* ....
              (0.5*j*QVY/Bx).^(-kx1).*(0.5*QAx/(Bx*Fsx*alfasx^2))^(kx1/2)/(factorial(kx1)*factorial(Lnx1-kx1)).*sume;

end 
end
end
sumxa = cat(3,sumxa,sumx);
end
%%%%%%%%%%%%%%%% Summation over y %%%%%%%%%%%%%%%%%%%%%%%
VYarr = [Vy Yy];sumya = [];
for iVY = 1:length(VYarr)
    VY = VYarr(iVY);sumy = 0; 
    QVY = By*VY - k*py;% According to new convention
    [row, col] = size(py);
gy1 = [];gy2 = [];
for irow = 1:row
    gy1r = [];gy2r = [];
    for icol = 1:col
        pys = py(irow,icol);
        QVYs = By*VY - k*pys;% According to new convention
        qoverp = j*Fsy*alfasy^2*QVYs/QAy;
        if qoverp <= ty1 & qoverp < ty2 & ty1 < ty2 ; gy1s = 1;gy2s = 1;end
        if ty1 < qoverp & ty1 < ty2 & qoverp < ty2 ; gy1s = -1;gy2s = 1;end
        if ty1 < ty2 & ty1 < qoverp & ty2 <= qoverp ; gy1s = -1;gy2s = -1;end
        gy1r = [gy1r gy1s];gy2r = [gy2r gy2s];
    end
    gy1 = [gy1;gy1r];gy2 = [gy2;gy2r];
end          
for Ly1 = 0:fix(m/2)
for Lmy1 = 0:m-2*Ly1
for ky1 = 0:Lmy1
    
    sume = -sqrt(pi)*(Errfuncs(sqrt(0.5/(By*Fsy*alfasy^2*QAy))*(ty1*QAy-j*Fsy*alfasy^2*QVY)) - ....
                      Errfuncs(sqrt(0.5/(By*Fsy*alfasy^2*QAy))*(ty2*QAy-j*Fsy*alfasy^2*QVY)));
    if ky1 > 0
    sume = 0;
    for ky2 = 0:fix((ky1-1)/2)    
    termexp = gy1.^(ky1-1).*exp(-0.5/(By*Fsy*alfasy^2*QAy)*(ty1*QAy-j*Fsy*alfasy^2*QVY).^2).* ....
                                (0.5/(By*Fsy*alfasy^2*QAy)*(ty1*QAy-j*Fsy*alfasy^2*QVY).^2).^((ky1-1)/2-ky2) - ....
              gy2.^(ky1-1).*exp(-0.5/(By*Fsy*alfasy^2*QAy)*(ty2*QAy-j*Fsy*alfasy^2*QVY).^2).* ....
                                (0.5/(By*Fsy*alfasy^2*QAy)*(ty2*QAy-j*Fsy*alfasy^2*QVY).^2).^((ky1-1)/2-ky2);   
    sume = sume + funE(ky1-1)/funE(ky1-1-2*ky2)/2^ky2*termexp;
    end
    termrem = sqrt(pi)*funE(ky1-1)*fix(abs(rem(ky1-1,2)))*(Errfuncs(sqrt(0.5/(By*Fsy*alfasy^2*QAy))*(ty1*QAy-j*Fsy*alfasy^2*QVY)) - ....
                                                           Errfuncs(sqrt(0.5/(By*Fsy*alfasy^2*QAy))*(ty2*QAy-j*Fsy*alfasy^2*QVY)));
    sume = sume - termrem;
    end         
    
TLy1 = funT(Ly1);facLy1 = funb(m,2*Ly1);facLmy1 = funb(m-2*Ly1,Lmy1);
sumy = sumy + (-1)^Ly1*2^(m-Ly1)*TLy1*facLy1*facLmy1*ay^Lmy1*by^(m-2*Ly1-Lmy1)* ....
              (j*Fsy*alfasy^2*QVY/QAy).^Lmy1*0.5*factorial(Lmy1)/sqrt(0.5*QAy/(By*Fsy*alfasy^2)).* ....
              (0.5*j*QVY/By).^(-ky1).*(0.5*QAy/(By*Fsy*alfasy^2))^(ky1/2)/(factorial(ky1)*factorial(Lmy1-ky1)).*sume;
end 
end
end
sumya = cat(3,sumya,sumy);
end
%sumxa
%sumya
QVx = Bx*Vx - k*px;
QVy = By*Vy - k*py;
QYx = Bx*Yx - k*px;
QYy = By*Yy - k*py; % According to new convention
expL1x = exp(-0.5*Fsx*alfasx^2*QVx.^2/(Bx*QAx));
expL1y = exp(-0.5*Fsy*alfasy^2*QVy.^2/(By*QAy));
expL2x = exp(-0.5*Fsx*alfasx^2*QYx.^2/(Bx*QAx));
expL2y = exp(-0.5*Fsy*alfasy^2*QYy.^2/(By*QAy));
u_r = -j*k/(2*pi*sqrt(Bx*By))*exp(+0.5*j*k*(Dx*px.^2/Bx + Dy*py.^2/By)).* ....
                   (L1*expL1x.*expL1y.*sumxa(:,:,1).*sumya(:,:,1) + ....
                    L2*expL2x.*expL2y.*sumxa(:,:,2).*sumya(:,:,2));                
I_r = u_r.*conj(u_r);I_r = I_r / maxgen; 
I_arr = cat(3,I_arr,I_r);
end
%I_arr(:,:,3)
%toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I_r = I_arr(:,:,3);L = Ld1 + Ld2; 
R = sqrt(sxset.^2+syset.^2);RP = [-R(1:round(numel(R)/2)) R(round(numel(R)/2)+1:end)];RPn = RP / 1e-2;
Lstr = [' \itL = \rm\bf' num2str(L/1e3) ' km'];

%%%%%%%% Combined Intensity Profile before and after Propagation %%%%%%%%%%%%%%%%%
plot(RPn,diag((I_s)),'-k','LineWidth',2); hold on
plot(RPn,diag((I_r)),'--k','LineWidth',3); 
xlabel('Slanted axis in cm','FontSize',12,'FontWeight','bold');
ylabel('\fontname{Times}< \it\bfI_r_N\fontname{Times}\rm > \bf','FontSize',14,'FontWeight','bold')
axis ([min(RPn)*1.2 max(RPn)*1.2 0 1.1])
set(gca,'FontSize',13);set(gcf,'Color','White')
leg1 = {[nstr ',' mstr];' ';[alfastrsx ',' alfastrsy];' ';[Fsxstr ',' Fsystr];' ';[Vxstr ',' Vystr]};
leg2 = {['                \itL \rm\bf= 0 km '];' ';['                \itL \rm\bf= ' num2str(fix(L/1e3)) ' km ']};
text(max(RPn)*0.9,max(max(I_s))*.9,leg1,'FontSize',12,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');
text(min(RPn)*0.9,max(max(I_s))*.9,leg2,'FontSize',12,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');
hold off;pause;clf reset
%%%%% Generalized Beam Contour Plot before and after Propagation %%%%%%%%%%%%%%%%%%%%%%%%
I_r = I_arr(:,:,2);L = Ld1 + Ld2;[px,py]=meshgrid(pxsetarr(2,:),pysetarr(2,:));
pxn = px / 1e-2;pyn = py / 1e-2;
Lstr = [' \itL = \rm\bf' num2str(L/1e3) ' km'];
[C,h] = contour(sxn,syn,I_s,':k');axis square
set(h,'LineWidth',1.5);set(gcf,'Color',[1 1 1]);hold on
[C,h] = contour(pxn,pyn,I_r,'-.k'); %%% Fig 8
xlabel('\its_x\rm\bf and \itp_x\rm\bf axis in cm','FontSize',12,'FontWeight','bold');
ylabel('\its_y\rm\bf and \itp_y\rm\bf axis in cm','FontSize',12,'FontWeight','bold');fax = 0.95;
axis([min(min(pxn))*fax max(max(pxn))*fax min(min(pyn))*fax max(max(pyn))*fax])
set(h,'LineWidth',2.5);Hc = clabel(C,h);set(Hc,'FontSize',10,'LineWidth',4);set(gca,'FontSize',13);
leg1 = {[nstr ',' mstr ',' axstr ',' bxstr ',' alfastrsx ',' Vxstr];' ';[Lstr ' ' lamstr]};
leg2 = 'Contour plot of \fontname{Times}\fontsize{14}\itI_s_N';
leg3 = 'Contour plot of \fontname{Times}\fontsize{14}\itI_r_N';
text(min(min(pxn))*0.8,max(max(pyn))*0.1,leg2,'FontSize',14,'FontWeight','bold','Rotation',0);
text(max(max(pxn))*0.15,min(min(pyn))*0.65,leg3,'FontSize',14,'FontWeight','bold','Rotation',0);
text(13,0,'Slanted axis','FontSize',12,'FontWeight','bold','Rotation',0);fax = 1.15;
if abs(tx2) < 5*alfasx
    %%%%%%%% Aperture Lines %%%%%%%%%%%%%%
    line([tx1n tx2n],[ty1n ty1n],'LineWidth',6.0,'color','black','LineStyle','--');
    line([tx1n tx2n],[ty2n ty2n],'LineWidth',6.0,'color','black','LineStyle','--');
    line([tx1n tx1n],[ty1n ty2n],'LineWidth',6.0,'color','black','LineStyle','--');
    line([tx2n tx2n],[ty1n ty2n],'LineWidth',6.0,'color','black','LineStyle','--');
    %%%%%%%% Coordinate Lines %%%%%%%%%%%%%
    %line([min(min(pxn))*fax tx1n],[ty1n ty1n],'LineWidth',1.5,'color','black','LineStyle','--');
    %line([min(min(pxn))*fax tx1n],[ty2n ty2n],'LineWidth',1.5,'color','black','LineStyle','--');
    line([tx2n max(max(pxn))*fax],[ty1n ty1n],'LineWidth',1.5,'color','black','LineStyle','--');
    line([tx2n max(max(pxn))*fax],[ty2n ty2n],'LineWidth',1.5,'color','black','LineStyle','--');
    line([tx1n tx1n],[min(min(pyn))*fax ty1n],'LineWidth',1.5,'color','black','LineStyle','--');
    line([tx2n tx2n],[min(min(pyn))*fax ty1n],'LineWidth',1.5,'color','black','LineStyle','--');
    %line([tx1n tx1n],[ty2n max(max(pyn))*fax],'LineWidth',1.5,'color','black','LineStyle','--');
    %line([tx2n tx2n],[ty2n max(max(pyn))*fax],'LineWidth',1.5,'color','black','LineStyle','--');
    text(tx1n, min(min(pyn))*fax,'\it\bft_x\rm\bf_1','FontSize',12,'FontWeight','bold','BackgroundColor','white');
    text(tx2n, min(min(pyn))*fax,'\it\bft_x\rm\bf_2','FontSize',12,'FontWeight','bold','BackgroundColor','white');
    text(max(max(pxn))*fax, ty1n,'\it\bft_y\rm\bf_1','FontSize',12,'FontWeight','bold','BackgroundColor','white');
    text(max(max(pxn))*fax, ty2n,'\it\bft_y\rm\bf_2','FontSize',12,'FontWeight','bold','BackgroundColor','white');
end
%%%%%%% s_x, s_y, Slanted Axis Lines %%%%%%%%%
line([min(min(pxn))*fax max(max(pxn))*fax],[0 0],'LineWidth',1.8,'color','black','LineStyle','--');
line([0 0],[min(min(pyn))*fax max(max(pyn))*fax],'LineWidth',1.8,'color','black','LineStyle','--');
line([min(min(pxn))*fax max(max(pxn))*fax],[min(min(pyn))*fax max(max(pyn))*fax],'LineWidth',2.0,'color','black','LineStyle','-');
text(min(min(pxn))*0.95,max(max(pyn))*fax*0.9,leg1,'FontSize',12,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');
pause;clf
%%%%%%%% Individual Intensity Profile before and after Propagation %%%%%%%%%%%%%%%%%
subplot(2,2,1);meshc(sxn,syn,real(I_s));view([35 45]);colormap([0 0 0])
xlabel('\its_x\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',-14);
ylabel('\its_y\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',35)
zlabel('\fontname{Verdana}\fontsize{12}\itI_s_N','FontSize',12,'FontWeight','bold')
title('Intensity at source plane','FontSize',10)
Lstr = [' \itL = \rm\bf' num2str(0/1e3) ' km'];
leg1 = {[nstr ',' mstr ',' alfastrsx ',' Fsxstr];' ';[Vxstr ',' Lstr]};
text(min(min(sxn))*0.9,min(min(syn))*0.1,max(max(I_s))*.9,leg1,'FontSize',10,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');
set(gcf,'Renderer','Zbuffer');set(gcf,'Color','White');set(gca,'FontSize',11);
zdat = 0.5;lgen = 3.0;lsty = '--';% For cos beam
%zdat = 1.05;lgen = 1.0;lsty = '-';% For TEM beam
% line([tx1n tx2n],[ty1n ty1n],[zdat zdat],'LineWidth',lgen,'color','black','LineStyle',lsty);% For cos and TEM beam
% line([tx1n tx2n],[ty2n ty2n],[zdat zdat],'LineWidth',lgen,'color','black','LineStyle',lsty);% For cos and TEM beam
% line([tx1n tx1n],[ty1n ty2n],[zdat zdat],'LineWidth',lgen,'color','black','LineStyle',lsty);% For cos and TEM beam
% line([tx2n tx2n],[ty1n ty2n],[zdat zdat],'LineWidth',lgen,'color','black','LineStyle',lsty);% For cos and TEM beam
axis ([min(min(sxn))*1 max(max(sxn))*1 min(min(syn))*1 max(max(syn))*1 0 1]);
%%%%% Second Plot %%%%%%%%%
%I_r = I_arr(:,:,4);[px,py]=meshgrid(pxsetarr(4,:),pysetarr(4,:));pxn = px/1e-2;pyn = py/1e-2;L = 0.1;% For cos beam
I_r = I_arr(:,:,1);[px,py]=meshgrid(pxsetarr(1,:),pysetarr(1,:));pxn = px/1e-2;pyn = py/1e-2;L = Ld1;% For cosh beam
Lstr = [' \itL = \rm\bf' num2str(L/1e3) ' km'];lamstr = [' \it\lambda = \rm\bf' num2str(lamda/1e-6) ' \mum'];
subplot(2,2,2);meshc(pxn,pyn,real(I_r));view([35 45]);
xlabel('\itp_x\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',-14);
ylabel('\itp_y\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',35)
zlabel('\fontname{Verdana}\fontsize{12}\it I_r_N ','FontSize',12,'FontWeight','bold')
title('Intensity before the lens','FontSize',10)
leg1 = Ld1str;
text(min(min(pxn))*0.9,min(min(pyn))*0.1,max(max(I_r))*.9,leg1,'FontSize',10,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');
set(gcf,'Renderer','Zbuffer');set(gcf,'Color','White');set(gca,'FontSize',11);
axis ([min(min(pxn))*1 max(max(pxn))*1 min(min(pyn))*1 max(max(pyn))*1 0 max(max(I_r))]);
%%%%% Third Plot %%%%%%%%%
I_r = I_arr(:,:,1);[px,py]=meshgrid(pxsetarr(1,:),pysetarr(1,:));pxn = px/1e-2;pyn = py/1e-2;L = Ld1;% For cos beam
I_r = I_arr(:,:,2);[px,py]=meshgrid(pxsetarr(2,:),pysetarr(2,:));pxn = px/1e-2;pyn = py/1e-2;L = Ld1;% For cosh beam
Lstr = [' \itL = \rm\bf' num2str(L/1e3) ' km'];lamstr = [' \it\lambda = \rm\bf' num2str(lamda/1e-6) ' \mum'];
subplot(2,2,3);meshc(pxn,pyn,real(I_r));view([35 45]);
xlabel('\itp_x\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',-14);
ylabel('\itp_y\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',35)
zlabel('\fontname{Verdana}\fontsize{12}\it I_r_N ','FontSize',12,'FontWeight','bold')
title('Intensity after the lens','FontSize',10)
leg1 = [alfaaxstr ',' FLxstr];
text(min(min(pxn))*0.7,min(min(pyn))*0.1,max(max(I_r))*.9,leg1,'FontSize',10,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');
set(gcf,'Renderer','Zbuffer');set(gcf,'Color','White');set(gca,'FontSize',11);
axis ([min(min(pxn))*1 max(max(pxn))*1 min(min(pyn))*1 max(max(pyn))*1 0 max(max(I_r))]);
%%%%% Fourth Plot %%%%%%%%%
I_r = I_arr(:,:,3);[px,py]=meshgrid(pxsetarr(3,:),pysetarr(3,:));pxn = px/1e-2;pyn = py/1e-2;L = Ld1 + Ld2;
subplot(2,2,4);meshc(pxn,pyn,real(I_r));view([35 45]);
xlabel('\itp_x\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',-14);
ylabel('\itp_y\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',35)
zlabel('\fontname{Verdana}\fontsize{12}\it I_r_N ','FontSize',12,'FontWeight','bold')
title('Intensity at receiver plane','FontSize',10)
leg1 = [Ld1str ',' Ld2str];
text(min(min(pxn))*0.9,min(min(pyn))*0.1,max(max(I_r))*.9,leg1,'FontSize',10,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');
set(gcf,'Renderer','Zbuffer');set(gcf,'Color','White');set(gca,'FontSize',11);
axis ([min(min(pxn))*1 max(max(pxn))*1 min(min(pyn))*1 max(max(pyn))*1 0 max(max(I_r))]);
pause;clf reset
%%%%%%% Functions Used %%%%%%%%%%%%%%%
function y = funb(a,b)
y = factorial(a)/factorial(b)/factorial(a-b);
function y = funT(el)
T = 1;
for j = 1:el
    T = (2*j-1)*T;
end
y = T;
function y = funE(ez)
if rem(ez,2) == 0;y = 2:2:ez;end
if rem(ez,2) == 1;y = 3:2:ez;end
y = prod(y);
function yout = H(nn,xx)
mat = cat(4,ones(size(xx)),2*xx,4*xx.^2-2);
for say = 4:nn+1
    mat(:,:,:,say) = 2*xx.*mat(:,:,:,say-1) - 2*(say-2)*mat(:,:,:,say-2);
end
yout = mat(:,:,:,nn+1);
%%%%%%%%%%% Calculation of error fun via quade %%%%%%%%%%%%%%%%%
function yerr = Errfunc(zarg)
%clear global zarg
global zargs
[row, col] = size(zarg);yerr = [];
for irow = 1:row
    yerrr = [];
    for icol = 1:col
        zargs = zarg(irow,icol);
        yerrs = 2*zargs/sqrt(pi)*quade(@funerrc,0,1);
        yerrr = [yerrr yerrs];
    end
    yerr = [yerr;yerrr];
end
%zarg
%yerr
function y = funerrc(x) 
global zargs
y = exp(-zargs^2*x.^2);
%%%%%%%%%%% Calculation of error fun via summation %%%%%%%%%%%%%%%%%
function yerr = Errfuncs(zarg)
x = real(zarg);y = imag(zarg);
res1 = erf(x) + exp(-x.^2).*(1 - cos(2*x.*y) + j*sin(2*x.*y))./(2*pi*x); 
fnxy = 2*x - 2*x.*cosh(y).*cos(2*x.*y) + sinh(y).*sin(2*x.*y);
gnxy = 2*x.*cosh(y).*sin(2*x.*y) + sinh(y).*cos(2*x.*y);
term = exp(-0.25)*(fnxy + j*gnxy)./(1 + 4*x.^2);n = 1;res2 = term;
    while abs(max(max(term))) > 1e-15
    n = n + 1;
    fnxy = 2*x - 2*x.*cosh(n*y).*cos(2*x.*y) + n*sinh(n*y).*sin(2*x.*y);
    gnxy = 2*x.*cosh(n*y).*sin(2*x.*y) + n*sinh(n*y).*cos(2*x.*y);
    term = exp(-0.25*n^2)*(fnxy + j*gnxy)./(n^2 + 4*x.^2);
    res2 = res2 + term;
    end
yerr = res1 + 2*exp(-x.^2).*res2/pi;
%%%%% Calculation of error function according to Weideman SIAM 1994, Table 1 %%%%%
function yerr = Errfuncw(zarg)
global Nsum Linc Ancoef
z1 = j*zarg;sig = sign(real(zarg));zc = z1.*sig;
Z = (Linc + j*zc)./(Linc - j*zc); ply = polyval(Ancoef,Z);
wfun = 2*ply./(Linc - j*zc).^2 + (1/sqrt(pi))./(Linc - j*zc);
wfun = 0.5*(1-sig)*2.*exp(-z1.^2) + sig.*wfun;
yerr = 1 - wfun.*exp(z1.^2);