function value = EvalPolys3D(x,y,z,n_deg)

% Evaluate the orthonormal basis polynomials of degree less than or equal
% to n_deg. This is done using the triple recursion relation for these
% orthonormal polynomials.
%
% The constants needed when using the triple recursion relation are
% calculated here by the following functions, all of which are given here: 
%               Ax, Ay1, Ay2, Az1, Az2, Az3, Cx, Cz3.
% When the orthonormal polynomials are to be evaluated repeatedly, the
% constants from these function can be evaluated once and then stored for
% repeated use. This provision is not included here.
%
% INPUT:
% (x,y,z) - The points at which the orthonormal polynomials are to be
%           evaluated.
% n_deg -   This is the upper limit on the degrees of the orthonormal
%           polynomials to be calculated.
%
% OUTPUT:
% value -   A three dimensional array giving the values of the orthonormal
%           polynomials. The partition value(:,:,:,j) contains the values
%           of the orthonormal polynomial #j, using the lexicographical
%           ordering of the indices of the orthonormal polynomials 
%           calculated by the program NthFunction. The first three 
%           dimensions are those of the dimensions of x, y, and z. 

N = (n_deg+3)*(n_deg+2)*(n_deg+1)/6;
[nr, nc, nh]=size(x);
value = zeros(nr,nc,nh,N);
const = sqrt(3/(4*pi));
value(:,:,:,1)= const;

if n_deg == 0
    return
end

const2 = sqrt(5)*const;
value(:,:,:,2)= const2*x; 
value(:,:,:,3)= const2*y;  
value(:,:,:,4)= const2*z;

if n_deg == 1
    return
end

for n = 1: n_deg-1
    for sumJK = 0:n-1
        for k=0:sumJK
            j = sumJK - k;
            n2 = NthFunction(n+1,j,k);
            n1 = NthFunction(n,j,k);
            n0 = NthFunction(n-1,j,k);
            value(:,:,:,n2) = (x.*value(:,:,:,n1) - ...
                Cx(n,j,k)*value(:,:,:,n0))*(1/Ax(n,j,k));
        end
    end
    
    for k=0:n
        j= n - k;
        n2 = NthFunction(n+1,j,k);
        n1 = NthFunction(n,j,k);
        value(:,:,:,n2) = (x.*value(:,:,:,n1))*(1/Ax(n,j,k));
    end
    
    for k=0:n-1
        j = n - k;
        n3 = NthFunction(n+1,j+1,k);
        n2 = NthFunction(n+1,j-1,k);
        n1 = NthFunction(n,j,k);
        n0 = NthFunction(n-1,j-1,k);
        value(:,:,:,n3) = (y.*value(:,:,:,n1) - ...
            Ay2(n,j,k)*value(:,:,:,n2) - Cy1(n,j,k)*value(:,:,:,n0))...
            *(1/Ay1(n,j,k));
    end
    
    n2 = NthFunction(n+1,1,n);
    n1 = NthFunction(n,0,n);
    value(:,:,:,n2) = (y.*value(:,:,:,n1))*(1/Ay1(n,0,n));
    
    n4 = NthFunction(n+1,0,n+1);
    n3 = NthFunction(n+1,0,n-1);
    n2 = NthFunction(n+1,2,n-1);
    n1 = NthFunction(n,0,n);
    n0 = NthFunction(n-1,0,n-1);
    value(:,:,:,n4) = (z.*value(:,:,:,n1) - Az3(n,0,n)*value(:,:,:,n3) -...
        Az2(n,0,n)*value(:,:,:,n2) - Cz3(n,0,n)*value(:,:,:,n0))*...
        (1/Az1(n,0,n));
end

end % End Evalpolys3D

% =====================================================================

function value = Ax(n,j,k)

% Produce the coefficients a_(n,j,k) used in setting up the triple
% recursion relation for the classical orthonormal polynomials.

% Initialize
% Mu can be any number between 0 and 1, but this program uses 1/2.
Mu = 1/2;
Mu2 = 1;
value = 0.5*sqrt(((Mu2+j+k+n+2)*(n+1-j-k))/((Mu+n+2)*(Mu+n+1)));

end % Ax

% =====================================================================

function value = Ay1(n,j,k)

% This function is for a^[n,y]_(j,k;j+1,k)

Mu = 1/2;
Mu2 = 1;

numerator =(Mu2+j+2*k+1)*(j+1)*(Mu2+j+k+n+3)*(Mu2+j+k+n+2);
denominator =(Mu+j+k+ 0.5)*(Mu+j+k+1.5)*(Mu+n+2)*(Mu+n+1);

value = 0.25* sqrt(numerator/denominator);

end % Ay1

% =====================================================================

function value = Ay2(n,j,k)

% This function is for a^[n,y]_(j,k;j-1,k)

Mu = 1/2;
Mu2 = 1;

numerator = j*(Mu2+j+2*k)*(n+2-j-k)*(n+1-j-k);
denominator =(Mu+j+k+ 0.5)*(Mu+j+k-0.5)*(Mu+n+1)*(Mu+n+2);

value = -0.25* sqrt(numerator/denominator);

end % Ay2

% =====================================================================

function value = Az1(n,j,k)

% This function is for a^[n,z]_(j,k;j,k+1)

Mu = 1/2;
Mu2 = 1;

numerator = (k+1)*(Mu2+k)*(Mu2+j+2*k+2)*(Mu2+j+2*k+1)*(Mu2+j+k+n+3)* ...
    (Mu2+j+k+n+2);
denominator =(Mu+k)*(Mu+k+1)*(Mu+j+k+1.5)*(Mu+j+k+0.5)*(Mu+n+2)* ...
    (Mu+n+1);

value = 0.125* sqrt(numerator/denominator);

end % Az1

% =====================================================================

function value = Az2(n,j,k)

% This function is for a^[n,z]_(j,k;j+2,k-1)

Mu = 1/2;
Mu2 = 1;

numerator = k*(Mu2+k-1)*(j+2)*(j+1)*(Mu2+j+k+n+3)*(Mu2+j+k+n+2);
denominator =(Mu+k)*(Mu+k-1)*(Mu+j+k+0.5)*(Mu+j+k+1.5)*(Mu+n+2)* ...
    (Mu+n+1);

value = -0.125* sqrt(numerator/denominator);

end % Az2

% =====================================================================

function value = Az3(n,j,k)

% This function is for a^[n,z]_(j,k;j,k-1)

Mu = 1/2;
Mu2 = 1;

numerator = k*(Mu2+k-1)*(Mu2+j+2*k)*(Mu2+j+2*k-1)*(n+2-j-k)*(n+1-j-k);
denominator =(Mu+k)*(Mu+k-1)*(Mu+j+k+0.5)*(Mu+j+k-0.5)*(Mu+n+2)* ...
    (Mu+n+1);

value = -0.125* sqrt(numerator/denominator);

end % Az3

% =====================================================================

function value = Cx(n,j,k)

value = Ax(n-1,j,k);

end % Cx

function value = Cy1(n,j,k)

% This function correponds to a^[n,y]_(j,k;j+1,k)

value = Ay1(n-1,j-1,k);

end % Cy1

% =====================================================================

function value = Cz3(n,j,k)

% This function is for c^[n,z]_(j,k;j,k-1), which equals Az1(n-1,j, k-1)
% For our purpose, j=0 and k=n.

% Mu = 1/2;
% Mu2 = 1;
value = Az1(n-1,j,k-1);

end % Cz3