%--------------------------------------------------------------------------
% function [k1] = FT_aperture_coupling_of_LG_modes(p,l,p_,l_,w,R)
%
% A function for Matlab which calculates the coupling coefficients for an
% LG mode incident on a circular aperture.  This function returns an
% analytical solution of the integral:
%
%      k = int_0^R ( int_0^2*pi ( U_p,l . U*_p'l' dphi) r dr)
% 
% where U_p,l is the incidnet mode and U_p'l' is the coupled mode.
%
% p,l:      Incident mode indices
% p_,l_:    Coupled mode indices
% w:        Beam spot size [m]
% R:        Aperture radius [m]
%
% k1:       Returned coupling coefficient
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    11.07.2012
%--------------------------------------------------------------------------
%

function [k1] = FT_aperture_coupling_of_LG_modes(p,l,p_,l_,w,R)

    % Change of variable
    X = 2*R.^2./(w.^2);

    % Only get coupling into modes with the same l as the incident mode
    % (angular integration)
    if l==l_
        
        dll_=1;
    else
        dll_=0;
    end
    
    X1_sum = 0;
    
    % Radial integration
    for i=0:p
        for j=0:p_
            
            % Calculate factors
            A1 = 1/(factorial(p-i)*factorial(abs(l)+i)*factorial(i));
            A2 = 1/(factorial(p_-j)*factorial(abs(l)+j)*factorial(j));
            
            N_factor = factorial(abs(l)+i+j);

            % Calculate X terms
            X1_term = (-1)^(i+j)*A1*A2*N_factor*gammainc(X,abs(l)+i+j+1);
            
            % Sum terms
            X1_sum = X1_sum+X1_term;
        end
    end
    
    % Final coupling coefficient
    k1 = dll_ * sqrt(factorial(p)*factorial(p_)*factorial(p+abs(l))*factorial(p_+abs(l)));
    k1 = k1*X1_sum;
    
           

end

