function zv=rotvar(pol,covm)
%ROTVAR	computes the standard deviations of roots to polynomials
%
%	ZV = rotvar(POL,COVM)
%
%	POL is the entered polynomial and
%	COVM is the covariance matrix of its coefficients
%	ZV is returned as the roots of POL and their standard deviations.
%	The first column of ZV contains the roots and the corresponding
%	elements in the second column give the standard deviations.
%	For complex conjugate pairs the second entry is the correlation
%	coefficient between the real and imaginary parts.
%	
%	The routine is primarily a help subroutine to th2zp.

%	L. Ljung 87-7-2
%	Copyright (c) 1987 by the MathWorks, Inc.
%	All Rights Reserved.

[dr,dc]=size(covm);lp=length(pol);
r=roots(pol);
zv(:,1)=r;lr=length(r);
if isempty(covm),zv(:,2)=zeros(lr,1);return,end
k=1;
while k<lr+1
	if imag(r(k))==0
		D(k,:)=real(-poly(r([1:k-1,k+1:lr])));
		k=k+1;
	else
	
		D(k,:)=real(-2*poly([r([1:k-1,k+2:lr]);real(r(k))]));
		D(k+1,:)=real([0, 2*imag(r(k))*poly(r([1:k-1,k+2:lr]))]);
		k=k+2;
	end
end
D=inv(D)';
if lr~=dr,
     DZ(:,2:lp)=D/pol(1);
     DZ(:,1)=-D*pol(2:lp)'/(pol(1)^2);
else
     DZ=D;
end
PV=DZ*covm*DZ';
k=1;l=1;i=sqrt(-1);
while k<lp
     %if r(k)==0,zv(k,2)=0;k=k+1;else 
     if imag(r(k))==0,
     zv(k,2)=sqrt(PV(l,l));k=k+1;l=l+1;
          else zv(k,2)=sqrt(PV(l,l))+i*sqrt(PV(l+1,l+1));
          zv(k+1,2)=PV(l,l+1)/sqrt(PV(l,l)*PV(l+1,l+1));
	  k=k+2;l=l+2;
          end
     end
end
