[zr,zi,fnu,kode,n,yr,yi,nz,nlast,fnul,tol,elim,alim]=zuni1(zr,zi,fnu,kode,n,yr,yi,nz,nlast,fnul,tol,elim,alim);

function [zr,zi,fnu,kode,n,yr,yi,nz,nlast,fnul,tol,elim,alim]=zuni1(zr,zi,fnu,kode,n,yr,yi,nz,nlast,fnul,tol,elim,alim);
%***BEGIN PROLOGUE  ZUNI1
%***SUBSIDIARY
%***PURPOSE  Subsidiary to ZBESI and ZBESK
%***LIBRARY   SLATEC
%***TYPE      ALL (CUNI1-A, ZUNI1-A)
%***AUTHOR  Amos, D. E., (SNL)
%***DESCRIPTION
%
%     ZUNI1 COMPUTES I(FNU,Z)  BY MEANS OF THE UNIFORM ASYMPTOTIC
%     EXPANSION FOR I(FNU,Z) IN -PI/3.LE.ARG Z.LE.PI/3.
%
%     FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
%     EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
%     NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
%     FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
%     Y(I)=CZERO FOR I=NLAST+1,N
%
%***SEE ALSO  ZBESI, ZBESK
%***ROUTINES CALLED  D1MACH, ZABS, ZUCHK, ZUNIK, ZUOIK
%***REVISION HISTORY  (YYMMDD)
%   830501  DATE WRITTEN
%   910415  Prologue converted to Version 4.0 format.  (BAB)
%***end PROLOGUE  ZUNI1
%     COMPLEX CFN,CONE,CRSC,CSCL,CSR,CSS,CWRK,CZERO,C1,C2,PHI,RZ,SUM,S1,
%    *S2,Y,Z,ZETA1,ZETA2
persistent aphi ascle bry c1r c2i c2m c2r coner crsc cscl csrr cssr cwrki cwrkr cyi cyr firstCall fn i iflag igo init k m nd nn nuf nw phii phir rast rs1 rzi rzr s1i s1r s2i s2r sti str sumi sumr zeroi zeror zeta1i zeta1r zeta2i zeta2r ; if isempty(firstCall),firstCall=1;end; 

if isempty(aphi), aphi=0; end;
if isempty(ascle), ascle=0; end;
if isempty(bry), bry=zeros(1,3); end;
if isempty(coner), coner=0; end;
if isempty(crsc), crsc=0; end;
if isempty(cscl), cscl=0; end;
if isempty(csrr), csrr=zeros(1,3); end;
if isempty(cssr), cssr=zeros(1,3); end;
if isempty(cwrki), cwrki=zeros(1,16); end;
if isempty(cwrkr), cwrkr=zeros(1,16); end;
if isempty(c1r), c1r=0; end;
if isempty(c2i), c2i=0; end;
if isempty(c2m), c2m=0; end;
if isempty(c2r), c2r=0; end;
if isempty(fn), fn=0; end;
if isempty(phii), phii=0; end;
if isempty(phir), phir=0; end;
if isempty(rast), rast=0; end;
if isempty(rs1), rs1=0; end;
if isempty(rzi), rzi=0; end;
if isempty(rzr), rzr=0; end;
if isempty(sti), sti=0; end;
if isempty(str), str=0; end;
if isempty(sumi), sumi=0; end;
if isempty(sumr), sumr=0; end;
if isempty(s1i), s1i=0; end;
if isempty(s1r), s1r=0; end;
if isempty(s2i), s2i=0; end;
if isempty(s2r), s2r=0; end;
if isempty(zeroi), zeroi=0; end;
if isempty(zeror), zeror=0; end;
if isempty(zeta1i), zeta1i=0; end;
if isempty(zeta1r), zeta1r=0; end;
if isempty(zeta2i), zeta2i=0; end;
if isempty(zeta2r), zeta2r=0; end;
if isempty(cyr), cyr=zeros(1,2); end;
if isempty(cyi), cyi=zeros(1,2); end;
if isempty(i), i=0; end;
if isempty(iflag), iflag=0; end;
if isempty(init), init=0; end;
if isempty(k), k=0; end;
if isempty(m), m=0; end;
if isempty(nd), nd=0; end;
if isempty(nn), nn=0; end;
if isempty(nuf), nuf=0; end;
if isempty(nw), nw=0; end;
if isempty(igo), igo=0; end;
if firstCall,   zeror =[0.0d0];  end;
if firstCall,  zeroi =[0.0d0];  end;
if firstCall,  coner=[1.0d0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  ZUNI1
nz = 0;
nd = fix(n);
nlast = 0;
%-----------------------------------------------------------------------
%     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
%     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
%     EXP(ALIM)=EXP(ELIM)*TOL
%-----------------------------------------------------------------------
cscl = 1.0d0./tol;
crsc = tol;
cssr(1) = cscl;
cssr(2) = coner;
cssr(3) = crsc;
csrr(1) = crsc;
csrr(2) = coner;
csrr(3) = cscl;
bry(1) = 1.0d3.*d1mach(1)./tol;
%-----------------------------------------------------------------------
%     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
%-----------------------------------------------------------------------
fn = max(fnu,1.0d0);
init = 0;
[zr,zi,fn,dumvar4,dumvar5,tol,init,phir,phii,zeta1r,zeta1i,zeta2r,zeta2i,sumr,sumi,cwrkr,cwrki]=zunik(zr,zi,fn,1,1,tol,init,phir,phii,zeta1r,zeta1i,zeta2r,zeta2i,sumr,sumi,cwrkr,cwrki);
if( kode==1 )
s1r = -zeta1r + zeta2r;
s1i = -zeta1i + zeta2i;
else;
str = zr + zeta2r;
sti = zi + zeta2i;
rast = fn./zabs(str,sti);
str = str.*rast.*rast;
sti = -sti.*rast.*rast;
s1r = -zeta1r + str;
s1i = -zeta1i + sti;
end;
rs1 = s1r;
igo=0;
if( abs(rs1)<=elim )
while (1);
nn = fix(min(2,nd));
igo=0;
for i = 1 : nn;
fn = fnu +(nd-i);
init = 0;
[zr,zi,fn,dumvar4,dumvar5,tol,init,phir,phii,zeta1r,zeta1i,zeta2r,zeta2i,sumr,sumi,cwrkr,cwrki]=zunik(zr,zi,fn,1,0,tol,init,phir,phii,zeta1r,zeta1i,zeta2r,zeta2i,sumr,sumi,cwrkr,cwrki);
if( kode==1 )
s1r = -zeta1r + zeta2r;
s1i = -zeta1i + zeta2i;
else;
str = zr + zeta2r;
sti = zi + zeta2i;
rast = fn./zabs(str,sti);
str = str.*rast.*rast;
sti = -sti.*rast.*rast;
s1r = -zeta1r + str;
s1i = -zeta1i + sti + zi;
end;
%-----------------------------------------------------------------------
%     TEST FOR UNDERFLOW AND OVERFLOW
%-----------------------------------------------------------------------
rs1 = s1r;
if( abs(rs1)>elim )
igo=1;
break;
end;
if( i==1 )
iflag = 2;
end;
if( abs(rs1)>=alim )
%-----------------------------------------------------------------------
%     REFINE  TEST AND SCALE
%-----------------------------------------------------------------------
[aphi ,phir,phii]=zabs(phir,phii);
rs1 = rs1 + log(aphi);
if( abs(rs1)>elim )
igo=1;
break;
end;
if( i==1 )
iflag = 1;
end;
if( rs1>=0.0d0 )
if( i==1 )
iflag = 3;
end;
end;
end;
%-----------------------------------------------------------------------
%     SCALE S1 IF ABS(S1).LT.ASCLE
%-----------------------------------------------------------------------
s2r = phir.*sumr - phii.*sumi;
s2i = phir.*sumi + phii.*sumr;
str = exp(s1r).*cssr(iflag);
s1r = str.*cos(s1i);
s1i = str.*sin(s1i);
str = s2r.*s1r - s2i.*s1i;
s2i = s2r.*s1i + s2i.*s1r;
s2r = str;
if( iflag==1 )
[s2r,s2i,nw,bry(1),tol]=zuchk(s2r,s2i,nw,bry(1),tol);
if( nw~=0 )
igo=1;
break;
end;
end;
cyr(i) = s2r;
cyi(i) = s2i;
m = fix(nd - i + 1);
yr(m) = s2r.*csrr(iflag);
yi(m) = s2i.*csrr(iflag);
end;
if( nd>2 && igo==0 )
rast = 1.0d0./zabs(zr,zi);
str = zr.*rast;
sti = -zi.*rast;
rzr =(str+str).*rast;
rzi =(sti+sti).*rast;
bry(2) = 1.0d0./bry(1);
[bry(3) ]=d1mach(2);
s1r = cyr(1);
s1i = cyi(1);
s2r = cyr(2);
s2i = cyi(2);
c1r = csrr(iflag);
ascle = bry(iflag);
k = fix(nd - 2);
fn = k;
for i = 3 : nd;
c2r = s2r;
c2i = s2i;
s2r = s1r +(fnu+fn).*(rzr.*c2r-rzi.*c2i);
s2i = s1i +(fnu+fn).*(rzr.*c2i+rzi.*c2r);
s1r = c2r;
s1i = c2i;
c2r = s2r.*c1r;
c2i = s2i.*c1r;
yr(k) = c2r;
yi(k) = c2i;
k = fix(k - 1);
fn = fn - 1.0d0;
if( iflag<3 )
str = abs(c2r);
sti = abs(c2i);
c2m = max(str,sti);
if( c2m>ascle )
iflag = fix(iflag + 1);
ascle = bry(iflag);
s1r = s1r.*c1r;
s1i = s1i.*c1r;
s2r = c2r;
s2i = c2i;
s1r = s1r.*cssr(iflag);
s1i = s1i.*cssr(iflag);
s2r = s2r.*cssr(iflag);
s2i = s2i.*cssr(iflag);
c1r = csrr(iflag);
end;
end;
end; i = fix(nd+1);
end;
if(igo==0)
return;
end;
%-----------------------------------------------------------------------
%     SET UNDERFLOW AND UPDATE PARAMETERS
%-----------------------------------------------------------------------
if( rs1<=0.0d0 )
yr(nd) = zeror;
yi(nd) = zeroi;
nz = fix(nz + 1);
nd = fix(nd - 1);
if( nd==0 )
return;
end;
[zr,zi,fnu,kode,dumvar5,nd,yr,yi,nuf,tol,elim,alim]=zuoik(zr,zi,fnu,kode,1,nd,yr,yi,nuf,tol,elim,alim);
if( nuf>=0 )
nd = fix(nd - nuf);
nz = fix(nz + nuf);
if( nd==0 )
return;
end;
fn = fnu +(nd-1);
if( fn>=fnul )
continue;
end;
nlast = fix(nd);
return;
end;
end;
break;
end;
elseif( rs1<=0.0d0 ) ;
nz = fix(n);
for i = 1 : n;
yr(i) = zeror;
yi(i) = zeroi;
end; i = fix(n+1);
return;
end;
nz = -1;
end %subroutine zuni1
%DECK ZUNI2