% This is Dark Hollow beam in at source plane
% based on my source formulation dated 4.07.2006

function DHbeam_sour
clear;clc;clf reset;warning off MATLAB:divideByZero;close all
lamda = 1.55e-6;k = 2*pi/lamda;

%%%%%%%%%%%% Source Parameters %%%%%%%%%%%%%%%%
Fsxarr = [5e15 5e15];Fsyarr = [5e15 5e15];n = 0;m = 0;
Rarr = [1 2 5 10];Tarr = [1 1 1 1];
alfasxarr = [1e-2 0.8e-2];alfasyarr = [1e-2 0.8e-2]; A1 = 1;A2 = 1;Aarr = [A1 -A2];

alfaxarr = 1./(k*alfasxarr.^2) + j./Fsxarr;alfayarr = 1./(k*alfasyarr.^2) + j./Fsyarr;
alfax1 = alfaxarr(1);alfax2 = alfaxarr(2);alfay1 = alfayarr(1);alfay2 = alfayarr(2);
sxsetn = -1.e-1:0.2e-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.4;syset = sysetn*0.4;
% -------- Intensity Plot of General Hermite Gaussian Beam for ABCDGH System ------------
[sx,sy]=meshgrid(sxset,syset);sxn = sx / 1e-2;syn = sy / 1e-2;
ax =1./alfasxarr(1);ay = 1.0./alfasyarr(1);bx =0*ax; by = 0*ay;
sx1 = sx;sx2 = sx;sy1 = sy;sy2 = sy;
lamstr = [' \it\lambda = \rm\bf' num2str(lamda/1e-6) ' \mum'];

nstr = {' \itn \rm\bf= ' num2str(n,'%3.0f') ' '};mstr = {' \itm \rm\bf= ' num2str(m,'%3.0f') ' '};

alfastrsx = {' \it\alpha_s_x = \rm\bf' num2str(alfasxarr/1e-2,'%4.1f') ' cm '};
alfastrsy = {' \it\alpha_s_y = \rm\bf' num2str(alfasyarr/1e-2,'%4.1f') ' cm '};
Fsxstr = {' \itF_s_x = \rm\bf' num2str(Fsxarr,'%3.0f') ' m '};
Fsystr = {' \itF_s_y = \rm\bf' num2str(Fsyarr,'%3.0f') ' m '};
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')};

%%%%%%%%%%%%% Beam at Source Plane %%%%%%%%%%%%%%%%%%
I_sarr = [];
for ib = 1:length(Rarr)
    R = Rarr(ib);T = Tarr(ib);
    u_s = zeros(size(sx));
for r = 1:R
for tt = 1:T;for i = 1:2
    t = tt;if T == 1;t = r;end;Br = funb(R,r);Bt = 1;if T > 1;Bt = funb(T,t);end
%     u_sarr = (-1)^(r+t)*Br*Bt*(A1*exp(-0.5*k*(r*alfax1*sx.^2 + t*alfay1*sy.^2)) - ....
%     A2*exp(-0.5*k*(r*alfax2*sx.^2 + t*alfay2*sy.^2)));
    u_sarr = (-1)^(r+tt)*Br*Bt*Aarr(i)*exp(-0.5*k*(r*alfaxarr(i)*sx.^2 + t*alfayarr(i)*sy.^2));(-1)^(r+tt)*Br*Bt*Aarr(i)/(R*T);  
    %%% Alternative expression
u_s = u_s + u_sarr; 
end;end;end
u_s = u_s.*H(n,ax*sx+bx).*H(m,ay*sy+by)/(R*T); 
I_s = u_s.*conj(u_s);maxgen = max(max(I_s));I_s = I_s / maxgen;
I_sarr = cat(3,I_sarr,I_s);
end
figure(1)
%%%%%%%% Single Plot of First Beam %%%%%%%%%%%%%%%%%
for ib = 1:length(Rarr)
    I_s = I_sarr(:,:,ib);subplot(2,2,ib);
meshc(sxn,syn,real(I_s));view([35 45]);colormap([0 0 0]);view([15 30]) % For Fig 3
xlabel('\its_x\rm\bf axis in cm','FontSize',14,'FontWeight','bold','Rotation',-14);
ylabel('\its_y\rm\bf axis in cm','FontSize',14,'FontWeight','bold','Rotation',31)
xlabel('\its_x\rm\bf axis in cm','FontSize',14,'FontWeight','bold','Rotation',-4.5); % For Fig 3
ylabel('\its_y\rm\bf axis in cm','FontSize',14,'FontWeight','bold','Rotation',47)
zlabel('\fontname{Verdana}\fontsize{12}\itI_s_N','FontSize',12,'FontWeight','bold')

Rstr = [' \itR \rm\bf= ' num2str(Rarr(ib),'%2.0f') ' '];
Tstr = [' \itT \rm\bf= ' num2str(Tarr(ib),'%2.0f') ' '];
alfastrsx = [' \it\alpha_s_x = \rm\bf' num2str(alfasxarr(1,:)/1e-2,'%4.1f') ' cm '];
alfastrsy = [' \it\alpha_s_y = \rm\bf' num2str(alfasyarr(1,:)/1e-2,'%4.1f') ' cm '];

leg1 = [Rstr ',' Tstr ',' alfastrsx];
%title('Source intensity plot for a general beam with parameters','FontSize',10)
text(min(min(sxn))*0.9,min(min(syn))*0.1,max(max(I_s))*.9,leg1,'FontSize',14,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');
set(gcf,'Renderer','Zbuffer');set(gcf,'Color','White');set(gca,'FontSize',14);
axis ([min(min(sxn))*1 max(max(sxn))*1 min(min(syn))*1 max(max(syn))*1 0 1]);
end
figure(2)
%%%%%  Contour Plots  %%%%%%%%%%%%%%%%%%%%%%%%
for ib = 1:length(Rarr)
I_s = I_sarr(:,:,ib);subplot(2,2,ib);
[C,h] = contour(sxn,syn,I_s,':k');axis square
set(h,'LineWidth',2.5);set(gcf,'Color',[1 1 1])
xlabel('\its_x\rm\bf  axis in cm','FontSize',12,'FontWeight','bold');
ylabel('\its_y\rm\bf  axis in cm','FontSize',12,'FontWeight','bold');fax = 0.9;
axis([min(min(sxn))*fax max(max(sxn))*fax min(min(syn))*fax max(max(syn))*fax])
set(h,'LineWidth',2.0);set(gca,'FontSize',13);fax = 1.15;
Hc = clabel(C,h);set(Hc,'FontSize',10,'LineWidth',4);
Rstr = [' \itR \rm\bf= ' num2str(Rarr(ib),'%2.0f') ' '];
Tstr = [' \itT \rm\bf= ' num2str(Tarr(ib),'%2.0f') ' '];
leg1 = [Rstr ',' Tstr ',' alfastrsx];
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','--');
line([min(min(sxn))*fax max(max(sxn))*fax],[min(min(syn))*fax max(max(syn))*fax],'LineWidth',2.0,'color','black','LineStyle','-');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

text(min(min(sxn))*0.95,max(max(syn))*fax*0.5,leg1,'FontSize',11,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');

text(13,0,'Slanted axis','FontSize',12,'FontWeight','bold','Rotation',0);
end
%%%%%%% Functions Used %%%%%%%%%%%%%%%

function y = funb(a,b)
y = factorial(a)/factorial(b)/factorial(a-b);

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);