% polygon(A,b,M)
%
% Draws the polygon determined by the system Ax<=b
% For example, Modellen & Simulatie page 82:
%  Pagina 112 (2008 editie)
% A = [1,2,2;-2,1,0;1,-2,-1;0,-1,-1;-2,-3,-3]
% b = [3;0;0;0;-1]
% M = rand(1,3)
%
% Use matlab figure controls to rotate, zoom, etc.
%
% -Arno Swart [swart@math.uu.nl]


function polygon(A,b,M)

nreqs = size(A,1);
nrvar = size(A,2);
disp(sprintf('%d equations in %d unknowns.',nreqs,nrvar));
disp(sprintf('Determining potential corner points...'));

candidates = nchoosek([1:nreqs],nrvar);
cornerpoints = [];
for i=1:size(candidates,1)
    if rank( A(candidates(i,:),:))==nrvar   
        cornerpoints=[cornerpoints; (A(candidates(i,:),:)\b(candidates(i,:)))'];
    end
end

cornerpoints;

if size(cornerpoints,1)==0
    disp(sprintf('No Cornerpoints!'));
    return;
end

disp(sprintf('Determining corner points inside polytope...'));

evaluated = A*cornerpoints';
cps=[];
Mval=[];
tol=1e-5;
for i=1:size(evaluated,2)
    if sum(evaluated(:,i) <= b+tol)==nreqs;
        cps=[cps;cornerpoints(i,:)];
        Mval = [Mval,M*(cornerpoints(i,:)')];
    end
end

K= convhulln(cps);
h=trisurf(K,cps(:,1),cps(:,2),cps(:,3),Mval);
set(h,'linewidth',2);
set(h,'edgecolor',[0,0,0]);
set(h,'facecolor','interp');
colorbar;
