Code covered by the BSD License  

Highlights from
slatec

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

[zr,fnu,n,y,nz,rz,ascle,tol,elim]=ckscl(zr,fnu,n,y,nz,rz,ascle,tol,elim);
function [zr,fnu,n,y,nz,rz,ascle,tol,elim]=ckscl(zr,fnu,n,y,nz,rz,ascle,tol,elim);
%***BEGIN PROLOGUE  CKSCL
%***SUBSIDIARY
%***PURPOSE  Subsidiary to CBKNU, CUNK1 and CUNK2
%***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  CBKNU, CUNK1, CUNK2
%***ROUTINES CALLED  CUCHK
%***REVISION HISTORY  (YYMMDD)
%   ??????  DATE WRITTEN
%   910415  Prologue converted to Version 4.0 format.  (BAB)
%***end PROLOGUE  CKSCL
persistent aa acs alas as celm ck cs csi csr cy czero elm firstCall fn helim i ic k kk nn nw s1 s2 xx zd zri ; if isempty(firstCall),firstCall=1;end; 

if isempty(ck), ck=0; end;
if isempty(cs), cs=0; end;
if isempty(cy), cy=zeros(1,2); end;
if isempty(czero), czero=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(zd), zd=0; end;
if isempty(celm), celm=0; end;
if isempty(aa), aa=0; end;
if isempty(acs), acs=0; end;
if isempty(as), as=0; end;
if isempty(csi), csi=0; end;
if isempty(csr), csr=0; end;
if isempty(fn), fn=0; end;
if isempty(xx), xx=0; end;
if isempty(zri), zri=0; end;
if isempty(elm), elm=0; end;
if isempty(alas), alas=0; end;
if isempty(helim), helim=0; end;
if isempty(i), i=0; end;
if isempty(ic), ic=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(nn), nn=0; end;
if isempty(nw), nw=0; end;
if firstCall,   czero=[complex(0.0e0,0.0e0)];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  CUCHK
nz = 0;
ic = 0;
xx = real(zr);
nn = fix(min(2,n));
for i = 1 : nn;
s1 = y(i);
cy(i) = s1;
as = abs(s1);
acs = -xx + log(as);
nz = fix(nz + 1);
y(i) = czero;
if( acs>=(-elim) )
cs = -zr + log(s1);
csr = real(cs);
csi = imag(cs);
aa = exp(csr)./tol;
cs = complex(aa,0.0e0).*complex(cos(csi),sin(csi));
[cs,nw,ascle,tol]=cuchk(cs,nw,ascle,tol);
if( nw==0 )
y(i) = cs;
nz = fix(nz - 1);
ic = fix(i);
end;
end;
end; i = fix(nn+1);
if( n==1 )
return;
end;
if( ic<=1 )
y(1) = czero;
nz = 2;
end;
if( n==2 )
return;
end;
if( nz==0 )
return;
end;
fn = fnu + 1.0e0;
ck = complex(fn,0.0e0).*rz;
s1 = cy(1);
s2 = cy(2);
helim = 0.5e0.*elim;
elm = exp(-elim);
celm = complex(elm,0.0e0);
zri = imag(zr);
zd = zr;
%
%     FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
%     S2 GETS LARGER THAN EXP(ELIM/2)
%
for i = 3 : n;
kk = fix(i);
cs = s2;
s2 = ck.*s2 + s1;
s1 = cs;
ck = ck + rz;
as = abs(s2);
alas = log(as);
acs = -xx + alas;
nz = fix(nz + 1);
y(i) = czero;
if( acs>=(-elim) )
cs = -zd + log(s2);
csr = real(cs);
csi = imag(cs);
aa = exp(csr)./tol;
cs = complex(aa,0.0e0).*complex(cos(csi),sin(csi));
[cs,nw,ascle,tol]=cuchk(cs,nw,ascle,tol);
if( nw==0 )
y(i) = cs;
nz = fix(nz - 1);
if( ic==(kk-1) )
nz = fix(kk - 2);
for k = 1 : nz;
y(k) = czero;
end; k = fix(nz+1);
return;
else;
ic = fix(kk);
continue;
end;
end;
end;
if( alas>=helim )
xx = xx - elim;
s1 = s1.*celm;
s2 = s2.*celm;
zd = complex(xx,zri);
end;
end; i = fix(n+1);
nz = fix(n);
if( ic==n )
nz = fix(n - 1);
end;
for k = 1 : nz;
y(k) = czero;
end; k = fix(nz+1);
end
%DECK CLBETA

Contact us at files@mathworks.com