function[Q,Lstar,As,Bs]=polysolve(S,CCstar,B,N,rBeta,DT,m,k)
%
% function[Q,Lstar,As,Bs]=polysolve(S,CCstar,B,N,rBeta,DT,m,k)
%
% Solves the polynomial equation
%     -m+k
%    z    S(CC_*)B_*N_* = rBeta_* Q + z(DT)L_*
%
%  with respect to Q and L_* of degree
%
%   nQ=max(ns+nc+m-k,nd+nt-1), nL=max(nc+nb+nn-m+k,nbta)-1
%
%  The polynomials may have complex coefficients.
%  In the input, CCstar = C(z^-1)*C_*(z)
%                rBeta  = r*Beta(z^-1)
%                m,k can have any sign.  (k is transducer delay)
%  On exit, Q and L_* polynomials are delivered, as well as
%  linear system As*x = Bs.
%  USES: sylv
%
%  Author: Anders Ahlen  Revised: Mikael Sternad
%  1988-10-13            1993-05-12
%
ns=length(S)-1;       nb=length(B)-1;
nc=(length(CCstar)-1)/2;  nd=length(DT)-1;
nn=length(N)-1;       BTA =rBeta ;
nbta=length(BTA)-1;   BN = conv(B,N);

nal=nc+nb+nn-m+k;     x=[nal nbta];
nll=max(x);           y=[ns+nc+m-k,nd-1];
nQ=max(y);
%
if nQ==-1
  Q=0; return;
end
%
nrow=nQ+nll+1;                hl=zeros(1,nrow);
ni=ns+nc+nb+nn+nc;            CCBBN=zeros(1,nrow);
BBN=conj(BN(nb+nn+1:-1:1));   CCBBN=conv(BBN,conv(CCstar,S));
%
  if nal>=nbta
    BBTA=conj(BTA(nbta+1:-1:1));
    BBTA=[zeros(1,nal-nbta) BBTA];
    s=sylv(DT,nll,BBTA,nQ+1,nrow);
    hl(1:ni+1)=CCBBN;
    th=s\hl';
  else
    BBTA=conj(BTA(nbta+1:-1:1));
    CCBBN=[zeros(1,nbta-nal) CCBBN];
    s=sylv(DT,nll,BBTA,nQ+1,nrow);
    hl(1:ni+nbta-nal+1)=CCBBN;
    th=s\hl';
  end
Q=th(nll+1:nll+nQ+1,1)';
LL=th(1:nll,1)';        Lstar=LL(nll:-1:1);
Bs=hl(nrow:-1:1)';      As=s(nrow:-1:1,nrow:-1:1);