%--------------------------------------------------------------------------
% [Rnm] = FT_Rnm_recurrence(n,m,rho,n0)
%
% A Matlab function which calculates the radial polynomial, Rnm, part of 
% the Zernike polynomials via a recurrence method, rather than from the 
% original formula (FT_Rnm.m).  Used for calculating Zernike polynomials/ 
% Rnm for large n, i.e. n>25, as the factorials involved in FT_Rnm become 
% very large, producing errors at lareg n.  Rnm is calculated using Chong 
% and Kintner relations.
%
% n:    Radial index
% m:    Azimuthal index
% rho:  Normalised radial position
%
% Rnm:  Returned values for radial polynomial.
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    02.05.2012
%--------------------------------------------------------------------------
%

function [Rnm] = FT_Rnm_recurrence(n,m,rho)

    % Use magnitude of m
    m = abs(m);
    
    % Calculate Rnm recurrence using Chong and Kintner relations
    % Modified Kintner method:
    Rnm_4 = rho.^n;
    Rnm_2 = n*rho.^n - (n-1)*rho.^(n-2);
    
    % If m=n or m=n-2 Rnm is given by results above (i.e. modified Kintner
    % method)
    if n==m
        Rnm = Rnm_4;
    elseif n-2==m
        Rnm = Rnm_2;
    else    
        % Otherwise calculate from recurrence (Chong) relation
        for m_=n-4:-2:m
        
            h3 = -4*(m_+2)*(m_+1)/((n+m_+2)*(n-m_));
            h2 = h3*(n+m_+4)*(n-m_-2)/(4*(m_+3)) + (m_+2);
            h1 = (m_+4)*(m_+3)/2 - (m_+4)*h2 + h3*(n+m_+6)*(n-m_-4)/8;
        
            Rnm = h1*Rnm_4 + (h2+(h3./(rho.^2))).*Rnm_2;
        
            Rnm_4 = Rnm_2;
            Rnm_2 = Rnm;        
        end
    end

end

