%--------------------------------------------------------------------------
% function [clipping_loss] = FT_LG_clipping_loss(p,l,w,R)
%
% A function for Matlab which calculates the clipping loss experienced by a
% particular Laguerre-Gauss beam incident an aperture for a given beam and 
% aperture radius.  This function calculates the clipping from an 
% analytical solution of the equation:
%
% p_loss = 1 - int_0^R U_pl x U_pl* 2*pi r dr
%
% where p_loss is the power clipped, R is the radius of the aperture and 
% U_pl is the incident LG mode 
%
%  p,l:             Indices of the LG beam
%  w:               Beam radius at the aperture [m]
%  R:               Aperture radius [m]
%
%  clipping_loss:	The resulting loss of power due to clipping at the
%		            aperture
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
%  Charlotte Bond	18.12.2010 
%--------------------------------------------------------------------------
%

function [clipping_loss] = FT_LG_clipping_loss(p,l,w,R)
 
    % Ratio of aperture to beam radius 
    R_w=R./w;

    % Set x limits
    X=2*R_w.^2;
    
    P_ref=0;    
    l=abs(l);
    
    % Calculate the power reflected through the aperture
    for m=0:p
        for n=0:p
            N=(-1)^(n+m)/(factorial(p-n)*factorial(p-m)*factorial(l+n)*factorial(l+m)*factorial(n)*factorial(m));
            P_term=factorial(p)*factorial(p+l)*N*factorial(l+n+m).*gammainc(X,l+n+m+1);
            P_ref=P_ref+P_term;
        end
    end
    
    % The clipping loss
    clipping_loss=1-P_ref;            
end

