[n,m,a,z,c,b,kbd]=cpevl(n,m,a,z,c,b,kbd);

function [n,m,a,z,c,b,kbd]=cpevl(n,m,a,z,c,b,kbd);
persistent bi bim1 ci cim1 d1 firstCall i j mini np1 q r s t ; if isempty(firstCall),firstCall=1;end; 

if isempty(d1), d1=0; end;
if isempty(r), r=0; end;
if isempty(s), s=0; end;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(mini), mini=0; end;
if isempty(np1), np1=0; end;
%***BEGIN PROLOGUE  CPEVL
%***SUBSIDIARY
%***PURPOSE  Subsidiary to CPZERO
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (CPEVL-S)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%        Evaluate a complex polynomial and its derivatives.
%        Optionally compute error bounds for these values.
%
%   INPUT...
%        N = Degree of the polynomial
%        M = Number of derivatives to be calculated,
%            M=0 evaluates only the function
%            M=1 evaluates the function and first derivative, etc.
%             if M .GT. N+1 function and all N derivatives will be
%                calculated.
%       A = Complex vector containing the N+1 coefficients of polynomial
%               A(I)= coefficient of Z**(N+1-I)
%        Z = Complex point at which the evaluation is to take place.
%        C = Array of 2(M+1) words into which values are placed.
%        B = Array of 2(M+1) words only needed if bounds are to be
%              calculated.  It is not used otherwise.
%        KBD = A logical variable, e.g. true or false which is
%              to be set true if bounds are to be computed.
%
%  OUTPUT...
%        C =  C(I+1) contains the complex value of the I-th
%              derivative at Z, I=0,...,M
%        B =  B(I) contains the bounds on the real and imaginary parts
%              of C(I) if they were requested.
%
%***SEE ALSO  CPZERO
%***ROUTINES CALLED  I1MACH
%***REVISION HISTORY  (YYMMDD)
%   810223  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900402  Added TYPE section.  (WRB)
%***end PROLOGUE  CPEVL
%
a_shape=size(a);a=reshape(a,1,[]);
c_shape=size(c);c=reshape(c,1,[]);
if isempty(ci), ci=0; end;
if isempty(cim1), cim1=0; end;
b_shape=size(b);b=reshape(b,1,[]);
if isempty(bi), bi=0; end;
if isempty(bim1), bim1=0; end;
if isempty(t), t=0; end;
% za= @(q)  cmplx(abs(real(q)),abs(aimag(q)));complex :: za ;
if isempty(q), q=0; end;
if firstCall,   d1=[0.0];  end;
firstCall=0;
za= @(q)  complex(abs(real(q)),abs(imag(q)));
%***FIRST EXECUTABLE STATEMENT  CPEVL
if( d1==0.0 )
d1 = real(i1mach(10)).^(1-i1mach(11));
end;
np1 = fix(n + 1);
for j = 1 : np1;
ci = 0.0;
cim1 = a(j);
bi = 0.0;
bim1 = 0.0;
mini = fix(min(m+1,n+2-j));
for i = 1 : mini;
if( j~=1 )
ci = c(i);
end;
if( i~=1 )
cim1 = c(i-1);
end;
c(i) = cim1 + z.*ci;
if( kbd )
if( j~=1 )
bi = b(i);
end;
if( i~=1 )
bim1 = b(i-1);
end;
za= @(q)  complex(abs(real(q)),abs(imag(q)));
t = bi +(3..*d1+4..*d1.*d1).*za(ci);
za= @(q)  complex(abs(real(q)),abs(imag(q)));
r = real(za(z).*complex(real(t),-imag(t)));
za= @(q)  complex(abs(real(q)),abs(imag(q)));
s = imag(za(z).*t);
za= @(q)  complex(abs(real(q)),abs(imag(q)));
b(i) =(1.+8..*d1).*(bim1+d1.*za(cim1)+complex(r,s));
if( j==1 )
b(i) = 0.0;
end;
end;
end; i = fix(mini+1);
end; j = fix(np1+1);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
end
%DECK CPEVLR