%
%--------------------------------------------------------------------------
% function [knn] = FT_zernike_coupling_of_LG_modes(p,l,p_,l_,lambda,w,R,A,n,m,no_terms)
%
% A function for Matlab which calculates an approximation of the 
% coupling coefficients from an input LG mode incident on a surface
% defined by a Zernike polynomial, into a new LG mode.  Use carefully for
% large surface distortions (i.e. A/lambda > 0.01) as this approximation
% depends on small A/lambda.  This function is based on the approximation:
%
%  U_p,l . exp(2ikZ_n,m) . U*_p_,l_  -->  U_p,l . (1 + 2ikZ_n,m - 2k^2(Z_nm)^2) . U*_p_,l_
%                                        Terms:    1st    2nd        3rd
%
% p,l:      input LG indices
% p_,l_:    LG mode coupled into
% lambda:   wavelength of incident beam [m]
% w:        beam radius at the surface [m]
% R:        radius of zernike polynomial [m]
% A:        amplitude of the zernike polynomial [nm]
% n,m:      indices of the Zernike surface causing the coupling.  -m
%           indicates odd polynomial (sin(m*phi)), +m indicates even
%           polynomial (cos(m*phi)).  The sign of m only effects the
%           phase of knn.
% no_terms: number of terms to be included in the approximation (2 or 3).
%           2 terms expands the exponential to 1st order, 3 terms is the 
%           maximum which expands the exponential to 2nd order.  
%
% knn:      returned coupling coefficent
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    18.05.2012            
%--------------------------------------------------------------------------
%

function [knn] = FT_zernike_coupling_of_LG_modes(p,l,p_,l_,lambda,w,R,A,n,m,no_terms)

    baseid = 'zernike_coupling_of_LG_modes';

    % Error messages
    % Mode/ zernike indices must be integers
    if (round(p)~=p || round(l)~=l || round(l_)~=l_ || round(p_)~=p_ || round(n)~=n || round(m)~=m)
        result='Invalid zernike/mode indices: p,l,p_,l_,n and m must be integers';
        msgid=[baseid,':checkarguments']
        error(msgid,result)
    end
    
    % p, p_ and n must be positive
    if (p<0 || p_<0 || n<0)
        result='Invalid zernike/mode indices: p,p_ and n must be positive';
        msgid=[baseid,':checkarguments']
        error(msgid,result)
    end
    
    % |m| <= n and n-m is even
    if (abs(m)>n || round((n-m)/2)~=(n-m)/2)
        result='Invalid zernike indices: m must be between -n and n, in steps of 2';
        msgid=[baseid,':checkarguments']
        error(msgid,result)
    end
    
    % no_terms must be a positive integer
    if (no_terms<1 || round(no_terms)~=no_terms)
        result='Invalid value: no_terms must be a positive integer';
        msgid=[baseid,':checkarguments']
        error(msgid,result)
    end

    % Amplitude into m
    A = 1e-9*A;
    
    % Wavenumber
    wn = 2*pi/lambda;
    
    % Change of variable
    R_w=R./w;
    X=2*R_w.^2;
    
    % First integration term
    if p==p_ && l==l_
        I1(1:length(X)) = 1;
    else
        I1(1:length(X)) = 0;
    end
    
    I2(1:length(X)) = 0;
    I3(1:length(X)) = 0;
    
    if no_terms>1
        % Second integration term
        %------------------------------------------------------------------
        I2 = FT_k_plpl_1st_order(p,l,p_,l_,wn,X,n,m,A);
        if no_terms>2
            if no_terms>3
                fprintf('Warning: approximation only includes 3 integration terms! \n')
            end
            % Third integration term
            %--------------------------------------------------------------
            I3 = FT_k_plpl_2nd_order(p,l,p_,l_,wn,X,n,m,A);
        end
    end
    
    % Final approximation of the coupling
    knn = I1+I2+I3;

  
end

