function s = find2dPIS(u,ABCD,options)%function s = find2dPIS(u,ABCD,options)%Find a positively invariant set for the 2nd-order binary modulator whose %loop filter is described by ABCD and whose input is a constant u.%Return s = [] if no invariant set is found.%%options = [ dbg=0 itnLimit=100 expFactor=0.005 N=1000 skip=100]%dbg=1 causes debugging information to be displayed.s=[];n = size(ABCD,1)-1;if n ~= 2    fprintf('%s: Error. The modulator must be second order.\n', mfilename);    return;end% Handle the optionsnlev=2;parameters = ['dbg      ';'itnLimit ';'expFactor';'N        ';'skip     '];defaults = [ 0 100 .005 1000 100];for i=1:length(defaults)    if i>length(options)    	eval([parameters(i,:) '=defaults(i);'])      elseif isnan(options(i))     	eval([parameters(i,:) '=defaults(i);'])      else    	eval([parameters(i,:) '=options(i);'])      endend% Compute a few iterations of difference equationsif size(u)==[1 1]    un = u(ones(1,N+skip));elseif size(u) == [2 1]	    if ABCD(n+1,n+1) ~= 0	% Require D1=0    	fprintf('%s: Limitation. D1 must be zero for u-ranges.\n', mfilename);    	return;    end    % Make 90% of the u-values take on the extremes    un = uvar(u,skip+N);else    fprintf('%s: Error. Argument 2 (u) has the wrong dimensions.\n', mfilename);    returnend[v x] = simulateDSM(un,ABCD,nlev);x = x(:,skip+1:skip+N);xmin = min(x(1,:)); xmax = max(x(1,:)); dx = xmax-xmin;ymin = min(x(2,:)); ymax = max(x(2,:)); dy = ymax-ymin;axis1 = [ xmin-dx/4 xmax+dx/4 ymin-dy/4 ymax+dy/4 ];% Take the convex hull of the results = hull2d(x')';ec = mean(x')';for i = 1:itnLimit    % Inflate the hull    shift = ec(:,ones(1,length(s)));    s = shift + (1+expFactor)*( s-shift );    % Split the set    [splus eplus sminus eminus] = dssplit2d(u,ABCD,s);    % Map the two halves     s1 = dsmap(u,ABCD,2,splus,eplus,1);    s2 = dsmap(u,ABCD,2,sminus,eminus,-1);    ns = [s1(:,1:size(s1,2)-1) s2(:,1:size(s2,2)-1)];    % Test for inclusion: ns inside s (the inflated hull)    out = outconvex2d(ns,s);    if dbg    	clf; hold on; grid; axis(axis1);        set(gca,'drawmode','fast');    	dotplot(x,'k.');    	dotplot(ec,'o');    	polyplot(s);    	polyplot(s1,'m');    	polyplot(s2,'c');    	outi = logical(sign(out));    	dotplot(ns(:,outi),'rs');    	str = sprintf('Iteration %d: %d image vertices outside',i, sum(outi));    	title(str);        drawnow;    end    if out == 0    	break;    end    % take the hull and repeat    s = hull2d(ns')';end