function value = EvalOrthoPolys(x,y,n_deg,a_tot,c_tot,d_tot)
%
% Evaluate the orthonormal basis polynomials of degree less than or equal
% to n_deg. The column vectors a_tot,c_tot,d_tot contain the coefficients 
% to be used in the evaluation, doing so using the triple recursion 
% relation for those polynomials. These constant vectors are produced using 
% 'TripleRecurCoeff'. See Dunkl % Xu (p. 88) for these polynomials.
%
% INPUT:
% (x,y) -   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.
% a_tot,c_tot,d_tot - These constants are to be used with the triple
%           recursion relation for the orthonormal polynomials.
%
% 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 standard
%           lexicographical ordering of the indices of the orthonormal
%           polynomials. The first two dimensions are those of the
%           dimensions of x and y. 

N = ((n_deg+1)*(n_deg+2))/2;
[nr,nc] = size(x);
value = zeros(nr,nc,N);

rsqrtpi = 1/sqrt(pi);
value(:,:,1) = rsqrtpi*ones(size(x));

if n_deg == 0 
    return
end

value(:,:,2) = 2*rsqrtpi*x; value(:,:,3) = 2*rsqrtpi*y;
if n_deg == 1
    return
end

value(:,:,4) = rsqrtpi*(4*x.^2 - 1);
value(:,:,5) = rsqrtpi*sqrt(24)*x.*y;
value(:,:,6) = rsqrtpi*sqrt(2)*(3*y.^2 + x.^2 -1);
if n_deg == 2
    return
end

p0 = value(:,:,2:3); p1 = value(:,:,4:6);
a0 = a_tot(1:2); d0 = d_tot(1:2); 
jump = 3;
step = 7;
for n=2:n_deg-1
    a1 = a_tot(jump:jump+n); 
    c1 = c_tot(jump:jump+n); 
    d1 = d_tot(jump:jump+n);
    p2 = zeros(nr,nc,n+2);
    for j=0:n-1
        p2(:,:,j+1) = (x.*p1(:,:,j+1)-a0(j+1)*p0(:,:,j+1))/a1(j+1);
    end
    p2(:,:,n+1) = x.*p1(:,:,n+1)/a1(n+1);
    p2(:,:,n+2) = (y.*p1(:,:,n+1)-c1(n+1)*p2(:,:,n)-d0(n)*p0(:,:,n))/d1(n+1);
    value(:,:,step:step+n+1) = p2;
    p0 = p1; p1 = p2;
    a0 = a1; d0 = d1;
    jump = jump + n + 1;
    step = step+n+2;
end

end % EvalOrthoPolys