[in,a,r,t,iflg,s]=cpzero(in,a,r,t,iflg,s);

function [in,a,r,t,iflg,s]=cpzero(in,a,r,t,iflg,s);
persistent i imax j n n1 nit nmax nr pn temp u v x ; 

if isempty(i), i=0; end;
if isempty(imax), imax=0; end;
if isempty(j), j=0; end;
if isempty(n), n=0; end;
if isempty(n1), n1=0; end;
if isempty(nit), nit=0; end;
if isempty(nmax), nmax=0; end;
if isempty(nr), nr=0; end;
if isempty(u), u=0; end;
if isempty(v), v=0; end;
if isempty(x), x=0; end;
%***BEGIN PROLOGUE  CPZERO
%***PURPOSE  Find the zeros of a polynomial with complex coefficients.
%***LIBRARY   SLATEC
%***CATEGORY  F1A1B
%***TYPE      COMPLEX (RPZERO-S, CPZERO-C)
%***KEYWORDS  POLYNOMIAL ROOTS, POLYNOMIAL ZEROS, REAL ROOTS
%***AUTHOR  Kahaner, D. K., (NBS)
%***DESCRIPTION
%
%      Find the zeros of the complex polynomial
%         P(Z)= A(1)*Z**N + A(2)*Z**(N-1) +...+ A(N+1)
%
%    Input...
%       IN = degree of P(Z)
%       A = complex vector containing coefficients of P(Z),
%            A(I) = coefficient of Z**(N+1-i)
%       R = N word complex vector containing initial estimates for zeros
%            if these are known.
%       T = 4(N+1) word array used for temporary storage
%       IFLG = flag to indicate if initial estimates of
%              zeros are input.
%            If IFLG .EQ. 0, no estimates are input.
%            If IFLG .NE. 0, the vector R contains estimates of
%               the zeros
%       ** WARNING ****** If estimates are input, they must
%                         be separated, that is, distinct or
%                         not repeated.
%       S = an N word array
%
%    Output...
%       R(I) = Ith zero,
%       S(I) = bound for R(I) .
%       IFLG = error diagnostic
%    Error Diagnostics...
%       If IFLG .EQ. 0 on return, all is well
%       If IFLG .EQ. 1 on return, A(1)=0.0 or N=0 on input
%       If IFLG .EQ. 2 on return, the program failed to converge
%                after 25*N iterations.  Best current estimates of the
%                zeros are in R(I).  Error bounds are not calculated.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  CPEVL
%***REVISION HISTORY  (YYMMDD)
%   810223  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%***end PROLOGUE  CPZERO
%
s_shape=size(s);s=reshape(s,1,[]);
r_shape=size(r);r=reshape(r,1,[]);
t_shape=size(t);t=reshape(t,1,[]);
a_shape=size(a);a=reshape(a,1,[]);
if isempty(pn), pn=0; end;
if isempty(temp), temp=0; end;
%***FIRST EXECUTABLE STATEMENT  CPZERO
if( in<=0 || abs(a(1))==0.0 )
iflg = 1;
else;
%
%       CHECK FOR EASILY OBTAINED ZEROS
%
n = fix(in);
n1 = fix(n + 1);
if( iflg==0 )
while( true );
n1 = fix(n + 1);
if( n<=1 )
r(1) = -a(2)./a(1);
s(1) = 0.0;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
return;
else;
if( abs(a(n1))~=0.0 )
break;
end;
r(n) = 0.0;
s(n) = 0.0;
n = fix(n - 1);
end;
end;
%
%          IF INITIAL ESTIMATES FOR ZEROS NOT GIVEN, FIND SOME
%
temp = -a(2)./(a(1).*n);
n_orig=n; t_orig=t;    [n,dumvar2,a,temp,t,dumvar6]=cpevl(n,n,a,temp,t,t,false);    t(dumvar6~=t_orig)=dumvar6(dumvar6~=t_orig); n(dumvar2~=n_orig)=dumvar2(dumvar2~=n_orig);
imax = fix(n + 2);
t(n1) = abs(t(n1));
for i = 2 : n1;
t(n+i) = -abs(t(n+2-i));
if( real(t(n+i))<real(t(imax)) )
imax = fix(n + i);
end;
end; i = fix(n1+1);
x =(-real(t(imax))./real(t(n1))).^(1../(imax-n1));
while( true );
x = 2..*x;
pn_orig=pn;    [n,dumvar2,t(sub2ind(size(t),max(n1,1)):end),dumvar4,pn,dumvar6]=cpevl(n,0,t(sub2ind(size(t),max(n1,1)):end),complex(x,0.0),pn,pn,false);    pn(dumvar6~=pn_orig)=dumvar6(dumvar6~=pn_orig);
if( real(pn)>=0. )
break;
end;
end;
u = .5.*x;
v = x;
while( true );
x = .5.*(u+v);
pn_orig=pn;    [n,dumvar2,t(sub2ind(size(t),max(n1,1)):end),dumvar4,pn,dumvar6]=cpevl(n,0,t(sub2ind(size(t),max(n1,1)):end),complex(x,0.0),pn,pn,false);    pn(dumvar6~=pn_orig)=dumvar6(dumvar6~=pn_orig);
if( real(pn)>0. )
v = x;
end;
if( real(pn)<=0. )
u = x;
end;
if((v-u)<=.001.*(1.+v) )
break;
end;
end;
for i = 1 : n;
u =(3.14159265./n).*(2.*i-1.5);
r(i) = max(x,.001.*abs(temp)).*complex(cos(u),sin(u)) + temp;
end; i = fix(n+1);
end;
%
%          MAIN ITERATION LOOP STARTS HERE
%
nr = 0;
nmax = fix(25.*n);
for nit = 1 : nmax;
for i = 1 : n;
if( nit==1 || abs(t(i))~=0. )
[n,dumvar2,a,r(i),pn,temp]=cpevl(n,0,a,r(i),pn,temp,true);
if( abs(real(pn))+abs(imag(pn))>real(temp)+imag(temp) )
temp = a(1);
for j = 1 : n;
if( j~=i )
temp = temp.*(r(i)-r(j));
end;
end; j = fix(n+1);
t(i) = pn./temp;
else;
t(i) = 0.0;
nr = fix(nr + 1);
end;
end;
end; i = fix(n+1);
for i = 1 : n;
r(i) = r(i) - t(i);
end; i = fix(n+1);
if( nr==n )
%
%          CALCULATE ERROR BOUNDS FOR ZEROS
%
for nr = 1 : n;
n_orig=n;    [n,dumvar2,a,r(nr),t,dumvar6]=cpevl(n,n,a,r(nr),t,t(sub2ind(size(t),max(n+2,1)):end),true);    n(dumvar2~=n_orig)=dumvar2(dumvar2~=n_orig);   dumvar6i=find((t(sub2ind(size(t),max(n+2,1)):end))~=(dumvar6));   t(n+2-1+dumvar6i)=dumvar6(dumvar6i); 
x = abs(complex(abs(real(t(1))),abs(imag(t(1))))+t(n+2));
s(nr) = 0.0;
for i = 1 : n;
x = x.*real(n1-i)./i;
temp = complex(max(abs(real(t(i+1)))-real(t(n1+i)),0.0),max(abs(imag(t(i+1)))-imag(t(n1+i)),0.0));
s(nr) = max(s(nr),(abs(temp)./x).^(1../i));
end; i = fix(n+1);
s(nr) = 1../s(nr);
end; nr = fix(n+1);
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
return;
end;
end; nit = fix(nmax+1);
%        ERROR EXITS
iflg = 2;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
return;
end;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
end %subroutine cpzero
%DECK CQRDC