function [J,Jinv] = rpdjacobian(A,B,r,p,d,normalised)

% RPDJACOBIAN  Rational fraction to residue-pole-rest Jacobian
%	RPDJACOBIAN(A,B,R,P,D) returns the Jacobian of the mapping between the
%	polynomial coefficient set [A B] of the rational fraction B/A and its 
%	residues, poles and division rest [R P D]. The inverse Jacobian is
%	also returned as second output. A non zero additional input argument 
%	indicates that A(1) is exactly 1 and that the first column
%	in the Jacobian can be omitted.
%
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.


if nargin < 6, normalised = 0; end
M = length(B) - 1; N = length(p); K = length(d);

Q = zeros(N,N); QQ = zeros(N,2*N-1);
for k = 1:N
 Q(k,:) = A(1) * poly(p([1:k-1,k+1:N]));
 QQ(k,:) = r(k) * conv(Q(k,:),Q(k,:));
end

TB = toepconv(B,N);
if normalised, TB(1,:) = []; end

CAB = [zeros(N+M+1+(~normalised),max(N-M-1,0)) [-TB;toepconv(A,M)]];
CRPD = [[zeros(N,K) Q*toepconv(A,N-1)];[zeros(N,K+1) QQ];toepconv(conv(A,A),K-1)];

J = CRPD' \ CAB';
if nargout > 1, Jinv = CAB' \ CRPD'; end