
function pts=disk_wam(n)

%--------------------------------------------------------------------------
% OBJECT.
%----------
%
% COMPUTE WAM POINTS FOR UNIT DISK "B([0,0],1)".
%
%--------------------------------------------------------------------------
% INPUT.
%----------
%
% n:  WAM DEGREE.
%
%--------------------------------------------------------------------------
% OUTPUT.
%----------
%
% pts: WAM POINTS ON THE UNIT DISK "B([0,0],1)".
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%% Copyright (C) 2010 
%% Len Bos, Stefano De Marchi, Alvise Sommariva, Marco Vianello.
%%
%% This program is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation; either version 2 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; if not, write to the Free Software
%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
%%
%% Author:  
%%          Len Bos           <leonardpeter.bos@univr.it>
%%          Stefano De Marchi <demarchi@math.unipd.it>
%%          Alvise Sommariva  <alvise@math.unipd.it>
%%          Marco Vianello    <marcov@math.unipd.it>    
%%
%% Date: October 18, 2010.
%--------------------------------------------------------------------------

if (mod(n,2) == 0)
    % SOME ZEROS COULD BE REPEATED IN THE ORIGINAL FORMULATION THAT
    % HAS BEEN CONSEQUENTLY MODIFIED.
    j=0:n; j=setdiff(j,n/2); % WE ADD [0 0] LATTER.
    k=0:n;
    [rho,theta]=meshgrid(cos(j*pi/n),k*pi/(n+1));
    meshA=[rho(:).*cos(theta(:)) rho(:).*sin(theta(:))];
    % ADDING [0 0]        
    nodes_x=[meshA(:,1); 0];
    nodes_y=[meshA(:,2); 0];
else
    % ALL POINTS ARE DISTINCT.
    j=0:n; 
    k=0:n;
    [rho,theta]=meshgrid(cos(j*pi/n),k*pi/(n+1));
    meshA=[rho(:).*cos(theta(:)) rho(:).*sin(theta(:))];   
    nodes_x=meshA(:,1);
    nodes_y=meshA(:,2);
end

pts=[nodes_x nodes_y];

plot (pts(:,1),pts(:,2), '+');
