% Diffraction pattern on the detector in Q-2Q geometry 
% from planar tapered waveguide illuminated by elementary 
% radiator with parabolic wave front. (Distribution of 
% intensity is a composite of diffraction patterns on the 
% detector corresponding to measurements with different 
% grazing incidence angles of the incoming beam).   
%
clear 
global k ri Ri DD Q0 m Alpha Din L Delta Beta
%tic

%---------------------Experiment ANKA--------------------------------------
lambda=0.113;                            % wavelength (nm) 
delta = 4.0336E-06;  beta = 5.0448E-08;  % Si n=1-delta+i*beta
L=4.e6;                                  % length of WG (nm)
Dex=1600;                                % exit aperture (nm) EXPER.1
Alpha=3.1e-4;                            % angle of the tilt (rad) EXPER.1

%--------------------Experiment ELETTRA------------------------------------
lambda=0.161;                            % wavelength (nm) 
delta = 8.2856E-06;  beta = 2.0523E-07;  % Si n=1-delta+i*beta
L=4.6e6;                                 % length of WG (nm)
Dex=430;                                 % exit aperture (nm)  
Alpha=1.e-7;                             % angle of the tilt (rad) 

%--------------------Experiment ESRF---------------------------------------
%lambda=0.113;                            % wavelength (nm)
%delta = 4.0336E-06;  beta = 5.0448E-08;  % Si n=1-delta+i*beta
%L=4e6;                                   % length of WG (nm)
%Dex=670;                                 % exit aperture (nm)
%Alpha=1.3e-4;                            % angle of the tilt (rad)

%--------------------Experiment Lab.source1--------------------------------
%lambda=0.154;                            % wavelength (nm)
%delta = 7.5735E-06;  beta = 1.7245E-07;  % wavelength (nm) 
%L=4.5e6;                                 % length Lab.source1
%Dex=330;                                 % exit aperture (nm) 
%Alpha=4e-5;                              % angle of the tilt (rad)

%--------------------Experiment Lab.source2--------------------------------
%lambda=0.154;                            % wavelength (nm)
%delta = 7.5735E-06;  beta = 1.7245E-07;  % wavelength (nm) 
%L=4.9e6;                                 % length Lab.source1
%Dex=200;                                 % exit aperture (nm) 
%Alpha=-2.16e-5;                          % angle of the tilt (rad)
%--------------------------------------------------------------------------                               

Ri=5e6;                     %distance Source - WG (nm) 
N0=100;                     %number of measurements on grazing incid.angles 
Nfi=500;                    %step on Fi - to take integral in polar coordinate (r,Fi)
Nd=1000;                    %number of pixels on the detector

k=2*pi/lambda;     
Delta=delta*2;  Beta=beta*2;  Qc=sqrt(Delta);    %eps=1-Delta+i*Beta
Din=Dex-Alpha*L;                                 %entrance aperture
Lex=Din/tan(Alpha)+L;                            %coordinate of exit aperture of the WG

ri=abs(Din/tan(Alpha));              
if Alpha>0 DD=Ri-ri; else DD=Ri+ri; end

Dmin=min(Din,Dex);
m_max=ceil(2*Dmin*Qc/lambda);            %maximum number of supported modes          
Qm(1:m_max)=(1:m_max)*lambda/2/Din;      %resonance angles

Q0_min=0; Q0_max=Qc;                     %limits grazing incidence angles(rad) 
dQ0=(Q0_max-Q0_min)/N0;

Fimin=0; Fimax=abs(Alpha);               %limits for Fi - to take integral
Fi=(0:Nfi-1)*Fimax/Nfi;

Qmax=4*Qc+Q0_max; dQd=Qmax/Nd;           %limits of angles in far field zone (rad)
NN=floor(N0*dQ0/dQd-abs(Alpha)/2/dQd);   % addition shift due to Q-2Q geometry

Qd=(-Nd/2+NN:Nd/2+NN-1)/Nd*Qmax;         %angle on detector in Q-2Q
yy=(0:Nd-1)*lambda/Qmax;
%-------------------------------------------------------------------------                    
    U_m=zeros(Nd,m_max); U=zeros(Nd); FF=zeros(Nd,N0);
    NNo=0;
    Q0=Q0_min-dQ0;

    for n=1:N0                               %number of measurements on angles
    process=floor(n/N0*100)
    Q0=Q0+dQ0; Qo(n)=Q0;  
                flag=0;if Q0<0|Q0>Qc*1.5;flag=1;else NNo=NNo+1;end
    for mm=1:m_max                       %number of modes   
                if flag==1break;end
       m=mm;
                  han1=besselh(m*pi/abs(Alpha),2,k*ri);
                  han2=besselh(k*Q0*Ri,1,k*Ri);
                  han2=han2/abs(han2);
                  norm=abs(Alpha)/pi/k/ri;
       A=han1*han2/norm*quadl(@Funh,Fimin,Fimax,1.e-11);
       Q_m=Qm(mm);              
       Ref=Refh(Q_m);
            for ii=1:Nd               %field on the exit aperture of the WG for the n-th measurment
              FI=abs(atan(yy(ii)/(Lex)));
                if FI>abs(Alpha)break;end      
              r=sqrt(Lex^2+yy(ii)^2);               
              han(ii)=besselh(mm*pi/abs(Alpha),1,k*r);
              U_m(ii,mm)=A*sin(m*pi*FI/abs(Alpha))*han(ii)*Ref;
            end% Nd
    end% mm
          U=sum(U_m,2);                %interf. patt. of all supported modes at the exit for the n-th mesurment
          Uf=fft(U);
          Ufs=fftshift(Uf);            %diffraction pattern in far field zone  
          F=Ufs;Uf=1;Ufs=1;  
    if flag==0
      FF(1:Nd,n)=sum(abs(F).^2,2);
      n1=floor(n*dQ0/dQd);
      FFF(1+n1:Nd+n1,n)=FF(1:Nd,n);    %diffraction pattern in far field zone in Q-2Q geometry 
    U_m=zeros(Nd,m_max); U=zeros(Nd); FF=zeros(Nd,N0);
    
    end
      flag=0;
end% N0
%--------------------------------------------------------------------------
NN=floor(N0*dQ0/dQd);
F1(1:Nd,:)=FFF(NN+1:Nd+NN,:);
figure
pcolor(Qo,Qd,F1)
shading interp
colorbar
%toc