[difference,U0]=twirl2(rho,varargin);

% twirl2   How much a state changes under unitaries of the form mkron(U,U,U,...)
%   twirl2(rho) gives the maximal difference between the original state,
%   rho and the state obtained from it by a unitary of the form
%   mkron(U,U,U,...), where U is a single qubit unitary.
%   The difference is computed through the norm 
%   ||A||=sum_kl |A_kl|^2. The difference 
%   is zero for Werner states. The form twirl(rho,d) makes it
%   possible to twirl a register of qudits with dimension d. 
%   Using the form twirl(rho,d,Nit) we can determine how many random
%   unitaries are used for twirling. The default value for Nit is 100.
%   The form  [difference,U0]=twirl2(rho) gives also back
%   the unitary U0 for which the difference is the largest
%   between the original and the rotated state.

function [difference,U0]=twirl2(rho,varargin);

if length(varargin)==0,
    % Dimension of quidits
    d=2;
    % Number of random unitaries used
    Nit=100;
elseif length(varargin)==1,
    d=varargin{1};
    Nit=100;
elseif length(varargin)==2,
    d=varargin{1};   
    Nit=varargin{2};
else
    error('Wrong number of input arguments');
end %if
   
x=[0 1;1 0];
z=[1 0;0 -1];
y=i*x*z;

[sy,sx]=size(rho);
N=log2(sx)/log2(d); 

difference=0;
U=zeros(d,d);
for n=1:Nit     
    
    % Create a random dxd unitary
    % from d orthogonal vectors
    for k=1:d
        vv=randn(d,1)+i*randn(d,1); 
        for m=1:k-1
            vv=vv-U(:,m)*(U(:,m)'*vv);
        end %for
        U(:,k)=vv/sqrt(vv'*vv);
    end %for 
    
    UU=U;
    for n=2:N
        UU=kron(UU,U);   
    end %for
    r=UU*rho*UU';
    % real() is important since MATLAB gives results with small
    % imaginary part; this spoils the use of > and <
    dd=real(trace((r-rho)*(r-rho)'));
    if dd>difference,
        difference=dd;
        U0=U;
    end %if
end %for