%--------------------------------------------------------------------------
% function [I2] = FT_k_plpl_1st_order(p,l,p_,l_,wn,X,n,m,A)
%
% A function for Matlab which analytically calculates the second term (the 
% 1st order term) in an expansion of the integral of:
% 
% U_pl . exp(2.i.k.Z_nm) . U*_p'l' -> 1 + 2.i.k.Z_nm - 2.k^2.(Z_nm)^2 . . .
%
% i.e. the term containg the first order expansion of the exponential,
% 2.i.k.Znm.  This integral represents the coupling from an LG mode U_pl to
% a mode U_p'l' via a surface defined by a Zernike polynomial Z_nm.
%
% p,l:      Indices of the incident mode
% p_,l_:    Indices of the reflected mode
% wn:       Wavenumber [m^-1]
% X:        X = 2*R^2/(w^2) where R is the radius of the Zernike polynomial
%           and w is the beam radius.
% n,m:      Indices of the Zernike polynomial
% A:        Amplitude of the zernike polynomial [m]
%
% I2:       Returned result for the second term (1st order term).
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    03.07.2012
%--------------------------------------------------------------------------
%

function [I2] = FT_k_plpl_1st_order(p,l,p_,l_,wn,X,n,m,A)

    % Only need |m|,|l|,|l'| for radial integration
    mm = abs(m);
    ll=abs(l);
    ll_=abs(l_);
    
    % Angular integration
    if abs(l-l_)==abs(m)
        if m>0
            I2_phi = pi;
        elseif m==0
            I2_phi = 2*pi;
        else
            I2_phi = (l-l_)/abs(m) * 1i * pi;
        end
    % The case where |l-l'| is not equal to m the coupling from this term
    % is 0
    else
        I2_phi = 0;
    end
    
    if I2_phi==0
        I2(1:length(X)) = 0;
    else
        % Radial integration
        X2_sum = 0;
        for i=0:p
            for j=0:p_
                for h=0:(n-mm)/2
                    
                    % Calculate factors
                    B1=1/(factorial(p-i)*factorial(ll+i)*factorial(i));
                    B2=1/(factorial(p_-j)*factorial(ll_+j)*factorial(j));
                    B3=factorial(n-h)/(factorial(h)*factorial((n+mm)/2-h)*factorial((n-mm)/2-h));
                    
                    % Need to un-normalise the Matlab incomplete gamma
                    % function
                    N_factor=factorial(i+j-h+(ll+ll_+n)/2);
                    
                    % Calculate X terms
                    X2_term=(1./(X.^((n-2*h)/2))).*gammainc(X,i+j-h+1+(ll+ll_+n)/2);
                    X2_term=(-1)^(i+j+h)*B1*B2*B3*N_factor*X2_term;
                    
                    % Add up terms
                    X2_sum=X2_sum+X2_term;
                end
            end
        end
        
        % Final coupling coefficient
        k2 = (1i * I2_phi/pi) * A * wn * sqrt(factorial(p)*factorial(p_)*factorial(p+ll)*factorial(p_+ll_));
        I2 = k2 * X2_sum;
    end

end

