function [Deltafull,Deltadiag,fullcount,diagcount] = finddelta(M,iter)

%
%	This short program is designed to provide a numerical verification 
%	of the difference between unstructured and structured uncertainty.
%	Robust control tells us that there exist matrices M for which a 
%	full matrix Delta can lead to an unstable M-Delta pair, but a diagonal 
%	matrix Delta will always give a stable M-Delta pair.
%
%	This program does not PROVE this theorem, but it DOES allow you 
%	to build up evidence by testing many random examples.
%
%	The test for stability comes from the spectral radius condition in 
%	Theorem 8.2 on page 305 in  Skogestad and Postlethwaite.  The 
%	assumptions regarding the uncertainty structure are found there.
%
%	Remark:  Note that it seems that as we go to larger matrices M (say about 
%	4-by-4) that it frequently happens that we can find more destabilizing diagonal 
%	matrices Delta, than full matrices Delta.  Strictly speaking this should 
%	not be possible - diagonal matrices are a subset of full matrices.
%
%	My theory is that the generation of random matrices currently being used
%	below, does not "find" the soft spot in M due to the many degrees of 
%	freedom in a full-matrix Delta.  It may be beneficial to generate Delta 
%	using the svd of M.  In other words 
%
%			Delta = V*diag(...random complex elements ...)*U';
%
%	but that is only a hypothesis and the testing is left as an exercise for 
%	the interested student.  In other words, I don't really know the answer!
%
%		Author:	Greg Stewart
%		Date:		October 21, 2001.
%

%
%	Determine the size of M.
%
n = length(M);

%
%	Search for a FULL complex matrix Delta that will result in an 
%	unstable M-Delta pair.
%
fullcount = 0;
for k = 1:iter
	Delta_temp = rand(n,n) + sqrt(-1)*rand(n,n);
   Delta_temp = Delta_temp/max(svd(Delta_temp));	% ensuring that max(svd(Delta))=1
   if max(abs(eig(M*Delta_temp))) >= 1
      fullcount = fullcount + 1;
   	Deltafull = Delta_temp;
   end
end

%
%	Search for DIAGONAL complex matrix Delta that will result in an 
%	unstable M-Delta pair.
%
diagcount = 0;
for k = 1:iter
	Delta_temp = diag(rand(1,n)) + sqrt(-1)*diag(rand(1,n));
   Delta_temp = Delta_temp/max(svd(Delta_temp));
   if max(abs(eig(M*Delta_temp))) >= 1
      diagcount = diagcount + 1;
      Deltadiag = Delta_temp;
   end
end