Code covered by the BSD License  

Highlights from
slatec

from slatec by Ben Barrowes
The slatec library converted into matlab functions.

[z,fnu,kode,ikflg,n,y,nuf,tol,elim,alim]=cuoik(z,fnu,kode,ikflg,n,y,nuf,tol,elim,alim);
function [z,fnu,kode,ikflg,n,y,nuf,tol,elim,alim]=cuoik(z,fnu,kode,ikflg,n,y,nuf,tol,elim,alim);
%***BEGIN PROLOGUE  CUOIK
%***SUBSIDIARY
%***PURPOSE  Subsidiary to CBESH, CBESI and CBESK
%***LIBRARY   SLATEC
%***TYPE      ALL (CUOIK-A, ZUOIK-A)
%***AUTHOR  Amos, D. E., (SNL)
%***DESCRIPTION
%
%     CUOIK COMPUTES THE LEADING TERMS OF THE UNIFORM ASYMPTOTIC
%     EXPANSIONS FOR THE I AND K FUNCTIONS AND COMPARES THEM
%     (IN LOGARITHMIC FORM) TO ALIM AND ELIM FOR OVER AND UNDERFLOW
%     WHERE ALIM.LT.ELIM. IF THE MAGNITUDE, BASED ON THE LEADING
%     EXPONENTIAL, IS LESS THAN ALIM OR GREATER THAN -ALIM, THEN
%     THE RESULT IS ON SCALE. IF NOT, THEN A REFINED TEST USING OTHER
%     MULTIPLIERS (IN LOGARITHMIC FORM) IS MADE BASED ON ELIM. HERE
%     EXP(-ELIM)=SMALLEST MACHINE NUMBER*1.0E+3 AND EXP(-ALIM)=
%     EXP(-ELIM)/TOL
%
%     IKFLG=1 MEANS THE I SEQUENCE IS TESTED
%          =2 MEANS THE K SEQUENCE IS TESTED
%     NUF = 0 MEANS THE LAST MEMBER OF THE SEQUENCE IS ON SCALE
%         =-1 MEANS AN OVERFLOW WOULD OCCUR
%     IKFLG=1 AND NUF.GT.0 MEANS THE LAST NUF Y VALUES WERE SET TO ZERO
%             THE FIRST N-NUF VALUES MUST BE SET BY ANOTHER ROUTINE
%     IKFLG=2 AND NUF.EQ.N MEANS ALL Y VALUES WERE SET TO ZERO
%     IKFLG=2 AND 0.LT.NUF.LT.N NOT CONSIDERED. Y MUST BE SET BY
%             ANOTHER ROUTINE
%
%***SEE ALSO  CBESH, CBESI, CBESK
%***ROUTINES CALLED  CUCHK, CUNHJ, CUNIK, R1MACH
%***REVISION HISTORY  (YYMMDD)
%   830501  DATE WRITTEN
%   910415  Prologue converted to Version 4.0 format.  (BAB)
%***end PROLOGUE  CUOIK
persistent aarg aic aphi arg ascle asum ax ay bsum cwrk cz czero firstCall fnn gnn gnu i iform init nn nw phi rcz summlv x yy zb zeta1 zeta2 zn zr ; if isempty(firstCall),firstCall=1;end; 

if isempty(arg), arg=0; end;
if isempty(asum), asum=0; end;
if isempty(bsum), bsum=0; end;
if isempty(cwrk), cwrk=zeros(1,16); end;
if isempty(cz), cz=0; end;
if isempty(czero), czero=0; end;
if isempty(phi), phi=0; end;
if isempty(summlv), summlv=0; end;
if isempty(zb), zb=0; end;
if isempty(zeta1), zeta1=0; end;
if isempty(zeta2), zeta2=0; end;
if isempty(zn), zn=0; end;
if isempty(zr), zr=0; end;
if isempty(aarg), aarg=0; end;
if isempty(aic), aic=0; end;
if isempty(aphi), aphi=0; end;
if isempty(ascle), ascle=0; end;
if isempty(ax), ax=0; end;
if isempty(ay), ay=0; end;
if isempty(fnn), fnn=0; end;
if isempty(gnn), gnn=0; end;
if isempty(gnu), gnu=0; end;
if isempty(rcz), rcz=0; end;
if isempty(x), x=0; end;
if isempty(yy), yy=0; end;
if isempty(i), i=0; end;
if isempty(iform), iform=0; end;
if isempty(init), init=0; end;
if isempty(nn), nn=0; end;
if isempty(nw), nw=0; end;
if firstCall,   czero=[complex(0.0e0,0.0e0)];  end;
if firstCall,   aic=[1.265512123484645396e+00];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  CUOIK
nuf = 0;
nn = fix(n);
x = real(z);
zr = z;
if( x<0.0e0 )
zr = -z;
end;
zb = zr;
yy = imag(zr);
ax = abs(x).*1.7321e0;
ay = abs(yy);
iform = 1;
if( ay>ax )
iform = 2;
end;
gnu = max(fnu,1.0e0);
if( ikflg~=1 )
fnn = nn;
gnn = fnu + fnn - 1.0e0;
gnu = max(gnn,fnn);
end;
%-----------------------------------------------------------------------
%     ONLY THE MAGNITUDE OF ARG AND PHI ARE NEEDED ALONG WITH THE
%     REAL PARTS OF ZETA1, ZETA2 AND ZB. NO ATTEMPT IS MADE TO GET
%     THE SIGN OF THE IMAGINARY PART CORRECT.
%-----------------------------------------------------------------------
if( iform==2 )
zn = -zr.*complex(0.0e0,1.0e0);
if( yy<=0.0e0 )
zn = conj(-zn);
end;
[zn,gnu,dumvar3,tol,phi,arg,zeta1,zeta2,asum,bsum]=cunhj(zn,gnu,1,tol,phi,arg,zeta1,zeta2,asum,bsum);
cz = -zeta1 + zeta2;
aarg = abs(arg);
else;
init = 0;
[zr,gnu,ikflg,dumvar4,tol,init,phi,zeta1,zeta2,summlv,cwrk]=cunik(zr,gnu,ikflg,1,tol,init,phi,zeta1,zeta2,summlv,cwrk);
cz = -zeta1 + zeta2;
end;
if( kode==2 )
cz = cz - zb;
end;
if( ikflg==2 )
cz = -cz;
end;
aphi = abs(phi);
rcz = real(cz);
%-----------------------------------------------------------------------
%     OVERFLOW TEST
%-----------------------------------------------------------------------
if( rcz>elim )
nuf = -1;
else;
while (1);
if( rcz<alim )
%-----------------------------------------------------------------------
%     UNDERFLOW TEST
%-----------------------------------------------------------------------
if( rcz>=(-elim) )
if( rcz>(-alim) )
break;
end;
rcz = rcz + log(aphi);
if( iform==2 )
rcz = rcz - 0.25e0.*log(aarg) - aic;
end;
if( rcz>(-elim) )
ascle = 1.0e3.*r1mach(1)./tol;
cz = cz + log(phi);
if( iform~=1 )
cz = cz - complex(0.25e0,0.0e0).*log(arg)- complex(aic,0.0e0);
end;
ax = exp(rcz)./tol;
ay = imag(cz);
cz = complex(ax,0.0e0).*complex(cos(ay),sin(ay));
[cz,nw,ascle,tol]=cuchk(cz,nw,ascle,tol);
if( nw~=1 )
break;
end;
end;
end;
for i = 1 : nn;
y(i) = czero;
end; i = fix(nn+1);
nuf = fix(nn);
return;
else;
rcz = rcz + log(aphi);
if( iform==2 )
rcz = rcz - 0.25e0.*log(aarg) - aic;
end;
if( rcz>elim )
nuf = -1;
return;
end;
end;
break;
end;
if( ikflg~=2 )
if( n~=1 )
%-----------------------------------------------------------------------
%     SET UNDERFLOWS ON I SEQUENCE
%-----------------------------------------------------------------------
while( true );
gnu = fnu +(nn-1);
if( iform==2 )
[zn,gnu,dumvar3,tol,phi,arg,zeta1,zeta2,asum,bsum]=cunhj(zn,gnu,1,tol,phi,arg,zeta1,zeta2,asum,bsum);
cz = -zeta1 + zeta2;
aarg = abs(arg);
else;
init = 0;
[zr,gnu,ikflg,dumvar4,tol,init,phi,zeta1,zeta2,summlv,cwrk]=cunik(zr,gnu,ikflg,1,tol,init,phi,zeta1,zeta2,summlv,cwrk);
cz = -zeta1 + zeta2;
end;
if( kode==2 )
cz = cz - zb;
end;
aphi = abs(phi);
rcz = real(cz);
if( rcz>=(-elim) )
if( rcz>(-alim) )
break;
end;
rcz = rcz + log(aphi);
if( iform==2 )
rcz = rcz - 0.25e0.*log(aarg) - aic;
end;
if( rcz>(-elim) )
ascle = 1.0e3.*r1mach(1)./tol;
cz = cz + log(phi);
if( iform~=1 )
cz = cz - complex(0.25e0,0.0e0).*log(arg)- complex(aic,0.0e0);
end;
ax = exp(rcz)./tol;
ay = imag(cz);
cz = complex(ax,0.0e0).*complex(cos(ay),sin(ay));
[cz,nw,ascle,tol]=cuchk(cz,nw,ascle,tol);
if( nw~=1 )
break;
end;
end;
end;
y(nn) = czero;
nn = fix(nn - 1);
nuf = fix(nuf + 1);
if( nn==0 )
break;
end;
end;
end;
end;
end;
end %subroutine cuoik
%DECK CV

Contact us at files@mathworks.com