clear all
N=4;
M=N/2;
piN = pi/N;
Phi   = (-N :N )'*piN;
Theta = (-M :M )'*piN;
Ni = ceil(N*cos(Theta)); Ni(1)=1; Ni(2*M+1)=1;
piNi = pi./Ni; 
% define full and reduced grid
for i = -M:M
   Phii{i+M+1}  = (-Ni(M+1+i):Ni(M+i+1)-1)'*piNi(i+M+1);
end
% enter function on the reduced grid, then interpolate to the fixed
% grid Phi = -N:N
for i = 1:2*M+1
   F{i}=cos(Theta(i))^2.*sin(Phii{i});
end
for i = 1:2*M+1
   FT{i}=CSInterp(F{i});
end
for i=1:2*M+1
   Fvals(:,i) = CSeval(FT{i},pi,Phi);
end

