function [ r, seed ] = user ( dim_num, n, seed, usertype, varargin )

%% USER samples points in a user-specified region with given density.
%
%  Discussion:
%
%    This routine can be used to 
%
%    * specify an interesting initial configuration for the data,
%      by specifing that USER be used for initialization (INIT = 3);
%
%    * specify the shape of the computational region, by specifying
%      that sample points are to be generated by this routine, 
%      (SAMPLE = 3) and then returning sample points uniformly at random.
%
%    * specify the distribution or density function, by specifying
%      that sample points are to be generated by this routine, 
%      (SAMPLE = 3 ) and then returning sample points according to a 
%      given probability density function.
%
%  Author:
%
%    John Burkardt
%
%  Modified:
%
%    23 June 2005
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer N, the number of sample points desired.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Input, string USERTYPE, specifying what sort of user data are supplied.
%
%    Input, struct USERDATA, user data.
%
%    Output, real R(DIM_NUM,N), a set of N sample values.
%
%    Output, integer SEED, the updated seed.  
%
    
    switch usertype
        case 'ellipse'
            r = randInEllipse( n, 0.3, 0.5 );
            if dim_num > 2
                r = [ r zeros( n, dim_num-2 ) ];
            end
        case { '', 'unitcircle' }
            r = randInCircle( n, 1 );
            if dim_num > 2
                r = [ r zeros( n, dim_num-2 ) ];
            end
        case 'meshcells'
            % Following arguments are the mesh and an array of cells.
            mesh = varargin{1};
            cellIndexes = varargin{2};
            r = randInFEMcells( mesh, cellIndexes, n )';
            r = r(1:dim_num,:);
        %    cellAreas = mesh.cellareas( cellIndexes );
        %    [cells,bcs] = randInTriangles( cellAreas, n );
        %    r = meshBaryToGlobalCoords( mesh, cellIndexes(cells), bcs )';
        case 'polygon'
            % Following argument is a convex 2D polygon.
            if isempty(varargin)
                userdata = [ [0 0];[2 -1];[2.5 1];[4 0.5];[3 2]; [1.5 4];[1 0]];
              % userdata = [ [0 0];[1 0];[1 1];[0 1]];
                r = randInPoly( n, userdata )';
            else
                r = randInPoly( n, varargin{1} )';
            end
            if dim_num > 2
                r = [ r zeros( n, dim_num-2 ) ];
            end
        case 'triangles'
            % Following arguments represent a list of triangles.
            % varargin{1} is a list of node coordinates.
            % varargin{2} is a list of triples of node indexes.
            % varargin{3} is a list of areas of the triangles.
            % varargin{4} is a list of indexes into varargin{2}.
            r = randInTriangles( n, varargin{:} );
            if dim_num > 2
                r = [ r zeros( n, dim_num-2 ) ];
            end
    end
end

