Code covered by the BSD License  

Highlights from
slatec

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

[zrr,zri,fnu,n,yr,yi,nz,rzr,rzi,ascle,tol,elim]=zkscl(zrr,zri,fnu,n,yr,yi,nz,rzr,rzi,ascle,tol,elim);
function [zrr,zri,fnu,n,yr,yi,nz,rzr,rzi,ascle,tol,elim]=zkscl(zrr,zri,fnu,n,yr,yi,nz,rzr,rzi,ascle,tol,elim);
%***BEGIN PROLOGUE  ZKSCL
%***SUBSIDIARY
%***PURPOSE  Subsidiary to ZBESK
%***LIBRARY   SLATEC
%***TYPE      ALL (CKSCL-A, ZKSCL-A)
%***AUTHOR  Amos, D. E., (SNL)
%***DESCRIPTION
%
%     SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE
%     ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN
%     RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL.
%
%***SEE ALSO  ZBESK
%***ROUTINES CALLED  ZABS, ZLOG, ZUCHK
%***REVISION HISTORY  (YYMMDD)
%   830501  DATE WRITTEN
%   910415  Prologue converted to Version 4.0 format.  (BAB)
%   930122  Added ZLOG to EXTERNAL statement.  (RWC)
%***end PROLOGUE  ZKSCL
%     COMPLEX CK,CS,CY,CZERO,RZ,S1,S2,Y,ZR,ZD,CELM
persistent acs alas as celmr cki ckr csi csr cyi cyr elm firstCall fn helim i ic idum igo kk nn nw s1i s1r s2i s2r str zdi zdr zeroi zeror ; if isempty(firstCall),firstCall=1;end; 

if isempty(acs), acs=0; end;
if isempty(as), as=0; end;
if isempty(cki), cki=0; end;
if isempty(ckr), ckr=0; end;
if isempty(csi), csi=0; end;
if isempty(csr), csr=0; end;
if isempty(cyi), cyi=zeros(1,2); end;
if isempty(cyr), cyr=zeros(1,2); end;
if isempty(fn), fn=0; end;
if isempty(str), str=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(zdr), zdr=0; end;
if isempty(zdi), zdi=0; end;
if isempty(celmr), celmr=0; end;
if isempty(elm), elm=0; end;
if isempty(helim), helim=0; end;
if isempty(alas), alas=0; end;
if isempty(i), i=0; end;
if isempty(ic), ic=0; end;
if isempty(idum), idum=0; end;
if isempty(kk), kk=0; end;
if isempty(nn), nn=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;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  ZKSCL
nz = 0;
ic = 0;
nn = fix(min(2,n));
for i = 1 : nn;
s1r = yr(i);
s1i = yi(i);
cyr(i) = s1r;
cyi(i) = s1i;
[as ,s1r,s1i]=zabs(s1r,s1i);
acs = -zrr + log(as);
nz = fix(nz + 1);
yr(i) = zeror;
yi(i) = zeroi;
if( acs>=(-elim) )
[s1r,s1i,csr,csi,idum]=zlog(s1r,s1i,csr,csi,idum);
csr = csr - zrr;
csi = csi - zri;
str = exp(csr)./tol;
csr = str.*cos(csi);
csi = str.*sin(csi);
[csr,csi,nw,ascle,tol]=zuchk(csr,csi,nw,ascle,tol);
if( nw==0 )
yr(i) = csr;
yi(i) = csi;
ic = fix(i);
nz = fix(nz - 1);
end;
end;
end; i = fix(nn+1);
if( n==1 )
return;
end;
if( ic<=1 )
yr(1) = zeror;
yi(1) = zeroi;
nz = 2;
end;
if( n==2 )
return;
end;
if( nz==0 )
return;
end;
fn = fnu + 1.0d0;
ckr = fn.*rzr;
cki = fn.*rzi;
s1r = cyr(1);
s1i = cyi(1);
s2r = cyr(2);
s2i = cyi(2);
helim = 0.5d0.*elim;
elm = exp(-elim);
celmr = elm;
zdr = zrr;
zdi = zri;
%
%     FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
%     S2 GETS LARGER THAN EXP(ELIM/2)
%
igo=0;
for i = 3 : n;
kk = fix(i);
csr = s2r;
csi = s2i;
s2r = ckr.*csr - cki.*csi + s1r;
s2i = cki.*csr + ckr.*csi + s1i;
s1r = csr;
s1i = csi;
ckr = ckr + rzr;
cki = cki + rzi;
[as ,s2r,s2i]=zabs(s2r,s2i);
alas = log(as);
acs = -zdr + alas;
nz = fix(nz + 1);
yr(i) = zeror;
yi(i) = zeroi;
if( acs>=(-elim) )
[s2r,s2i,csr,csi,idum]=zlog(s2r,s2i,csr,csi,idum);
csr = csr - zdr;
csi = csi - zdi;
str = exp(csr)./tol;
csr = str.*cos(csi);
csi = str.*sin(csi);
[csr,csi,nw,ascle,tol]=zuchk(csr,csi,nw,ascle,tol);
if( nw==0 )
yr(i) = csr;
yi(i) = csi;
nz = fix(nz - 1);
if( ic==kk-1 )
nz = fix(kk - 2);
igo=1;
break;
else;
ic = fix(kk);
continue;
end;
end;
end;
if( alas>=helim )
zdr = zdr - elim;
s1r = s1r.*celmr;
s1i = s1i.*celmr;
s2r = s2r.*celmr;
s2i = s2i.*celmr;
end;
end;
if(igo==0)
nz = fix(n);
if( ic==n )
nz = fix(n - 1);
end;
end;
for i = 1 : nz;
yr(i) = zeror;
yi(i) = zeroi;
end; i = fix(nz+1);
end %subroutine zkscl
%DECK ZLOG

Contact us at files@mathworks.com