%--------------------------------------------------------------------------
% function [I3] = FT_k_plpl_2nd_order(p,l,p_,l_,wn,X,n,m,A)
%
% A function for Matlab which analytically calculates the third term 
% (second order) 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 second order expansion of the exponential,
% -2.k^2.(Z_nm)^2.  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]
%
% I3:       Returned result for the third term (2nd order term).
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    03.07.2012
%--------------------------------------------------------------------------
%

function [I3] = FT_k_plpl_2nd_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
    % Different values for different combinations of l, l' and m
    if abs(2*m)==abs(l-l_)
        if m==0
            I3_phi = 2*pi;
        else
            I3_phi = (m/mm) * pi/2;
        end
    elseif l==l_
        I3_phi = pi;
    else
        I3_phi = 0;
    end
    
    if I3_phi==0
        I3(1:length(X)) = 0;
    else
        % Radial integration
        X3_sum = 0;
        for i=0:p
            for j=0:p_
                for h=0:(n-mm)/2
                    for g=0:(n-mm)/2
                        
                        % Factorial factors
                        C1 = 1/(factorial(p-i)*factorial(ll+i)*factorial(i));
                        C2 = 1/(factorial(p_-j)*factorial(ll_+j)*factorial(j));
                        C3 = factorial(n-h)/(factorial(h)*factorial((n+mm)/2-h)*factorial((n-mm)/2-h));
                        C4 = factorial(n-g)/(factorial(g)*factorial((n+mm)/2-g)*factorial((n-mm)/2-g));
                        
                        % Unormalise incmplete gamma function
                        N_factor=factorial(i+j+n-h-g+(ll+ll_)/2);
                        
                        % X terms
                        X3_term = X.^(h+g-n).*gammainc(X,i+j+n-h-g+(ll+ll_)/2+1);
                        X3_term = (-1)^(i+j+h+g)*C1*C2*C3*C4*N_factor*X3_term;
                        
                        % Sum
                        X3_sum = X3_sum + X3_term;
                    end
                end
            end
        end
        
        % Final third term integration result
        k3 = -(I3_phi/pi)*wn^2*A^2*sqrt(factorial(p)*factorial(p_)*factorial(p+ll)*factorial(p_+ll_));
        I3 = k3 * X3_sum;
    end


end

