function showmodel(model,dim,p,q,maxrad)

% showmodel(model,dim,p,q,max)
% 2-D or 3-D plotting of clusters of points
% The diameter of spheres is the p-percentile distance between the nearest neighbor
% within cluster members.
%
% Input:
% model: generated model
% dim: is [2] or [3]
% p: percentile (0<p<1) (p=0.5 stands for the median)
% q: quality of the 3-D image [use 20 for default]
% max: maximal radius
%
% 0:blue 1:red 2:magenta 3:cyan 4:green 5:yellow 6:dark purple 7:orange 8:gray 9:dark green
% 10: 11: 12:
%
% Kevin Dawson 11/02/3

if nargin<1
    help showmodel;
    return
end


if nargin <4
    q=50;
end
if nargin <3
    p=.95;
end
if nargin <2
    dim=3;
elseif dim<2 | dim>3 | dim ~=  round(dim)
    error ('dim should be 2 or 3');
end


catN= size(model.groupstart,2)-1;

if catN > 60
    error('Too many categories');
end

if nargin <5
    maxrad=0;
else
    maxd=maxrad*2;
end


 cols={[0,0,1] [1,0,0] [1,0,1] [0,1,1] [0,1,0] [1,1,0] [.4,.1,.2] [1,.5,0] [.7,.7,.7] [.2,.5,.2] [.2 .2 .2] [.3 .4 .5] [0 .3 .4] [.1 .1 .1] [.1 .1 .3] [.1 .1 .5] [.1 .1 .7] [.1 .1 .9] [.1 .3 .1] [.1 .3 .3] [.1 .3 .5] [.1 .3 .7] [.1 .3 .9] [.1 .5 .1] [.1 .5 .3] [.1 .5 .5] [.1 .5 .7] [.1 .5 .9] [.1 .7 .1] [.1 .7 .3] [.1 .7 .5] [.1 .7 .7] [1. .7 .9] [.1 .9 .1] [.1 .9 .3] [.1 .9 .5] [.1 .9 .7] [.1 .9 .9] [.3 .1 .1] [.3 .1 .3] [.3 .1 .5] [.3 .1 .7] [.3 .1 .9] [.3 .3 .1] [.3 .3 .3] [.3 .3 .5] [.3 .3 .7] [.3 .3 .9] [.3 .5 .1] [.3 .5 .3] [.3 .5 .5] [.3 .5 .7] [.3 .5 .7] [.3 .5 .9] [.3 .7 .1] [.3 .7 .3] [.3 .7 .5] [.3 .7 .7] [.3 .7 .9] [.3 .9 .1] [.3 .9 .3]};
% 1:blue 2:red 3:magenta 4:cyan 5:green 6:yellow 7:dark purple 8:orange 9:gray 10:dark green
%cols={[0,0,1] [0,1,1] [1 1 0] [1,0,0]}; %blue-cyan-yellow-red
%cols={[0,0,1] [0,1,1] [1,0,0]}; %blue-cyan-red
%cols={[0 0 1] [ 0 1 1] [1 0 0] [1 0 0]} %2s
%cols={[0 0 0] [.5 .5 .5] [1 1 1] [1 1 1]} %2sb
%cols={[0,0,.5] [1,0,0] [0,1,1]} %2l
%cols={[0 0 0] [1 1 1] [.05 .05 .05]} %2lb

% cols={[1 1 1] [.7 .7 .7] [.6 .6 .6] [.5 .5 .5] [.4 .4 .4] [.3 .3 .3] [.2 .2 .2] [.1 .1 .1] [0 0 0]};
% cols={[1 1 1] [.8 .8 .8][.6 .6 .6] [.4 .4 .4] [.2 .2 .2] [0 0 0]};
% cols={[1 1 1] [.7 .7 .7] [.1 .1 .1]};
% cols={[1 1 1] [.6 .6 .6] [.4 .4 .4] [0 0 0]};
%for i=1:60;c=(i-1)/59;cols{i}=[c c c];end

figure;
hold on;
grid on;
for i=1:catN
    from=model.groupstart(i);
    to=model.groupstart(i+1)-1;
    subM=model.coords{dim}(:,from:to);
    subD=L2_distance(subM,subM,1);
    n=size(subD,1);
    for j=1:n
        subD(j,j)=NaN;
    end
    d=sort(min(subD));
    j=size(d,2);
    k=round(j*p);
    if k<1
        k=1;
    end
    if k>j
        k=j;
    end
    d=d(k);
    r=d/2;
    if maxrad~=0,
        r=min(r,maxd);
        d=2*r;
    end
    
    
    
    if dim==2
        for j=1:n
            rectangle('Position',[subM(1,j)-r,subM(2,j)-r,d,d],'Curvature',[1,1],'EdgeColor','k','FaceColor',cols{model.phenos(i)});
        end
    else
        light('Color','w');
        for j=1:n
            [Xs,Ys,Zs]=ellipsoid(subM(1,j),subM(2,j),subM(3,j),r,r,r,q);
            surf(Xs,Ys,Zs,'FaceLight','gouraud','EdgeColor','none','FaceColor',cols{model.phenos(i)},'DiffuseStrength',.2,'AmbientStrength',.2);
        end
    
    end
end
hold off;


return