%--------------------------------------------------------------------------
% function [Rnm] = FT_Rnm(n,m,rho)
%
% A function for Matlab which calculates the radial polynomial, Rnm, part 
% of the Zernike polynomials, for a given n and m.  For n>25 the radial
% polynomials are calculated via a recurrence formula to prevent errors
% resulting from large numbers generated in the original calculation of
% Rnm.
%
% n:    Radial coordinate
% m:    Azimuthal coordinate
% rho:  Normalised radial coordinate (Zernike polynomials are only valid in
%       region with rho<=1).
%
% Rnm:  Returned radial polynomial values
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    02.05.2012
%--------------------------------------------------------------------------
%

function [Rnm] = FT_Rnm(n,m,rho)
    
    baseid = 'Rnm';

    % Error messages
    % Indices must be integers
    if (round(n)~=n || round(m)~=m)
        result='Invalid zernike indices: n and m must be integers';
        msgid=[baseid,':checkarguments']
        error(msgid,result);
    end
    % n must be positive or 0
    if (n<0)
        result='Invalid zernike indices: n must be 0 or positive';
        msgid=[baseid,':checkarguments']
        error(msgid,result);
    end
    % Integer difference between n and m must be even (otherwise polynomial is 
    % 0, but we don't want to include these).
    if (abs(m)>n || round((m-n)/2)~=(m-n)/2)
        result='Invalid zernike indices: m must be between -n and +n, in steps of 2';
        msgid=[baseid,':checkarguments']
        error(msgid,result);
    end

    Rnm=0;

    % Only use magnitude of m
    m = abs(m);
    
    % If n<=25 we calculate the radial polynomial with original formula,
    % otherwise use recurrence formula as factorials lead to extreme large
    % numbers.
    n0 = 25;
    
    if n<=n0
        % Calculating radial term
        for k=0:((n-m)/2)
    
            a=((-1)^k)*factorial(n-k)/(factorial(k)*factorial((n+m)*0.5-k)*factorial((n-m)*0.5-k));
            rho_power=rho.^(n-(2*k));
    
            rad_term=a*rho_power;
            Rnm=Rnm+rad_term;
        end
    else
        % calculating using recurrence relation
        Rnm = FT_Rnm_recurrence(n,m,rho);
    end
    
    


end

