function [a_tot,c_tot,d_tot] = TripleRecurCoeff(n_deg)

% Produce the coefficients a_{n,j},c_{n,j},d_{n,j} used in setting up the
% triple recursion relation for the classical orthonormal polynomials. 
% This routine produces these coefficients from the matrices A_{n,1}, 
% A_{n,2} of Dunkl & Xu (p. 88) for n=1,2,...,n_deg. These will allow 
% calculating the orthonormal polynomials of degree less than or equal to
% n_deg. This is done in the routine 'EvalOrthoPolys'.
%
% WARNING: Use only with n_deg > 0. 

if  n_deg < 1
    disp('In TripleRecurCoeff, the value of n_deg is too small.')
    stop
end

% Initialize
N = (n_deg*(n_deg+3))/2;
a_tot = zeros(N,1); c_tot = a_tot; d_tot = a_tot;
jump = 1;

% Calculate the entries in a_tot, b_tot, and c_tot.
for n=1:n_deg
    [a,c,d] = TripleRecurCoeff_unit(n);
    a_tot(jump:jump+n) = a;
    c_tot(jump:jump+n) = c;
    d_tot(jump:jump+n) = d;
    jump = jump+n+1;
end

end % TripleRecurCoeff_total

function [a,c,d] = TripleRecurCoeff_unit(n)

% Produce the coefficients {a_n} for the matrix A_1 and the coefficients
% {c_n} and {d_n} for the matrix A_2. See page 88 of Dunkl & Xu.

a = zeros(n+1,1); c = a; d = a;
% a(k) will contain the coefficient a_{k-1,n}, k=1,...,n+1;
% c(k+1) will contain the coefficient c_{k,n}, k=1,...,n;
% d(k) will contain the coefficient d_{k-1,n}, k=1,...,n+1.
%
% WARNING: Use only with n > 0. 

if  n < 1
    disp('In TripleRecurCoeff_unit, the value of n is too small.')
    stop
end

for k=0:n
    a(k+1) = 0.5*sqrt(((n-k+1)/(n+1))*((n+k+2)/(n+2)));
    d(k+1) = 0.5*sqrt(((k+1)/(2*k+1))*((k+1)/(2*k+3))...
                                   *((n+k+3)/(n+1))*((n+k+2)/(n+2)));
end

for k=1:n
    c(k+1) = -0.5*sqrt((k/(2*k+1))*(k/(2*k-1))*((n-k+1)/(n+1))...
                                   *((n-k+2)/(n+2)));
end

end % TripleRecurCoeff_unit