function [hypval,cfs]=hyperlissa(deg,pts,f3);

% Stefano De Marchi and Marco Vianello, February 2015

% computes the hyperinterpolation polynomial of a trivariate function
% at a rank-1 Chebyshev lattice on a Lissajous curve of the cube

% a call has the form [hypval,cfs]=hyperlissa(deg,pts,@f3) where f3.m
% is a user defined trivariate Matlab function

% NOTE: to be used in connection with the Chebfun package
% http://www.chebfun.org

% input:
% deg: hyperinterpolation degree
% pts: 3-column array of evaluation points

% output:
% hypval: 1-column array of values of the hyperinterpolation polynomial
% at the points pts
% cfs: hyperinterpolation coefficients

n=deg;
if mod(n,2)==0
a=3/4*n^2+1/2*n;
b=3/4*n^2+n;
c=3/4*n^2+3/2*n+1;
else
a=3/4*n^2+1/4;
b=3/4*n^2+3/2*n-1/4;
c=3/4*n^2+3/2*n+3/4;
end
mu=n*c+1;


% chebfun along the Lissajous curve
g=chebfun(@(t) f3(cos(a*acos(t)),cos(b*acos(t)),cos(c*acos(t))),mu+1);
cc=chebpoly(g);
% cc=chebcoeffs(g) with Chebfun Version 5 
u=pi*cc;
u(2:end-1)=u(2:end-1)/2;
u=fliplr(u);


% indexes i,j,k: graded lexicographical order
ind=(0:1:n);
[j1,j2,j3]=meshgrid(ind);

for s=0:n
good=find(j1(:)+j2(:)+j3(:)==s);
i(1+s*(s+1)*(s+2)/6:(s+1)*(s+2)*(s+3)/6)=j1(good);
j(1+s*(s+1)*(s+2)/6:(s+1)*(s+2)*(s+3)/6)=j2(good);
k(1+s*(s+1)*(s+2)/6:(s+1)*(s+2)*(s+3)/6)=j3(good);
end;


% hyperinterpolation coeffs
alpha1=i*a+j*b+k*c+1;
alpha2=abs(i*a+j*b-k*c)+1;
alpha3=abs(i*a-j*b)+k*c+1;
alpha4=abs(abs(i*a-j*b)-k*c)+1;

cfs=pi^2/4*sigma(i*a).*sigma(j*b).*sigma(k*c).*(u(alpha1)+u(alpha2)+ ...
u(alpha3)+u(alpha4));

% hyperinterpolation values
x=pts(:,1);y=pts(:,2);z=pts(:,3);

for s=1:length(x)
tix=cos(i.*acos(x(s))).*sigma(i);
tjy=cos(j.*acos(y(s))).*sigma(j);
tkz=cos(k.*acos(z(s))).*sigma(k);
val(s)=cfs*(tix.*tjy.*tkz)';
end
hypval=val';

end


function sigmaval=sigma(m)

sigmaval=sqrt(1+sign(m))./sqrt(pi);

end


