function [r,p,d,Rrpd,rsd,psd,dsd] = rpdcov(A,B,RAB)

% RPDCOV  Rational fraction to residue-pole covariance transformation
%	[R,P,D] = RPDCOV(A,B) finds the residues R, the poles P
%	and the division rest D of the rational fraction B/A.
%	[R,P,D,RRPD] = RPDCOV(A,B,RAB) transforms the coefficient
%	covariance matrix RAB in the residue-pole-rest covariance
%	matrix RRPD. If A(1) is exaclty 1 the first line and column in
%	RAB can be ommited.
%	[R,P,D,RRPD,RSD,PSD,DSD] = RPDCOV(A,B,RAB) also returns the 
%	standard deviations.
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

nA = length(A) - 1;
[nu,nB] = size(B); nB = nB - 1; nd = max(nB-nA+1,0);

if nA == 0 % bug in residue(1,1)
 r = []; p = []; d = B/A;
else
 r = zeros(nA,nu); d = zeros(nu,nd);
 for ku = 1:nu
  if nd
   [r(:,ku),p,d(ku,:)] = residue(B(ku,:),A);
  else
   [r(:,ku),p] = residue(B(ku,:),A);
  end
 end
end

if nargin > 2
 normalised = length(RAB) == nA + nu*(nB+1);
 Rrpd = zeros(nA+nu*(nA+nd));
 for ku = 1:nu
  kAB = [1:nA nA+1+(ku-1)*(nB+1)+(0:nB)];
  if ~normalised, kAB = [1 kAB+1]; end
  krpd = [(ku-1)*nA+(1:nA) nu*nA+(1:nA) (nu+1)*nA+(ku-1)*nd+(1:nd)];
  if nA == 0
   J = rpdjacobian(A,B(ku,:),[],[],d(ku,:),normalised);
  elseif nd
   J = rpdjacobian(A,B(ku,:),r(:,ku),p,d(ku,:),normalised);
  else
   J = rpdjacobian(A,B(ku,:),r(:,ku),p,[],normalised);
  end
  Rrpd(krpd,krpd) = J * RAB(kAB,kAB) * J.';
 end
end

if nargout > 4 % should be improved for complex numbers
 [J,J,sd] = cpaircov([r(:);p;d(:)],Rrpd);
 rsd = reshape(sd(1:nA*nu),r);
 psd = sd(nA*nu+(1:nA));
 dsd = sd(nA*nu+nA+(1:nd*nu));
end
